Line data Source code
1 : !=======================================================================
2 :
3 : ! Water isotope tracers within sea ice and snow
4 : !
5 : ! authors: David Bailey, NCAR
6 : ! Jiang Zhu, UW-Madison
7 : !
8 : ! 2014: Added i2x evaporation flux
9 : ! Added fractionation options
10 : ! Fixed bugs
11 : !
12 : module icepack_isotope
13 :
14 : use icepack_kinds
15 : use icepack_parameters, only: c0, c1, c2, p001, p1, p5, puny
16 : use icepack_tracers, only: n_iso
17 : use icepack_warnings, only: icepack_warnings_add, icepack_warnings_setabort
18 :
19 : implicit none
20 : private
21 :
22 : public :: update_isotope, &
23 : isoice_alpha
24 :
25 : character(len=5), parameter, public :: &
26 : isotope_frac_method = 'cfrac' ! fractionation coefficient calculation method
27 : ! cfrac, constant fractionation
28 : ! nfrac, nonfractionation
29 : ! gfrac, growth-rate dependent for H2_18O
30 :
31 : ! Species indicies - public so thay can be seen by water_tracers
32 : integer, parameter, public :: ispundef = 0 ! Undefined
33 : integer, parameter, public :: isph2o = 1 ! H2O "regular" water
34 : integer, parameter, public :: isph216o = 2 ! H216O similar to "regular"
35 : integer, parameter, public :: isphdo = 3 ! HDO
36 : integer, parameter, public :: isph218o = 4 ! H218O
37 : integer, parameter, public :: pwtspec = 4 ! number of water species
38 : ! h2o,hdo,h218o,h216o
39 :
40 : !=======================================================================
41 :
42 : contains
43 :
44 : !=======================================================================
45 : !
46 : ! Increase isotope in ice or snow surface due to deposition and loss
47 : !
48 21948794 : subroutine update_isotope (dt, &
49 : nilyr, nslyr, &
50 : meltt, melts, &
51 : meltb, congel, &
52 : snoice, evap, &
53 : fsnow, Tsfc, &
54 21948794 : Qref_iso, &
55 21948794 : isosno, isoice, &
56 : aice_old, &
57 : vice_old, vsno_old, &
58 : vicen, vsnon, aicen, &
59 43897588 : fiso_atm, fiso_evapn, &
60 21948794 : fiso_ocnn, HDO_ocn, &
61 : H2_16O_ocn, H2_18O_ocn)
62 :
63 : ! use water_isotopes, only: wiso_alpi
64 : use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh
65 :
66 : integer (kind=int_kind), intent(in) :: &
67 : nilyr, nslyr
68 :
69 : real (kind=dbl_kind), intent(in) :: &
70 : dt ! time step
71 :
72 : real (kind=dbl_kind), intent(in) :: &
73 : Tsfc, & ! surface temperature
74 : meltt, & ! top melt
75 : melts, & ! snow melt
76 : meltb, & ! bottom melt
77 : congel, & ! congelation ice growth (m/step)
78 : snoice, & ! ice thickness increase (m/step)
79 : evap, & ! surface evaporation
80 : fsnow, & ! snowfall (kg/m^2/s of water)
81 : vicen, & ! volume per unit area of ice (m)
82 : vsnon, & ! volume per unit area of snow (m)
83 : aicen, & ! ice area
84 : aice_old, & ! beginning values
85 : vice_old, &
86 : vsno_old, &
87 : HDO_ocn, & !
88 : H2_16O_ocn, & !
89 : H2_18O_ocn !
90 :
91 : real (kind=dbl_kind), dimension(:), intent(in) :: &
92 : fiso_atm, & ! isotopic snowfall (kg/m^2/s of water)
93 : Qref_iso ! isotope reference humidity
94 :
95 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
96 : fiso_ocnn, & ! isotopic freshwater (kg/m^2/s)
97 : fiso_evapn ! evaporative water flux (kg/m^2/s)
98 :
99 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
100 : isosno, isoice ! mass of isotopes (kg)
101 :
102 : ! local variables
103 :
104 : real (kind=dbl_kind) :: &
105 205470 : dzsno, &
106 205470 : dzice, &
107 205470 : evaps, & ! evaporation over snow (m/step)
108 205470 : evapi, & ! evaporation over ice (m/step)
109 205470 : dhs_snoice,& ! snow thickness reduction (m/step)
110 205470 : hi, &
111 205470 : hs
112 :
113 : real (kind=dbl_kind), dimension(n_iso) :: &
114 : isotot, isotot0 ! for diagnostics
115 :
116 : real (kind=dbl_kind) :: &
117 616410 : hs_old, hi_old, dhs, dhi, sloss1, &
118 205470 : TsfK, & ! snow/ice temperature (K)
119 205470 : alphai, & ! ice/vapour fractionation coefficient
120 205470 : ratio, & ! isotopic ratio
121 205470 : work, & ! temporary variable
122 205470 : alpha
123 :
124 : integer (kind=int_kind) :: k
125 :
126 : character(len=*),parameter :: subname='(update_isotope)'
127 :
128 : ! initialize
129 :
130 21948794 : hs_old=vsno_old/aice_old
131 21948794 : hi_old=vice_old/aice_old
132 :
133 21948794 : dzsno = hs_old
134 21948794 : dzice = hi_old
135 :
136 21948794 : if (aicen > puny) then
137 21945663 : hs = vsnon/aicen
138 21945663 : hi = vicen/aicen
139 3131 : elseif (aice_old > puny) then
140 3131 : hs = vsnon/aice_old
141 3131 : hi = vicen/aice_old
142 : endif
143 :
144 21948794 : if (ktherm == 2) then
145 21948794 : dhs_snoice = snoice
146 : else
147 0 : dhs_snoice = snoice*rhoi/rhos
148 : endif
149 :
150 : ! if (hs > puny) then
151 : ! evaps = evap*dt/rhos
152 : ! else
153 : ! evapi = evap*dt/rhoi
154 : ! endif
155 21948794 : evaps = hs - (hs_old - melts - dhs_snoice + fsnow/rhos*dt)
156 21948794 : evapi = hi - (hi_old - meltt - meltb + congel + snoice)
157 :
158 : ! condensation of vapor onto snow and ice
159 :
160 21948794 : TsfK = Tsfc + Tffresh
161 :
162 21948794 : if (evaps > c0) then ! condensation to snow
163 33520540 : do k = 1, n_iso
164 25140405 : ratio = c1 ! ratio between 18O(HDO) and 16O in humidity
165 25140405 : alphai = c1 ! fractionation coefficient
166 25140405 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
167 0 : ratio = Qref_iso(k)/Qref_iso(2)
168 25140405 : if (Qref_iso(2) <= puny) ratio = c0
169 25140405 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
170 25140405 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
171 25140405 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
172 25140405 : work = alphai*ratio*rhos*evaps*aicen
173 25140405 : fiso_evapn(k) = fiso_evapn(k) + work/dt
174 33520540 : isosno(k) = isosno(k) + work
175 : enddo
176 8380135 : dzsno = dzsno + evaps
177 : endif
178 :
179 21948794 : if (evapi > c0) then ! condensation to ice
180 24545892 : do k = 1, n_iso
181 18409419 : ratio = c1 ! ratio between 18O(HDO) and 16O in ref humidity
182 18409419 : alphai = c1 ! fractionation coefficient
183 18409419 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
184 0 : ratio = Qref_iso(k)/Qref_iso(2)
185 18409419 : if (Qref_iso(2) <= puny) ratio = c0
186 18409419 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
187 18409419 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
188 18409419 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
189 18409419 : work = alphai*ratio*rhoi*evapi*aicen
190 18409419 : fiso_evapn(k) = fiso_evapn(k) + work/dt
191 24545892 : isoice(k) = isoice(k) + work
192 : enddo
193 6136473 : dzice = dzice + evapi
194 : endif
195 :
196 : ! basal ice growth and isotope uptake
197 :
198 21948794 : if (congel > c0) then
199 22384532 : do k = 1,n_iso
200 16788399 : if (k == 1) then
201 5596133 : alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method)
202 5596133 : work = alpha*HDO_ocn*rhoi*congel*aicen
203 11192266 : elseif (k == 2) then
204 5596133 : alpha = isoice_alpha(congel/dt,'H2_16O',isotope_frac_method)
205 5596133 : work = alpha*H2_16O_ocn*rhoi*congel*aicen
206 5596133 : elseif (k == 3) then
207 5596133 : alpha = isoice_alpha(congel/dt,'H2_18O',isotope_frac_method)
208 5596133 : work = alpha*H2_18O_ocn*rhoi*congel*aicen
209 : else
210 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
211 0 : call icepack_warnings_add(subname//' ERROR: n_iso > 3')
212 0 : return
213 : endif
214 16788399 : isoice(k) = isoice(k) + work
215 22384532 : fiso_ocnn(k) = fiso_ocnn(k) - work/dt
216 : enddo
217 :
218 5596133 : dzice = dzice + congel
219 : endif
220 :
221 : ! sublimation of snow and ice
222 :
223 21948794 : if (evaps < c0) then ! snow sublimation (no fractionation)
224 49758836 : do k = 1, n_iso
225 : !ratio = c1 ! ratio between 18O(HDO) and 16O in snow
226 : !if (isosno(2) > puny) ratio = isosno(k)/isosno(2)
227 : !if (ratio > c5) ratio = c1 !! remove latter?
228 : !work = ratio*rhos*evaps*aicen
229 : !fiso_evapn(k) = fiso_evapn(k)+work/dt
230 :
231 37319127 : sloss1 = c0
232 37319127 : if (dzsno > puny) sloss1 = isosno(k)*min(-evaps,dzsno)/dzsno
233 37319127 : if (isosno(k) >= sloss1) then
234 37304160 : isosno(k) = isosno(k)-sloss1
235 : else
236 14967 : sloss1 = isosno(k)
237 14967 : isosno(k) = c0
238 : endif
239 : ! if (isosno(k) < c0) then
240 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
241 : ! endif
242 49758836 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
243 : enddo
244 :
245 12439709 : dzsno = dzsno + evaps
246 12439709 : if (dzsno <= c0) then ! snow goes away
247 454856 : fiso_evapn(:) = fiso_evapn(:) - isosno(:)/dt
248 454856 : isosno(:) = c0
249 113714 : dzsno = c0
250 : endif
251 : endif
252 :
253 21948794 : if (evapi < c0) then ! ice sublimation (no fractionation)
254 21769128 : do k = 1, n_iso
255 : !!ratio = c1 ! ratio between 18O(HDO) and 16O in ice
256 : !!if (isoice(2) > puny) ratio = isoice(k)/isoice(2)
257 : !!if (ratio > c5) ratio = c1 ! remove later?
258 : !!work = ratio*rhoi*evapi*aicen
259 : !!fiso_evapn(k) = fiso_evapn(k)+work/dt
260 :
261 16326846 : sloss1 = c0
262 16326846 : if (dzice > puny) &
263 16326846 : sloss1 = isoice(k) * min(-evapi,dzice)/dzice
264 16326846 : if (isoice(k) >= sloss1) then
265 16326846 : isoice(k) = isoice(k)-sloss1
266 : else
267 0 : sloss1 = isoice(k)
268 0 : isoice(k) = c0
269 : endif
270 21769128 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
271 : enddo
272 :
273 5442282 : dzice = dzice + evapi
274 5442282 : if (dzice <= c0) then ! ice goes away
275 0 : fiso_evapn(:) = fiso_evapn(:) - isoice(:)/dt
276 0 : isoice(:) = c0
277 0 : dzice = c0
278 : endif
279 : endif
280 :
281 : ! surface snow melt
282 :
283 21948794 : if (melts > c0) then
284 32290704 : do k=1,n_iso
285 24218028 : sloss1=c0
286 24218028 : if (dzsno > puny) &
287 24151944 : sloss1 = isosno(k) * min(melts,dzsno)/dzsno
288 24218028 : if (isosno(k) >= sloss1) then
289 24135906 : isosno(k) = isosno(k)-sloss1
290 : else
291 82122 : sloss1 = isosno(k)
292 82122 : isosno(k) = c0
293 : endif
294 : ! if (isosno(k) < c0) then
295 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
296 : ! endif
297 32290704 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
298 : enddo ! n_iso
299 :
300 8072676 : dzsno = dzsno - melts
301 8072676 : if (dzsno <= c0) then ! snow melts away
302 2364164 : fiso_ocnn(:) = fiso_ocnn(:) + isosno(:)/dt
303 2364164 : isosno(:) = c0
304 591041 : dzsno = c0
305 : endif
306 : endif
307 :
308 : ! surface ice melt
309 21948794 : if (meltt > c0) then
310 8172568 : do k=1,n_iso
311 6129426 : sloss1=c0
312 6129426 : if (dzice > puny) sloss1=isoice(k) * min(meltt,dzice)/dzice
313 6129426 : if (isoice(k) >= sloss1) then
314 6129396 : isoice(k) = isoice(k)-sloss1
315 : else
316 30 : sloss1 = isoice(k)
317 30 : isoice(k) = c0
318 : endif
319 8172568 : fiso_ocnn(k)=fiso_ocnn(k) + sloss1/dt
320 : enddo
321 :
322 2043142 : dzice = dzice - meltt
323 2043142 : if (dzice <= c0) then ! ice ice melts away
324 1192 : fiso_ocnn(:) = fiso_ocnn(:)+isoice(:)
325 1192 : isoice(:) = c0
326 298 : dzice = c0
327 : endif
328 : endif
329 :
330 : ! basal ice melt. Assume all isotopes lost in basal melt
331 :
332 21948794 : if (meltb > c0) then
333 65409456 : do k=1,n_iso
334 49057092 : sloss1=c0
335 49057092 : if (dzice > puny) sloss1=max(meltb-dzice,c0) * isoice(k)/dzice
336 49057092 : if (isoice(k) >= sloss1) then
337 49056912 : isoice(k) = isoice(k)-sloss1
338 : else
339 180 : sloss1 = isoice(k)
340 180 : isoice(k) = c0
341 : endif
342 65409456 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
343 : enddo
344 :
345 16352364 : dzice = dzice - meltb
346 16352364 : if (dzice <= c0) then ! ice ice melts away
347 20416 : fiso_ocnn(:) = fiso_ocnn(:) + isoice(:)
348 20416 : isoice(:) = c0
349 5104 : dzice = c0
350 : endif
351 : endif
352 :
353 : ! snowfall and isotope deposition
354 :
355 21948794 : if (fsnow > c0) then
356 49569768 : isosno(:) = isosno(:) + fiso_atm(:)*aicen*dt
357 12392442 : dzsno = dzsno + fsnow/rhos*dt
358 : endif
359 :
360 : ! snoice formation
361 :
362 21948794 : if (dhs_snoice > c0) then
363 4753332 : do k=1,n_iso
364 3564999 : sloss1=c0
365 3564999 : if (dzsno > puny) sloss1 = min(dhs_snoice,dzsno) * isosno(k)/dzsno
366 3564999 : if (isosno(k) >= sloss1) then
367 3564810 : isosno(k) = isosno(k)-sloss1
368 : else
369 189 : sloss1 = isosno(k)
370 189 : isosno(k) = c0
371 : endif
372 : ! if (isosno(k) < c0) then
373 : ! write(nu_diag,*) 'Snow-ice isosno(k)',isosno(k),sloss1
374 : ! endif
375 4753332 : isoice(k) = isoice(k) + sloss1
376 : enddo
377 :
378 1188333 : dzsno = dzsno - dhs_snoice
379 1188333 : dzice = dzice + snoice
380 1188333 : if (dzsno <= c0) then ! snow goes away
381 6180 : fiso_ocnn(:)= fiso_ocnn(:) + isosno(:)/dt
382 6180 : isosno(:) = c0
383 1545 : dzsno = c0
384 : endif
385 : endif
386 :
387 : ! do k=1,n_iso
388 : ! isotot(k) = isosno(k) + isoice(k)
389 :
390 : ! if ( (isotot(k)-isotot0(k)) &
391 : ! - fiso_atm (k)*dt*aicen &
392 : ! - fiso_evapn(k)*dt &
393 : ! + fiso_ocnn (k)*dt > 1e-6) then
394 : ! write(nu_diag,*) 'isotope tracer: ',k
395 : ! write(nu_diag,*) 'isotot-isotot0 ',isotot(k)-isotot0(k)
396 : ! write(nu_diag,*) 'fiso_atm-fiso_ocnn ',fiso_atm (k)*dt*aicen &
397 : ! + fiso_evapn(k)*dt &
398 : ! - fiso_ocnn (k)*dt
399 : ! endif
400 : ! enddo ! n_iso
401 :
402 : ! scale fiso_ocnn. It will be re-scaled by aicen later in merge_fluxes
403 21948794 : if (aicen > puny) then
404 87782652 : fiso_ocnn(:) = fiso_ocnn(:)/aicen
405 87782652 : fiso_evapn(:) = fiso_evapn(:)/aicen
406 : endif
407 :
408 : end subroutine update_isotope
409 :
410 : !=======================================================================
411 :
412 : ! calculate the fractionation coefficient for sea-ice formation
413 :
414 41939484 : function isoice_alpha(growth_rate, sp, frac)
415 : !
416 : ! authors: Jiang Zhu, UW-Madison
417 : !
418 : real (kind=dbl_kind), intent(in) :: &
419 : growth_rate ! sea-ice formation rate (m/s)
420 : character(*), intent(in) :: &
421 : sp,frac ! species: H2_16O, H2_18O, HDO
422 : ! calculation methods:
423 : ! cfrac, constant fractionation
424 : ! nfrac, nonfractionation
425 : ! gfrac, growth-rate dependent
426 : real (kind=dbl_kind) :: &
427 : isoice_alpha ! return fractionation
428 :
429 : character(len=*),parameter :: subname='(isoice_alpha)'
430 :
431 0 : if (frac == 'nfrac') isoice_alpha = c1
432 41939484 : if (sp == 'H2_16O') isoice_alpha = c1
433 :
434 : ! Lehmann and Siegenthaler, 1991
435 : !--------------------------------------------------
436 41939484 : if (frac == 'cfrac' .and. sp == 'HDO') &
437 13979828 : isoice_alpha = 1.02120_dbl_kind
438 41939484 : if (frac == 'cfrac' .and. sp == 'H2_18O') &
439 13979828 : isoice_alpha = 1.00291_dbl_kind
440 :
441 : ! Eq.9, Toyota et al., 2013
442 : ! For HDO, 7.2852 = 0.2120/0.00291
443 : !--------------------------------------------------
444 41939484 : if (frac == 'gfrac' .and. sp == 'HDO') &
445 : isoice_alpha = c1+7.2852_dbl_kind*1.2280E-3_dbl_kind+ &
446 : 0.7311E-3_dbl_kind*exp(-growth_rate/8.0100E8_dbl_kind)+ &
447 0 : 0.8441E-3_dbl_kind*exp(-growth_rate/0.7800E6_dbl_kind)
448 41939484 : if (frac == 'gfrac' .and. sp == 'H2_18O') &
449 : isoice_alpha = c1+1.2280E-3_dbl_kind+ &
450 : 0.7311E-3_dbl_kind*exp(-growth_rate/8.0100E8_dbl_kind)+ &
451 0 : 0.8441E-3_dbl_kind*exp(-growth_rate/0.7800E6_dbl_kind)
452 41939484 : return
453 :
454 41939484 : end function isoice_alpha
455 :
456 : !=======================================================================
457 :
458 43549824 : function wiso_alpi(isp,tk)
459 :
460 : !-----------------------------------------------------------------------
461 : ! Purpose: return ice/vapour fractionation from loop-up tables
462 : ! Author: David Noone <dcn@caltech.edu> - Tue Jul 1 12:02:24 MDT 2003
463 : !-----------------------------------------------------------------------
464 :
465 : integer , intent(in) :: isp ! species indes
466 : real(kind=dbl_kind), intent(in) :: tk ! temperature (k)
467 : real(kind=dbl_kind) :: wiso_alpi ! return fractionation
468 :
469 : character(len=*),parameter :: subname='(wiso_alpi)'
470 :
471 : !From Merlivat & Nief,1967 for HDO, and Majoube, 1971b for H218O:
472 : real(kind=dbl_kind), parameter, dimension(pwtspec) :: & ! ice/vapour
473 : alpai = (/ 0._dbl_kind, 0._dbl_kind, 16289._dbl_kind, 0._dbl_kind /), &
474 : alpbi = (/ 0._dbl_kind, 0._dbl_kind, 0._dbl_kind, 11.839_dbl_kind /), &
475 : alpci = (/ 0._dbl_kind, 0._dbl_kind, -9.45e-2_dbl_kind, -28.224e-3_dbl_kind /)
476 :
477 : !-----------------------------------------------------------------------
478 43549824 : if (isp == isph2o) then
479 0 : wiso_alpi = c1
480 0 : return
481 : end if
482 :
483 43549824 : wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp))
484 :
485 43549824 : return
486 : end function wiso_alpi
487 :
488 : !=======================================================================
489 :
490 : end module icepack_isotope
491 :
492 : !=======================================================================
|
