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 225452 : subroutine update_isotope(dt, &
49 : meltt, melts, &
50 : meltb, congel, &
51 : snoice, evap, &
52 : fsnow, Tsfc, &
53 225452 : Qref_iso, &
54 225452 : isosno, isoice, &
55 : aice_old, &
56 : vice_old, vsno_old, &
57 : vicen, vsnon, aicen, &
58 450904 : fiso_atm, fiso_evapn, &
59 225452 : fiso_ocnn, HDO_ocn, &
60 : H2_16O_ocn, H2_18O_ocn)
61 :
62 : ! use water_isotopes, only: wiso_alpi
63 : use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh
64 :
65 : real (kind=dbl_kind), intent(in) :: &
66 : dt ! time step
67 :
68 : real (kind=dbl_kind), intent(in) :: &
69 : Tsfc, & ! surface temperature
70 : meltt, & ! top melt
71 : melts, & ! snow melt
72 : meltb, & ! bottom melt
73 : congel, & ! congelation ice growth (m/step)
74 : snoice, & ! ice thickness increase (m/step)
75 : evap, & ! surface evaporation
76 : fsnow, & ! snowfall (kg/m^2/s of water)
77 : vicen, & ! volume per unit area of ice (m)
78 : vsnon, & ! volume per unit area of snow (m)
79 : aicen, & ! ice area
80 : aice_old, & ! beginning values
81 : vice_old, &
82 : vsno_old
83 :
84 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
85 : fiso_ocnn, & ! isotopic freshwater (kg/m^2/s)
86 : fiso_evapn ! evaporative water flux (kg/m^2/s)
87 :
88 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
89 : isosno, & ! mass of isotopes (kg)
90 : isoice
91 :
92 : real (kind=dbl_kind), dimension(:), intent(in) :: &
93 : fiso_atm, & ! isotopic snowfall (kg/m^2/s of water)
94 : Qref_iso ! isotope reference humidity
95 :
96 : real (kind=dbl_kind), intent(in) :: &
97 : HDO_ocn, & ! ocean concentration of HDO (kg/kg)
98 : H2_16O_ocn, & ! ocean concentration of H2_16O (kg/kg)
99 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
100 :
101 : ! local variables
102 :
103 : real (kind=dbl_kind) :: &
104 : dzsno, &
105 : dzice, &
106 : evaps, & ! evaporation over snow (m/step)
107 : evapi, & ! evaporation over ice (m/step)
108 : dhs_snoice,& ! snow thickness reduction (m/step)
109 : hi, &
110 : hs
111 :
112 : ! real (kind=dbl_kind), dimension(n_iso) :: &
113 : ! isotot, isotot0 ! for diagnostics
114 :
115 : real (kind=dbl_kind) :: &
116 : hs_old, hi_old, sloss1, &
117 : TsfK, & ! snow/ice temperature (K)
118 : alphai, & ! ice/vapour fractionation coefficient
119 : ratio, & ! isotopic ratio
120 : work, & ! temporary variable
121 : alpha
122 :
123 : integer (kind=int_kind) :: k
124 :
125 : character(len=*),parameter :: subname='(update_isotope)'
126 :
127 : ! initialize
128 :
129 225452 : hs_old=vsno_old/aice_old
130 225452 : hi_old=vice_old/aice_old
131 :
132 225452 : dzsno = hs_old
133 225452 : dzice = hi_old
134 :
135 225452 : if (aicen > puny) then
136 225452 : hs = vsnon/aicen
137 225452 : hi = vicen/aicen
138 0 : elseif (aice_old > puny) then
139 0 : hs = vsnon/aice_old
140 0 : hi = vicen/aice_old
141 : endif
142 :
143 225452 : if (ktherm == 2) then
144 225452 : dhs_snoice = snoice
145 : else
146 0 : dhs_snoice = snoice*rhoi/rhos
147 : endif
148 :
149 : ! if (hs > puny) then
150 : ! evaps = evap*dt/rhos
151 : ! else
152 : ! evapi = evap*dt/rhoi
153 : ! endif
154 225452 : evaps = hs - (hs_old - melts - dhs_snoice + fsnow/rhos*dt)
155 225452 : evapi = hi - (hi_old - meltt - meltb + congel + snoice)
156 :
157 : ! condensation of vapor onto snow and ice
158 :
159 225452 : TsfK = Tsfc + Tffresh
160 :
161 225452 : if (evaps > c0) then ! condensation to snow
162 3280 : do k = 1, n_iso
163 2460 : ratio = c1 ! ratio between 18O(HDO) and 16O in humidity
164 2460 : alphai = c1 ! fractionation coefficient
165 2460 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
166 0 : ratio = Qref_iso(k)/Qref_iso(2)
167 2460 : if (Qref_iso(2) <= puny) ratio = c0
168 2460 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
169 2460 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
170 2460 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
171 2460 : work = alphai*ratio*rhos*evaps*aicen
172 2460 : fiso_evapn(k) = fiso_evapn(k) + work/dt
173 3280 : isosno(k) = isosno(k) + work
174 : enddo
175 820 : dzsno = dzsno + evaps
176 : endif
177 :
178 225452 : if (evapi > c0) then ! condensation to ice
179 186496 : do k = 1, n_iso
180 139872 : ratio = c1 ! ratio between 18O(HDO) and 16O in ref humidity
181 139872 : alphai = c1 ! fractionation coefficient
182 139872 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
183 0 : ratio = Qref_iso(k)/Qref_iso(2)
184 139872 : if (Qref_iso(2) <= puny) ratio = c0
185 139872 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
186 139872 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
187 139872 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
188 139872 : work = alphai*ratio*rhoi*evapi*aicen
189 139872 : fiso_evapn(k) = fiso_evapn(k) + work/dt
190 186496 : isoice(k) = isoice(k) + work
191 : enddo
192 46624 : dzice = dzice + evapi
193 : endif
194 :
195 : ! basal ice growth and isotope uptake
196 :
197 225452 : if (congel > c0) then
198 487084 : do k = 1,n_iso
199 365313 : work = c0
200 365313 : if (k == 1) then
201 121771 : alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method)
202 121771 : work = alpha*HDO_ocn*rhoi*congel*aicen
203 243542 : elseif (k == 2) then
204 121771 : alpha = isoice_alpha(congel/dt,'H2_16O',isotope_frac_method)
205 121771 : work = alpha*H2_16O_ocn*rhoi*congel*aicen
206 121771 : elseif (k == 3) then
207 121771 : alpha = isoice_alpha(congel/dt,'H2_18O',isotope_frac_method)
208 121771 : 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 365313 : isoice(k) = isoice(k) + work
215 487084 : fiso_ocnn(k) = fiso_ocnn(k) - work/dt
216 : enddo
217 :
218 121771 : dzice = dzice + congel
219 : endif
220 :
221 : ! sublimation of snow and ice
222 :
223 225452 : if (evaps < c0) then ! snow sublimation (no fractionation)
224 891388 : 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 668541 : sloss1 = c0
232 668541 : if (dzsno > puny) sloss1 = isosno(k)*min(-evaps,dzsno)/dzsno
233 668541 : if (isosno(k) >= sloss1) then
234 663726 : isosno(k) = isosno(k)-sloss1
235 : else
236 4815 : sloss1 = isosno(k)
237 4815 : isosno(k) = c0
238 : endif
239 : ! if (isosno(k) < c0) then
240 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
241 : ! endif
242 891388 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
243 : enddo
244 :
245 222847 : dzsno = dzsno + evaps
246 222847 : if (dzsno <= c0) then ! snow goes away
247 95952 : fiso_evapn(:) = fiso_evapn(:) - isosno(:)/dt
248 95952 : isosno(:) = c0
249 23988 : dzsno = c0
250 : endif
251 : endif
252 :
253 225452 : if (evapi < c0) then ! ice sublimation (no fractionation)
254 323056 : 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 242292 : sloss1 = c0
262 242292 : if (dzice > puny) &
263 242292 : sloss1 = isoice(k) * min(-evapi,dzice)/dzice
264 242292 : if (isoice(k) >= sloss1) then
265 242292 : isoice(k) = isoice(k)-sloss1
266 : else
267 0 : sloss1 = isoice(k)
268 0 : isoice(k) = c0
269 : endif
270 323056 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
271 : enddo
272 :
273 80764 : dzice = dzice + evapi
274 80764 : 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 225452 : if (melts > c0) then
284 310200 : do k=1,n_iso
285 232650 : sloss1=c0
286 232650 : if (dzsno > puny) &
287 232650 : sloss1 = isosno(k) * min(melts,dzsno)/dzsno
288 232650 : if (isosno(k) >= sloss1) then
289 232621 : isosno(k) = isosno(k)-sloss1
290 : else
291 29 : sloss1 = isosno(k)
292 29 : isosno(k) = c0
293 : endif
294 : ! if (isosno(k) < c0) then
295 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
296 : ! endif
297 310200 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
298 : enddo ! n_iso
299 :
300 77550 : dzsno = dzsno - melts
301 77550 : if (dzsno <= c0) then ! snow melts away
302 816 : fiso_ocnn(:) = fiso_ocnn(:) + isosno(:)/dt
303 816 : isosno(:) = c0
304 204 : dzsno = c0
305 : endif
306 : endif
307 :
308 : ! surface ice melt
309 225452 : if (meltt > c0) then
310 74892 : do k=1,n_iso
311 56169 : sloss1=c0
312 56169 : if (dzice > puny) sloss1=isoice(k) * min(meltt,dzice)/dzice
313 56169 : if (isoice(k) >= sloss1) then
314 56169 : isoice(k) = isoice(k)-sloss1
315 : else
316 0 : sloss1 = isoice(k)
317 0 : isoice(k) = c0
318 : endif
319 74892 : fiso_ocnn(k)=fiso_ocnn(k) + sloss1/dt
320 : enddo
321 :
322 18723 : dzice = dzice - meltt
323 18723 : 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 225452 : if (meltb > c0) then
333 414724 : do k=1,n_iso
334 311043 : sloss1=c0
335 311043 : if (dzice > puny) sloss1=max(meltb-dzice,c0) * isoice(k)/dzice
336 311043 : if (isoice(k) >= sloss1) then
337 311043 : isoice(k) = isoice(k)-sloss1
338 : else
339 0 : sloss1 = isoice(k)
340 0 : isoice(k) = c0
341 : endif
342 414724 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
343 : enddo
344 :
345 103681 : dzice = dzice - meltb
346 103681 : if (dzice <= c0) then ! ice ice melts away
347 0 : fiso_ocnn(:) = fiso_ocnn(:) + isoice(:)
348 0 : isoice(:) = c0
349 0 : dzice = c0
350 : endif
351 : endif
352 :
353 : ! snowfall and isotope deposition
354 :
355 225452 : if (fsnow > c0) then
356 863740 : isosno(:) = isosno(:) + fiso_atm(:)*aicen*dt
357 215935 : dzsno = dzsno + fsnow/rhos*dt
358 : endif
359 :
360 : ! snoice formation
361 :
362 225452 : if (dhs_snoice > c0) then
363 15264 : do k=1,n_iso
364 11448 : sloss1=c0
365 11448 : if (dzsno > puny) sloss1 = min(dhs_snoice,dzsno) * isosno(k)/dzsno
366 11448 : if (isosno(k) >= sloss1) then
367 11448 : 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 15264 : isoice(k) = isoice(k) + sloss1
376 : enddo
377 :
378 3816 : dzsno = dzsno - dhs_snoice
379 3816 : dzice = dzice + snoice
380 3816 : 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 225452 : if (aicen > puny) then
404 901808 : fiso_ocnn(:) = fiso_ocnn(:)/aicen
405 901808 : 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 434865 : 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 434865 : if (sp == 'H2_16O') isoice_alpha = c1
433 :
434 : ! Lehmann and Siegenthaler, 1991
435 : !--------------------------------------------------
436 434865 : if (frac == 'cfrac' .and. sp == 'HDO') &
437 144955 : isoice_alpha = 1.02120_dbl_kind
438 434865 : if (frac == 'cfrac' .and. sp == 'H2_18O') &
439 144955 : 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 434865 : 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 434865 : 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 434865 : return
453 :
454 434865 : end function isoice_alpha
455 :
456 : !=======================================================================
457 :
458 142332 : 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 142332 : if (isp == isph2o) then
479 0 : wiso_alpi = c1
480 0 : return
481 : end if
482 :
483 142332 : wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp))
484 :
485 142332 : return
486 : end function wiso_alpi
487 :
488 : !=======================================================================
489 :
490 : end module icepack_isotope
491 :
492 : !=======================================================================
|
