Line data Source code
1 : !=======================================================================
2 :
3 : ! Level-ice meltpond parameterization
4 : !
5 : ! This meltpond parameterization was developed for use with the delta-
6 : ! Eddington radiation scheme, and only affects the radiation budget in
7 : ! the model. That is, although the pond volume is tracked, that liquid
8 : ! water is not used elsewhere in the model for mass budgets or other
9 : ! physical processes.
10 : !
11 : ! authors Elizabeth Hunke (LANL)
12 : ! David Hebert (NRL Stennis)
13 : ! Olivier Lecomte (Univ. Louvain)
14 :
15 : module icepack_meltpond_lvl
16 :
17 : use icepack_kinds
18 : use icepack_parameters, only: c0, c1, c2, c10, p01, p5, puny
19 : use icepack_parameters, only: viscosity_dyn, rhoi, rhos, rhow, Timelt, Tffresh, Lfresh
20 : use icepack_parameters, only: gravit, depressT, rhofresh, kice
21 : use icepack_warnings, only: warnstr, icepack_warnings_add
22 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
23 :
24 : implicit none
25 :
26 : private
27 : public :: compute_ponds_lvl
28 :
29 : !=======================================================================
30 :
31 : contains
32 :
33 : !=======================================================================
34 :
35 4305360 : subroutine compute_ponds_lvl(dt, nilyr, &
36 : ktherm, &
37 : hi_min, dpscale, &
38 0 : frzpnd, pndaspect, &
39 : rfrac, meltt, melts, &
40 : frain, Tair, fsurfn,&
41 : dhs, ffrac, &
42 : aicen, vicen, vsnon, &
43 4305360 : qicen, sicen, &
44 : Tsfcn, alvl, &
45 : apnd, hpnd, ipnd)
46 :
47 : integer (kind=int_kind), intent(in) :: &
48 : nilyr, & ! number of ice layers
49 : ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
50 :
51 : real (kind=dbl_kind), intent(in) :: &
52 : dt, & ! time step (s)
53 : hi_min, & ! minimum ice thickness allowed for thermo (m)
54 : dpscale, & ! alter e-folding time scale for flushing
55 : pndaspect ! ratio of pond depth to pond fraction
56 :
57 : character (len=char_len), intent(in) :: &
58 : frzpnd ! pond refreezing parameterization
59 :
60 : real (kind=dbl_kind), &
61 : intent(in) :: &
62 : Tsfcn, & ! surface temperature (C)
63 : alvl, & ! fraction of level ice
64 : rfrac, & ! water fraction retained for melt ponds
65 : meltt, & ! top melt rate (m/s)
66 : melts, & ! snow melt rate (m/s)
67 : frain, & ! rainfall rate (kg/m2/s)
68 : Tair, & ! air temperature (K)
69 : fsurfn,& ! atm-ice surface heat flux (W/m2)
70 : aicen, & ! ice area fraction
71 : vicen, & ! ice volume (m)
72 : vsnon ! snow volume (m)
73 :
74 : real (kind=dbl_kind), &
75 : intent(inout) :: &
76 : apnd, hpnd, ipnd
77 :
78 : real (kind=dbl_kind), dimension (:), intent(in) :: &
79 : qicen, & ! ice layer enthalpy (J m-3)
80 : sicen ! salinity (ppt)
81 :
82 : real (kind=dbl_kind), &
83 : intent(in) :: &
84 : dhs ! depth difference for snow on sea ice and pond ice
85 :
86 : real (kind=dbl_kind), &
87 : intent(out) :: &
88 : ffrac ! fraction of fsurfn over pond used to melt ipond
89 :
90 : ! local temporary variables
91 :
92 : real (kind=dbl_kind) :: &
93 1670640 : volpn ! pond volume per unit area (m)
94 :
95 : real (kind=dbl_kind), dimension (nilyr) :: &
96 19341120 : Tmlt ! melting temperature (C)
97 :
98 : real (kind=dbl_kind) :: &
99 1670640 : hi , & ! ice thickness (m)
100 1670640 : hs , & ! snow depth (m)
101 1670640 : dTs , & ! surface temperature diff for freeze-up (C)
102 1670640 : Tp , & ! pond freezing temperature (C)
103 1670640 : Ts , & ! surface air temperature (C)
104 1670640 : apondn , & ! local pond area
105 1670640 : hpondn , & ! local pond depth (m)
106 1670640 : dvn , & ! change in pond volume (m)
107 3341280 : hlid, alid , & ! refrozen lid thickness, area
108 1670640 : dhlid , & ! change in refrozen lid thickness
109 1670640 : bdt , & ! 2 kice dT dt / (rhoi Lfresh)
110 1670640 : alvl_tmp , & ! level ice fraction of ice area
111 5011920 : draft, deltah, pressure_head, perm, drain ! for permeability
112 :
113 : real (kind=dbl_kind), parameter :: &
114 : Td = c2 , & ! temperature difference for freeze-up (C)
115 : rexp = p01 ! pond contraction scaling
116 :
117 : character(len=*),parameter :: subname='(compute_ponds_lvl)'
118 :
119 : !-----------------------------------------------------------------
120 : ! Initialize
121 : !-----------------------------------------------------------------
122 :
123 4305360 : volpn = hpnd * aicen * alvl * apnd
124 4305360 : ffrac = c0
125 :
126 : !-----------------------------------------------------------------
127 : ! Identify grid cells where ponds can be
128 : !-----------------------------------------------------------------
129 :
130 4305360 : if (aicen*alvl > puny**2) then
131 :
132 2758086 : hi = vicen/aicen
133 2758086 : hs = vsnon/aicen
134 2758086 : alvl_tmp = alvl
135 :
136 2758086 : if (hi < hi_min) then
137 :
138 : !--------------------------------------------------------------
139 : ! Remove ponds on thin ice
140 : !--------------------------------------------------------------
141 8145 : apondn = c0
142 8145 : hpondn = c0
143 8145 : volpn = c0
144 8145 : hlid = c0
145 :
146 : else
147 :
148 : !-----------------------------------------------------------
149 : ! initialize pond area as fraction of ice
150 : !-----------------------------------------------------------
151 2749941 : apondn = apnd*alvl_tmp
152 :
153 : !-----------------------------------------------------------
154 : ! update pond volume
155 : !-----------------------------------------------------------
156 : ! add melt water
157 : dvn = rfrac/rhofresh*(meltt*rhoi &
158 : + melts*rhos &
159 2749941 : + frain* dt)*aicen
160 :
161 : ! shrink pond volume under freezing conditions
162 2749941 : if (trim(frzpnd) == 'cesm') then
163 261523 : Tp = Timelt - Td
164 261523 : dTs = max(Tp - Tsfcn,c0)
165 261523 : dvn = dvn - volpn * (c1 - exp(rexp*dTs/Tp))
166 :
167 : else
168 : ! trim(frzpnd) == 'hlid' Stefan approximation
169 : ! assumes pond is fresh (freezing temperature = 0 C)
170 : ! and ice grows from existing pond ice
171 2488418 : hlid = ipnd
172 2488418 : if (dvn == c0) then ! freeze pond
173 1463019 : Ts = Tair - Tffresh
174 1463019 : if (Ts < c0) then
175 : ! if (Ts < -c2) then ! as in meltpond_cesm
176 1463019 : bdt = -c2*Ts*kice*dt/(rhoi*Lfresh)
177 1463019 : dhlid = p5*sqrt(bdt) ! open water freezing
178 1463019 : if (hlid > dhlid) dhlid = p5*bdt/hlid ! existing ice
179 1463019 : dhlid = min(dhlid, hpnd*rhofresh/rhoi)
180 1463019 : hlid = hlid + dhlid
181 : else
182 0 : dhlid = c0 ! to account for surface inversions
183 : endif
184 : else ! convert refrozen pond ice back to water
185 1025399 : dhlid = max(fsurfn*dt / (rhoi*Lfresh), c0) ! > 0
186 1025399 : dhlid = -min(dhlid, hlid) ! < 0
187 1025399 : hlid = max(hlid + dhlid, c0)
188 1025399 : if (hs - dhs < puny) then ! pond ice is snow-free
189 52858 : ffrac = c1 ! fraction of fsurfn over pond used to melt ipond
190 52858 : if (fsurfn > puny) &
191 44163 : ffrac = min(-dhlid*rhoi*Lfresh/(dt*fsurfn), c1)
192 : endif
193 : endif
194 2488418 : alid = apondn * aicen
195 2488418 : dvn = dvn - dhlid*alid*rhoi/rhofresh
196 : endif
197 :
198 2749941 : volpn = volpn + dvn
199 :
200 : !-----------------------------------------------------------
201 : ! update pond area and depth
202 : !-----------------------------------------------------------
203 2749941 : if (volpn <= c0) then
204 1457026 : volpn = c0
205 1457026 : apondn = c0
206 : endif
207 :
208 2749941 : if (apondn*aicen > puny) then ! existing ponds
209 : apondn = max(c0, min(alvl_tmp, &
210 774853 : apondn + 0.5*dvn/(pndaspect*apondn*aicen)))
211 774853 : hpondn = c0
212 774853 : if (apondn > puny) &
213 772955 : hpondn = volpn/(apondn*aicen)
214 :
215 1975088 : elseif (alvl_tmp*aicen > c10*puny) then ! new ponds
216 1961388 : apondn = min (sqrt(volpn/(pndaspect*aicen)), alvl_tmp)
217 1961388 : hpondn = pndaspect * apondn
218 :
219 : else ! melt water runs off deformed ice
220 13700 : apondn = c0
221 13700 : hpondn = c0
222 : endif
223 2749941 : apondn = max(apondn, c0)
224 :
225 : ! limit pond depth to maintain nonnegative freeboard
226 2749941 : hpondn = min(hpondn, ((rhow-rhoi)*hi - rhos*hs)/rhofresh)
227 :
228 : ! fraction of grid cell covered by ponds
229 2749941 : apondn = apondn * aicen
230 :
231 2749941 : volpn = hpondn*apondn
232 2749941 : if (volpn <= c0) then
233 1468440 : volpn = c0
234 1468440 : apondn = c0
235 1468440 : hpondn = c0
236 1468440 : hlid = c0
237 : endif
238 :
239 : !-----------------------------------------------------------
240 : ! drainage due to permeability (flushing)
241 : ! setting dpscale = 0 turns this off
242 : ! NOTE this uses the initial salinity and melting T profiles
243 : !-----------------------------------------------------------
244 :
245 2749941 : if (ktherm /= 2 .and. hpondn > c0 .and. dpscale > puny) then
246 361363 : draft = (rhos*hs + rhoi*hi)/rhow + hpondn
247 361363 : deltah = hpondn + hi - draft
248 361363 : pressure_head = gravit * rhow * max(deltah, c0)
249 2890904 : Tmlt(:) = -sicen(:) * depressT
250 0 : call brine_permeability(nilyr, qicen, &
251 361363 : sicen, Tmlt, perm)
252 361363 : if (icepack_warnings_aborted(subname)) return
253 361363 : drain = perm*pressure_head*dt / (viscosity_dyn*hi) * dpscale
254 361363 : deltah = min(drain, hpondn)
255 361363 : dvn = -deltah*apondn
256 361363 : volpn = volpn + dvn
257 : apondn = max(c0, min(apondn &
258 361363 : + 0.5*dvn/(pndaspect*apondn), alvl_tmp*aicen))
259 361363 : hpondn = c0
260 361363 : if (apondn > puny) hpondn = volpn/apondn
261 : endif
262 :
263 : endif
264 :
265 : !-----------------------------------------------------------
266 : ! Reload tracer array
267 : !-----------------------------------------------------------
268 :
269 2758086 : hpnd = hpondn
270 2758086 : apnd = apondn / (aicen*alvl_tmp)
271 2758086 : if (trim(frzpnd) == 'hlid') ipnd = hlid
272 :
273 : endif
274 :
275 : end subroutine compute_ponds_lvl
276 :
277 : !=======================================================================
278 :
279 : ! determine the liquid fraction of brine in the ice and the permeability
280 :
281 361363 : subroutine brine_permeability(nilyr, qicen, salin, Tmlt, perm)
282 :
283 : use icepack_therm_shared, only: calculate_Tin_from_qin
284 :
285 : integer (kind=int_kind), intent(in) :: &
286 : nilyr ! number of ice layers
287 :
288 : real (kind=dbl_kind), dimension(:), intent(in) :: &
289 : qicen, & ! enthalpy for each ice layer (J m-3)
290 : salin, & ! salinity (ppt)
291 : Tmlt ! melting temperature (C)
292 :
293 : real (kind=dbl_kind), intent(out) :: &
294 : perm ! permeability (m^2)
295 :
296 : ! local variables
297 :
298 : real (kind=dbl_kind) :: &
299 129283 : Sbr ! brine salinity
300 :
301 : real (kind=dbl_kind), dimension(nilyr) :: &
302 1498424 : Tin, & ! ice temperature (C)
303 1627707 : phi ! liquid fraction
304 :
305 : integer (kind=int_kind) :: k
306 :
307 : character(len=*),parameter :: subname='(brine_permeability)'
308 :
309 : !-----------------------------------------------------------------
310 : ! Compute ice temperatures from enthalpies using quadratic formula
311 : !-----------------------------------------------------------------
312 :
313 2890904 : do k = 1,nilyr
314 2890904 : Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt(k))
315 : enddo
316 :
317 : !-----------------------------------------------------------------
318 : ! brine salinity and liquid fraction
319 : !-----------------------------------------------------------------
320 :
321 2890904 : do k = 1,nilyr
322 2529541 : Sbr = c1/(1.e-3_dbl_kind - depressT/Tin(k)) ! Notz thesis eq 3.6
323 2529541 : phi(k) = salin(k)/Sbr ! liquid fraction
324 2890904 : if (phi(k) < 0.05) phi(k) = c0 ! impermeable
325 : enddo
326 :
327 : !-----------------------------------------------------------------
328 : ! permeability
329 : !-----------------------------------------------------------------
330 :
331 3252267 : perm = 3.0e-8_dbl_kind * (minval(phi))**3
332 :
333 361363 : end subroutine brine_permeability
334 :
335 : !=======================================================================
336 :
337 : end module icepack_meltpond_lvl
338 :
339 : !=======================================================================
|
