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 229760 : subroutine update_isotope (dt, &
49 : nilyr, nslyr, &
50 : meltt, melts, &
51 : meltb, congel, &
52 : snoice, evap, &
53 : fsnow, Tsfc, &
54 229760 : Qref_iso, &
55 229760 : isosno, isoice, &
56 : aice_old, &
57 : vice_old, vsno_old, &
58 : vicen, vsnon, aicen, &
59 459520 : fiso_atm, fiso_evapn, &
60 229760 : 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 0 : dzsno, &
106 0 : dzice, &
107 0 : evaps, & ! evaporation over snow (m/step)
108 0 : evapi, & ! evaporation over ice (m/step)
109 0 : dhs_snoice,& ! snow thickness reduction (m/step)
110 0 : hi, &
111 0 : hs
112 :
113 : real (kind=dbl_kind), dimension(n_iso) :: &
114 : isotot, isotot0 ! for diagnostics
115 :
116 : real (kind=dbl_kind) :: &
117 0 : hs_old, hi_old, dhs, dhi, sloss1, &
118 0 : TsfK, & ! snow/ice temperature (K)
119 0 : alphai, & ! ice/vapour fractionation coefficient
120 0 : ratio, & ! isotopic ratio
121 0 : work, & ! temporary variable
122 0 : alpha
123 :
124 : integer (kind=int_kind) :: k
125 :
126 : character(len=*),parameter :: subname='(update_isotope)'
127 :
128 : ! initialize
129 :
130 229760 : hs_old=vsno_old/aice_old
131 229760 : hi_old=vice_old/aice_old
132 :
133 229760 : dzsno = hs_old
134 229760 : dzice = hi_old
135 :
136 229760 : if (aicen > puny) then
137 229759 : hs = vsnon/aicen
138 229759 : hi = vicen/aicen
139 1 : elseif (aice_old > puny) then
140 1 : hs = vsnon/aice_old
141 1 : hi = vicen/aice_old
142 : endif
143 :
144 229760 : if (ktherm == 2) then
145 229760 : 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 229760 : evaps = hs - (hs_old - melts - dhs_snoice + fsnow/rhos*dt)
156 229760 : evapi = hi - (hi_old - meltt - meltb + congel + snoice)
157 :
158 : ! condensation of vapor onto snow and ice
159 :
160 229760 : TsfK = Tsfc + Tffresh
161 :
162 229760 : if (evaps > c0) then ! condensation to snow
163 2208 : do k = 1, n_iso
164 1656 : ratio = c1 ! ratio between 18O(HDO) and 16O in humidity
165 1656 : alphai = c1 ! fractionation coefficient
166 1656 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
167 0 : ratio = Qref_iso(k)/Qref_iso(2)
168 1656 : if (Qref_iso(2) <= puny) ratio = c0
169 1656 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
170 1656 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
171 1656 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
172 1656 : work = alphai*ratio*rhos*evaps*aicen
173 1656 : fiso_evapn(k) = fiso_evapn(k) + work/dt
174 2208 : isosno(k) = isosno(k) + work
175 : enddo
176 552 : dzsno = dzsno + evaps
177 : endif
178 :
179 229760 : if (evapi > c0) then ! condensation to ice
180 180120 : do k = 1, n_iso
181 135090 : ratio = c1 ! ratio between 18O(HDO) and 16O in ref humidity
182 135090 : alphai = c1 ! fractionation coefficient
183 135090 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
184 0 : ratio = Qref_iso(k)/Qref_iso(2)
185 135090 : if (Qref_iso(2) <= puny) ratio = c0
186 135090 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
187 135090 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
188 135090 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
189 135090 : work = alphai*ratio*rhoi*evapi*aicen
190 135090 : fiso_evapn(k) = fiso_evapn(k) + work/dt
191 180120 : isoice(k) = isoice(k) + work
192 : enddo
193 45030 : dzice = dzice + evapi
194 : endif
195 :
196 : ! basal ice growth and isotope uptake
197 :
198 229760 : if (congel > c0) then
199 491948 : do k = 1,n_iso
200 368961 : if (k == 1) then
201 122987 : alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method)
202 122987 : work = alpha*HDO_ocn*rhoi*congel*aicen
203 245974 : elseif (k == 2) then
204 122987 : alpha = isoice_alpha(congel/dt,'H2_16O',isotope_frac_method)
205 122987 : work = alpha*H2_16O_ocn*rhoi*congel*aicen
206 122987 : elseif (k == 3) then
207 122987 : alpha = isoice_alpha(congel/dt,'H2_18O',isotope_frac_method)
208 122987 : 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 368961 : isoice(k) = isoice(k) + work
215 491948 : fiso_ocnn(k) = fiso_ocnn(k) - work/dt
216 : enddo
217 :
218 122987 : dzice = dzice + congel
219 : endif
220 :
221 : ! sublimation of snow and ice
222 :
223 229760 : if (evaps < c0) then ! snow sublimation (no fractionation)
224 906952 : 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 680214 : sloss1 = c0
232 680214 : if (dzsno > puny) sloss1 = isosno(k)*min(-evaps,dzsno)/dzsno
233 680214 : if (isosno(k) >= sloss1) then
234 674783 : isosno(k) = isosno(k)-sloss1
235 : else
236 5431 : sloss1 = isosno(k)
237 5431 : isosno(k) = c0
238 : endif
239 : ! if (isosno(k) < c0) then
240 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
241 : ! endif
242 906952 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
243 : enddo
244 :
245 226738 : dzsno = dzsno + evaps
246 226738 : if (dzsno <= c0) then ! snow goes away
247 115300 : fiso_evapn(:) = fiso_evapn(:) - isosno(:)/dt
248 115300 : isosno(:) = c0
249 28825 : dzsno = c0
250 : endif
251 : endif
252 :
253 229760 : if (evapi < c0) then ! ice sublimation (no fractionation)
254 347308 : 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 260481 : sloss1 = c0
262 260481 : if (dzice > puny) &
263 260481 : sloss1 = isoice(k) * min(-evapi,dzice)/dzice
264 260481 : if (isoice(k) >= sloss1) then
265 260481 : isoice(k) = isoice(k)-sloss1
266 : else
267 0 : sloss1 = isoice(k)
268 0 : isoice(k) = c0
269 : endif
270 347308 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
271 : enddo
272 :
273 86827 : dzice = dzice + evapi
274 86827 : 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 229760 : if (melts > c0) then
284 299520 : do k=1,n_iso
285 224640 : sloss1=c0
286 224640 : if (dzsno > puny) &
287 224640 : sloss1 = isosno(k) * min(melts,dzsno)/dzsno
288 224640 : if (isosno(k) >= sloss1) then
289 224639 : isosno(k) = isosno(k)-sloss1
290 : else
291 1 : sloss1 = isosno(k)
292 1 : isosno(k) = c0
293 : endif
294 : ! if (isosno(k) < c0) then
295 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
296 : ! endif
297 299520 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
298 : enddo ! n_iso
299 :
300 74880 : dzsno = dzsno - melts
301 74880 : if (dzsno <= c0) then ! snow melts away
302 60 : fiso_ocnn(:) = fiso_ocnn(:) + isosno(:)/dt
303 60 : isosno(:) = c0
304 15 : dzsno = c0
305 : endif
306 : endif
307 :
308 : ! surface ice melt
309 229760 : if (meltt > c0) then
310 80648 : do k=1,n_iso
311 60486 : sloss1=c0
312 60486 : if (dzice > puny) sloss1=isoice(k) * min(meltt,dzice)/dzice
313 60486 : if (isoice(k) >= sloss1) then
314 60486 : isoice(k) = isoice(k)-sloss1
315 : else
316 0 : sloss1 = isoice(k)
317 0 : isoice(k) = c0
318 : endif
319 80648 : fiso_ocnn(k)=fiso_ocnn(k) + sloss1/dt
320 : enddo
321 :
322 20162 : dzice = dzice - meltt
323 20162 : if (dzice <= c0) then ! ice ice melts away
324 0 : fiso_ocnn(:) = fiso_ocnn(:)+isoice(:)
325 0 : isoice(:) = c0
326 0 : dzice = c0
327 : endif
328 : endif
329 :
330 : ! basal ice melt. Assume all isotopes lost in basal melt
331 :
332 229760 : if (meltb > c0) then
333 427092 : do k=1,n_iso
334 320319 : sloss1=c0
335 320319 : if (dzice > puny) sloss1=max(meltb-dzice,c0) * isoice(k)/dzice
336 320319 : if (isoice(k) >= sloss1) then
337 320319 : isoice(k) = isoice(k)-sloss1
338 : else
339 0 : sloss1 = isoice(k)
340 0 : isoice(k) = c0
341 : endif
342 427092 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
343 : enddo
344 :
345 106773 : dzice = dzice - meltb
346 106773 : if (dzice <= c0) then ! ice ice melts away
347 4 : fiso_ocnn(:) = fiso_ocnn(:) + isoice(:)
348 4 : isoice(:) = c0
349 1 : dzice = c0
350 : endif
351 : endif
352 :
353 : ! snowfall and isotope deposition
354 :
355 229760 : if (fsnow > c0) then
356 876956 : isosno(:) = isosno(:) + fiso_atm(:)*aicen*dt
357 219239 : dzsno = dzsno + fsnow/rhos*dt
358 : endif
359 :
360 : ! snoice formation
361 :
362 229760 : if (dhs_snoice > c0) then
363 12536 : do k=1,n_iso
364 9402 : sloss1=c0
365 9402 : if (dzsno > puny) sloss1 = min(dhs_snoice,dzsno) * isosno(k)/dzsno
366 9402 : if (isosno(k) >= sloss1) then
367 9402 : isosno(k) = isosno(k)-sloss1
368 : else
369 0 : sloss1 = isosno(k)
370 0 : 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 12536 : isoice(k) = isoice(k) + sloss1
376 : enddo
377 :
378 3134 : dzsno = dzsno - dhs_snoice
379 3134 : dzice = dzice + snoice
380 3134 : if (dzsno <= c0) then ! snow goes away
381 0 : fiso_ocnn(:)= fiso_ocnn(:) + isosno(:)/dt
382 0 : isosno(:) = c0
383 0 : 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 229760 : if (aicen > puny) then
404 919036 : fiso_ocnn(:) = fiso_ocnn(:)/aicen
405 919036 : 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 462654 : 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 462654 : if (sp == 'H2_16O') isoice_alpha = c1
433 :
434 : ! Lehmann and Siegenthaler, 1991
435 : !--------------------------------------------------
436 462654 : if (frac == 'cfrac' .and. sp == 'HDO') &
437 154218 : isoice_alpha = 1.02120_dbl_kind
438 462654 : if (frac == 'cfrac' .and. sp == 'H2_18O') &
439 154218 : 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 462654 : 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 462654 : 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 462654 : return
453 :
454 462654 : end function isoice_alpha
455 :
456 : !=======================================================================
457 :
458 136746 : 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 136746 : if (isp == isph2o) then
479 0 : wiso_alpi = c1
480 0 : return
481 : end if
482 :
483 136746 : wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp))
484 :
485 : #ifdef NOFRAC
486 : wiso_alpi = c1
487 : #endif
488 : !
489 136746 : return
490 : end function wiso_alpi
491 :
492 : !=======================================================================
493 :
494 : end module icepack_isotope
495 :
496 : !=======================================================================
|
