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