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 15774367 : subroutine update_isotope(dt, &
49 : meltt, melts, & ! LCOV_EXCL_LINE
50 : meltb, congel, & ! LCOV_EXCL_LINE
51 : snoice, evap, & ! LCOV_EXCL_LINE
52 : fsnow, Tsfc, & ! LCOV_EXCL_LINE
53 15774367 : Qref_iso, & ! LCOV_EXCL_LINE
54 15774367 : isosno, isoice, & ! LCOV_EXCL_LINE
55 : aice_old, & ! LCOV_EXCL_LINE
56 : vice_old, vsno_old, & ! LCOV_EXCL_LINE
57 : vicen, vsnon, aicen, & ! LCOV_EXCL_LINE
58 31548734 : fiso_atm, fiso_evapn, & ! LCOV_EXCL_LINE
59 15774367 : fiso_ocnn, HDO_ocn, & ! LCOV_EXCL_LINE
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 ! LCOV_EXCL_LINE
70 : meltt, & ! top melt ! LCOV_EXCL_LINE
71 : melts, & ! snow melt ! LCOV_EXCL_LINE
72 : meltb, & ! bottom melt ! LCOV_EXCL_LINE
73 : congel, & ! congelation ice growth (m/step) ! LCOV_EXCL_LINE
74 : snoice, & ! ice thickness increase (m/step) ! LCOV_EXCL_LINE
75 : evap, & ! surface evaporation ! LCOV_EXCL_LINE
76 : fsnow, & ! snowfall (kg/m^2/s of water) ! LCOV_EXCL_LINE
77 : vicen, & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
78 : vsnon, & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
79 : aicen, & ! ice area ! LCOV_EXCL_LINE
80 : aice_old, & ! beginning values ! LCOV_EXCL_LINE
81 : vice_old, & ! LCOV_EXCL_LINE
82 : vsno_old
83 :
84 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
85 : fiso_ocnn, & ! isotopic freshwater (kg/m^2/s) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
90 : isoice
91 :
92 : real (kind=dbl_kind), dimension(:), intent(in) :: &
93 : fiso_atm, & ! isotopic snowfall (kg/m^2/s of water) ! LCOV_EXCL_LINE
94 : Qref_iso ! isotope reference humidity
95 :
96 : real (kind=dbl_kind), intent(in) :: &
97 : HDO_ocn, & ! ocean concentration of HDO (kg/kg) ! LCOV_EXCL_LINE
98 : H2_16O_ocn, & ! ocean concentration of H2_16O (kg/kg) ! LCOV_EXCL_LINE
99 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
100 :
101 : ! local variables
102 :
103 : real (kind=dbl_kind) :: &
104 : dzsno, & ! LCOV_EXCL_LINE
105 : dzice, & ! LCOV_EXCL_LINE
106 : evaps, & ! evaporation over snow (m/step) ! LCOV_EXCL_LINE
107 : evapi, & ! evaporation over ice (m/step) ! LCOV_EXCL_LINE
108 : dhs_snoice,& ! snow thickness reduction (m/step) ! LCOV_EXCL_LINE
109 : hi, & ! LCOV_EXCL_LINE
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, & ! LCOV_EXCL_LINE
117 : TsfK, & ! snow/ice temperature (K) ! LCOV_EXCL_LINE
118 : alphai, & ! ice/vapour fractionation coefficient ! LCOV_EXCL_LINE
119 : ratio, & ! isotopic ratio ! LCOV_EXCL_LINE
120 : work, & ! temporary variable ! LCOV_EXCL_LINE
121 : alpha
122 :
123 : integer (kind=int_kind) :: k
124 :
125 : character(len=*),parameter :: subname='(update_isotope)'
126 :
127 : ! initialize
128 :
129 15774367 : hs_old=vsno_old/aice_old
130 15774367 : hi_old=vice_old/aice_old
131 :
132 15774367 : dzsno = hs_old
133 15774367 : dzice = hi_old
134 :
135 15774367 : if (aicen > puny) then
136 15771831 : hs = vsnon/aicen
137 15771831 : hi = vicen/aicen
138 2536 : elseif (aice_old > puny) then
139 2536 : hs = vsnon/aice_old
140 2536 : hi = vicen/aice_old
141 : endif
142 :
143 15774367 : if (ktherm == 2) then
144 15774367 : 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 15774367 : evaps = hs - (hs_old - melts - dhs_snoice + fsnow/rhos*dt)
155 15774367 : evapi = hi - (hi_old - meltt - meltb + congel + snoice)
156 :
157 : ! condensation of vapor onto snow and ice
158 :
159 15774367 : TsfK = Tsfc + Tffresh
160 :
161 15774367 : if (evaps > c0) then ! condensation to snow
162 24540240 : do k = 1, n_iso
163 18405180 : ratio = c1 ! ratio between 18O(HDO) and 16O in humidity
164 18405180 : alphai = c1 ! fractionation coefficient
165 18405180 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
166 0 : ratio = Qref_iso(k)/Qref_iso(2)
167 18405180 : if (Qref_iso(2) <= puny) ratio = c0
168 18405180 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
169 18405180 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
170 18405180 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
171 18405180 : work = alphai*ratio*rhos*evaps*aicen
172 18405180 : fiso_evapn(k) = fiso_evapn(k) + work/dt
173 24540240 : isosno(k) = isosno(k) + work
174 : enddo
175 6135060 : dzsno = dzsno + evaps
176 : endif
177 :
178 15774367 : if (evapi > c0) then ! condensation to ice
179 17684080 : do k = 1, n_iso
180 13263060 : ratio = c1 ! ratio between 18O(HDO) and 16O in ref humidity
181 13263060 : alphai = c1 ! fractionation coefficient
182 13263060 : if (isotope_frac_method.ne.'nfrac' .and. Qref_iso(2)>puny) &
183 0 : ratio = Qref_iso(k)/Qref_iso(2)
184 13263060 : if (Qref_iso(2) <= puny) ratio = c0
185 13263060 : if (isotope_frac_method.ne.'nfrac' .and. k==1) alphai = wiso_alpi(3,TsfK)
186 13263060 : if (isotope_frac_method.ne.'nfrac' .and. k==2) alphai = wiso_alpi(2,TsfK)
187 13263060 : if (isotope_frac_method.ne.'nfrac' .and. k==3) alphai = wiso_alpi(4,TsfK)
188 13263060 : work = alphai*ratio*rhoi*evapi*aicen
189 13263060 : fiso_evapn(k) = fiso_evapn(k) + work/dt
190 17684080 : isoice(k) = isoice(k) + work
191 : enddo
192 4421020 : dzice = dzice + evapi
193 : endif
194 :
195 : ! basal ice growth and isotope uptake
196 :
197 15774367 : if (congel > c0) then
198 17997604 : do k = 1,n_iso
199 13498203 : work = c0
200 13498203 : if (k == 1) then
201 4499401 : alpha = isoice_alpha(congel/dt,'HDO',isotope_frac_method)
202 4499401 : work = alpha*HDO_ocn*rhoi*congel*aicen
203 8998802 : elseif (k == 2) then
204 4499401 : alpha = isoice_alpha(congel/dt,'H2_16O',isotope_frac_method)
205 4499401 : work = alpha*H2_16O_ocn*rhoi*congel*aicen
206 4499401 : elseif (k == 3) then
207 4499401 : alpha = isoice_alpha(congel/dt,'H2_18O',isotope_frac_method)
208 4499401 : 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 13498203 : isoice(k) = isoice(k) + work
215 17997604 : fiso_ocnn(k) = fiso_ocnn(k) - work/dt
216 : enddo
217 :
218 4499401 : dzice = dzice + congel
219 : endif
220 :
221 : ! sublimation of snow and ice
222 :
223 15774367 : if (evaps < c0) then ! snow sublimation (no fractionation)
224 35331424 : 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 26498568 : sloss1 = c0
232 26498568 : if (dzsno > puny) sloss1 = isosno(k)*min(-evaps,dzsno)/dzsno
233 26498568 : if (isosno(k) >= sloss1) then
234 26489166 : isosno(k) = isosno(k)-sloss1
235 : else
236 9402 : sloss1 = isosno(k)
237 9402 : isosno(k) = c0
238 : endif
239 : ! if (isosno(k) < c0) then
240 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
241 : ! endif
242 35331424 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
243 : enddo
244 :
245 8832856 : dzsno = dzsno + evaps
246 8832856 : if (dzsno <= c0) then ! snow goes away
247 307924 : fiso_evapn(:) = fiso_evapn(:) - isosno(:)/dt
248 307924 : isosno(:) = c0
249 76981 : dzsno = c0
250 : endif
251 : endif
252 :
253 15774367 : if (evapi < c0) then ! ice sublimation (no fractionation)
254 15625236 : 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 11718927 : sloss1 = c0
262 11718927 : if (dzice > puny) &
263 11718927 : sloss1 = isoice(k) * min(-evapi,dzice)/dzice
264 11718927 : if (isoice(k) >= sloss1) then
265 11718927 : isoice(k) = isoice(k)-sloss1
266 : else
267 0 : sloss1 = isoice(k)
268 0 : isoice(k) = c0
269 : endif
270 15625236 : fiso_evapn(k) = fiso_evapn(k) - sloss1/dt
271 : enddo
272 :
273 3906309 : dzice = dzice + evapi
274 3906309 : 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 15774367 : if (melts > c0) then
284 22724240 : do k=1,n_iso
285 17043180 : sloss1=c0
286 17043180 : if (dzsno > puny) &
287 16995957 : sloss1 = isosno(k) * min(melts,dzsno)/dzsno
288 17043180 : if (isosno(k) >= sloss1) then
289 16984053 : isosno(k) = isosno(k)-sloss1
290 : else
291 59127 : sloss1 = isosno(k)
292 59127 : isosno(k) = c0
293 : endif
294 : ! if (isosno(k) < c0) then
295 : ! write(nu_diag,*) 'Neg isosno(k)',isosno(k),sloss1
296 : ! endif
297 22724240 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
298 : enddo ! n_iso
299 :
300 5681060 : dzsno = dzsno - melts
301 5681060 : if (dzsno <= c0) then ! snow melts away
302 1733064 : fiso_ocnn(:) = fiso_ocnn(:) + isosno(:)/dt
303 1733064 : isosno(:) = c0
304 433266 : dzsno = c0
305 : endif
306 : endif
307 :
308 : ! surface ice melt
309 15774367 : if (meltt > c0) then
310 5856908 : do k=1,n_iso
311 4392681 : sloss1=c0
312 4392681 : if (dzice > puny) sloss1=isoice(k) * min(meltt,dzice)/dzice
313 4392681 : if (isoice(k) >= sloss1) then
314 4392645 : isoice(k) = isoice(k)-sloss1
315 : else
316 36 : sloss1 = isoice(k)
317 36 : isoice(k) = c0
318 : endif
319 5856908 : fiso_ocnn(k)=fiso_ocnn(k) + sloss1/dt
320 : enddo
321 :
322 1464227 : dzice = dzice - meltt
323 1464227 : if (dzice <= c0) then ! ice ice melts away
324 1188 : fiso_ocnn(:) = fiso_ocnn(:)+isoice(:)
325 1188 : isoice(:) = c0
326 297 : dzice = c0
327 : endif
328 : endif
329 :
330 : ! basal ice melt. Assume all isotopes lost in basal melt
331 :
332 15774367 : if (meltb > c0) then
333 45098680 : do k=1,n_iso
334 33824010 : sloss1=c0
335 33824010 : if (dzice > puny) sloss1=max(meltb-dzice,c0) * isoice(k)/dzice
336 33824010 : if (isoice(k) >= sloss1) then
337 33823902 : isoice(k) = isoice(k)-sloss1
338 : else
339 108 : sloss1 = isoice(k)
340 108 : isoice(k) = c0
341 : endif
342 45098680 : fiso_ocnn(k) = fiso_ocnn(k) + sloss1/dt
343 : enddo
344 :
345 11274670 : dzice = dzice - meltb
346 11274670 : if (dzice <= c0) then ! ice ice melts away
347 14772 : fiso_ocnn(:) = fiso_ocnn(:) + isoice(:)
348 14772 : isoice(:) = c0
349 3693 : dzice = c0
350 : endif
351 : endif
352 :
353 : ! snowfall and isotope deposition
354 :
355 15774367 : if (fsnow > c0) then
356 36049516 : isosno(:) = isosno(:) + fiso_atm(:)*aicen*dt
357 9012379 : dzsno = dzsno + fsnow/rhos*dt
358 : endif
359 :
360 : ! snoice formation
361 :
362 15774367 : if (dhs_snoice > c0) then
363 3288144 : do k=1,n_iso
364 2466108 : sloss1=c0
365 2466108 : if (dzsno > puny) sloss1 = min(dhs_snoice,dzsno) * isosno(k)/dzsno
366 2466108 : if (isosno(k) >= sloss1) then
367 2465946 : isosno(k) = isosno(k)-sloss1
368 : else
369 162 : sloss1 = isosno(k)
370 162 : 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 3288144 : isoice(k) = isoice(k) + sloss1
376 : enddo
377 :
378 822036 : dzsno = dzsno - dhs_snoice
379 822036 : dzice = dzice + snoice
380 822036 : if (dzsno <= c0) then ! snow goes away
381 4260 : fiso_ocnn(:)= fiso_ocnn(:) + isosno(:)/dt
382 4260 : isosno(:) = c0
383 1065 : 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 & ! LCOV_EXCL_LINE
392 : ! - fiso_evapn(k)*dt & ! LCOV_EXCL_LINE
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 & ! LCOV_EXCL_LINE
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 15774367 : if (aicen > puny) then
404 63087324 : fiso_ocnn(:) = fiso_ocnn(:)/aicen
405 63087324 : 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 32709558 : 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 32709558 : if (sp == 'H2_16O') isoice_alpha = c1
433 :
434 : ! Lehmann and Siegenthaler, 1991
435 : !--------------------------------------------------
436 32709558 : if (frac == 'cfrac' .and. sp == 'HDO') &
437 10903186 : isoice_alpha = 1.02120_dbl_kind
438 32709558 : if (frac == 'cfrac' .and. sp == 'H2_18O') &
439 10903186 : 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 32709558 : if (frac == 'gfrac' .and. sp == 'HDO') &
445 : isoice_alpha = c1+7.2852_dbl_kind*1.2280E-3_dbl_kind+ & ! LCOV_EXCL_LINE
446 : 0.7311E-3_dbl_kind*exp(-growth_rate/8.0100E8_dbl_kind)+ & ! LCOV_EXCL_LINE
447 0 : 0.8441E-3_dbl_kind*exp(-growth_rate/0.7800E6_dbl_kind)
448 32709558 : if (frac == 'gfrac' .and. sp == 'H2_18O') &
449 : isoice_alpha = c1+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 32709558 : return
453 :
454 32709558 : end function isoice_alpha
455 :
456 : !=======================================================================
457 :
458 31668240 : 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 ! LCOV_EXCL_LINE
473 : alpai = (/ 0._dbl_kind, 0._dbl_kind, 16289._dbl_kind, 0._dbl_kind /), & ! LCOV_EXCL_LINE
474 : alpbi = (/ 0._dbl_kind, 0._dbl_kind, 0._dbl_kind, 11.839_dbl_kind /), & ! LCOV_EXCL_LINE
475 : alpci = (/ 0._dbl_kind, 0._dbl_kind, -9.45e-2_dbl_kind, -28.224e-3_dbl_kind /)
476 :
477 : !-----------------------------------------------------------------------
478 31668240 : if (isp == isph2o) then
479 0 : wiso_alpi = c1
480 0 : return
481 : end if
482 :
483 31668240 : wiso_alpi = exp(alpai(isp)/tk**2 + alpbi(isp)/tk + alpci(isp))
484 :
485 31668240 : return
486 : end function wiso_alpi
487 :
488 : !=======================================================================
489 :
490 : end module icepack_isotope
491 :
492 : !=======================================================================
|
