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