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 7046304408 : subroutine compute_ponds_lvl(dt, &
38 : rfrac, meltt, melts, & ! LCOV_EXCL_LINE
39 : frain, Tair, fsurfn,& ! LCOV_EXCL_LINE
40 : dhs, ffrac, & ! LCOV_EXCL_LINE
41 : aicen, vicen, vsnon, & ! LCOV_EXCL_LINE
42 7046304408 : qicen, sicen, & ! LCOV_EXCL_LINE
43 : Tsfcn, alvl, & ! LCOV_EXCL_LINE
44 : apnd, hpnd, ipnd, & ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
52 : alvl, & ! fraction of level ice ! LCOV_EXCL_LINE
53 : rfrac, & ! water fraction retained for melt ponds ! LCOV_EXCL_LINE
54 : meltt, & ! top melt rate (m/s) ! LCOV_EXCL_LINE
55 : melts, & ! snow melt rate (m/s) ! LCOV_EXCL_LINE
56 : frain, & ! rainfall rate (kg/m2/s) ! LCOV_EXCL_LINE
57 : Tair, & ! air temperature (K) ! LCOV_EXCL_LINE
58 : fsurfn,& ! atm-ice surface heat flux (W/m2) ! LCOV_EXCL_LINE
59 : aicen, & ! ice area fraction ! LCOV_EXCL_LINE
60 : vicen, & ! ice volume (m) ! LCOV_EXCL_LINE
61 : vsnon, & ! snow volume (m) ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
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 7046304408 : Tmlt ! melting temperature (C)
84 :
85 : real (kind=dbl_kind) :: &
86 : hi , & ! ice thickness (m) ! LCOV_EXCL_LINE
87 : hs , & ! snow depth (m) ! LCOV_EXCL_LINE
88 : dTs , & ! surface temperature diff for freeze-up (C) ! LCOV_EXCL_LINE
89 : Tp , & ! pond freezing temperature (C) ! LCOV_EXCL_LINE
90 : Ts , & ! surface air temperature (C) ! LCOV_EXCL_LINE
91 : apondn , & ! local pond area ! LCOV_EXCL_LINE
92 : hpondn , & ! local pond depth (m) ! LCOV_EXCL_LINE
93 : dvn , & ! change in pond volume (m) ! LCOV_EXCL_LINE
94 : hlid, alid , & ! refrozen lid thickness, area ! LCOV_EXCL_LINE
95 : dhlid , & ! change in refrozen lid thickness ! LCOV_EXCL_LINE
96 : bdt , & ! 2 kice dT dt / (rhoi Lfresh) ! LCOV_EXCL_LINE
97 : alvl_tmp , & ! level ice fraction of ice area ! LCOV_EXCL_LINE
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) ! LCOV_EXCL_LINE
102 : rexp = p01 ! pond contraction scaling
103 :
104 : character(len=*),parameter :: subname='(compute_ponds_lvl)'
105 :
106 : !-----------------------------------------------------------------
107 : ! Initialize
108 : !-----------------------------------------------------------------
109 :
110 7046304408 : volpn = hpnd * aicen * alvl * apnd
111 7046304408 : ffrac = c0
112 :
113 : !-----------------------------------------------------------------
114 : ! Identify grid cells where ponds can be
115 : !-----------------------------------------------------------------
116 :
117 7046304408 : if (aicen*alvl > puny**2) then
118 :
119 1352965794 : hi = vicen/aicen
120 1352965794 : hs = vsnon/aicen
121 1352965794 : alvl_tmp = alvl
122 :
123 1352965794 : if (hi < hi_min) then
124 :
125 : !--------------------------------------------------------------
126 : ! Remove ponds on thin ice
127 : !--------------------------------------------------------------
128 1815468 : apondn = c0
129 1815468 : hpondn = c0
130 1815468 : volpn = c0
131 1815468 : hlid = c0
132 :
133 : else
134 :
135 : !-----------------------------------------------------------
136 : ! initialize pond area as fraction of ice
137 : !-----------------------------------------------------------
138 1351150326 : apondn = apnd*alvl_tmp
139 :
140 : !-----------------------------------------------------------
141 : ! update pond volume
142 : !-----------------------------------------------------------
143 : ! add melt water
144 1351150326 : if (use_smliq_pnd) then
145 : dvn = rfrac/rhofresh*(meltt*rhoi &
146 0 : + meltsliqn)*aicen
147 : else
148 : dvn = rfrac/rhofresh*(meltt*rhoi &
149 : + melts*rhos & ! LCOV_EXCL_LINE
150 1351150326 : + frain* dt)*aicen
151 : endif
152 :
153 : ! shrink pond volume under freezing conditions
154 1351150326 : if (trim(frzpnd) == 'cesm') then
155 0 : Tp = Timelt - Td
156 0 : dTs = max(Tp - Tsfcn,c0)
157 0 : 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 1351150326 : hlid = ipnd
164 1351150326 : if (dvn == c0) then ! freeze pond
165 777379893 : Ts = Tair - Tffresh
166 777379893 : if (Ts < c0) then
167 : ! if (Ts < -c2) then ! as in meltpond_cesm
168 769102220 : bdt = -c2*Ts*kice*dt/(rhoi*Lfresh)
169 769102220 : dhlid = p5*sqrt(bdt) ! open water freezing
170 769102220 : if (hlid > dhlid) dhlid = p5*bdt/hlid ! existing ice
171 769102220 : dhlid = min(dhlid, hpnd*rhofresh/rhoi)
172 769102220 : hlid = hlid + dhlid
173 : else
174 8277673 : dhlid = c0 ! to account for surface inversions
175 : endif
176 : else ! convert refrozen pond ice back to water
177 573770433 : dhlid = max(fsurfn*dt / (rhoi*Lfresh), c0) ! > 0
178 573770433 : dhlid = -min(dhlid, hlid) ! < 0
179 573770433 : hlid = max(hlid + dhlid, c0)
180 573770433 : if (hs - dhs < puny) then ! pond ice is snow-free
181 122057866 : ffrac = c1 ! fraction of fsurfn over pond used to melt ipond
182 122057866 : if (fsurfn > puny) &
183 99975530 : ffrac = min(-dhlid*rhoi*Lfresh/(dt*fsurfn), c1)
184 : endif
185 : endif
186 1351150326 : alid = apondn * aicen
187 1351150326 : dvn = dvn - dhlid*alid*rhoi/rhofresh
188 : endif
189 :
190 1351150326 : volpn = volpn + dvn
191 :
192 : !-----------------------------------------------------------
193 : ! update pond area and depth
194 : !-----------------------------------------------------------
195 1351150326 : if (volpn <= c0) then
196 750440662 : volpn = c0
197 750440662 : apondn = c0
198 : endif
199 :
200 1351150326 : if (apondn*aicen > puny) then ! existing ponds
201 : apondn = max(c0, min(alvl_tmp, &
202 593322561 : apondn + 0.5*dvn/(pndaspect*apondn*aicen)))
203 593322561 : hpondn = c0
204 593322561 : if (apondn > puny) &
205 592806212 : hpondn = volpn/(apondn*aicen)
206 :
207 757827765 : elseif (alvl_tmp*aicen > c10*puny) then ! new ponds
208 741791655 : apondn = min (sqrt(volpn/(pndaspect*aicen)), alvl_tmp)
209 741791655 : hpondn = pndaspect * apondn
210 :
211 : else ! melt water runs off deformed ice
212 16036110 : apondn = c0
213 16036110 : hpondn = c0
214 : endif
215 1351150326 : apondn = max(apondn, c0)
216 :
217 : ! limit pond depth to maintain nonnegative freeboard
218 1351150326 : hpondn = min(hpondn, ((rhow-rhoi)*hi - rhos*hs)/rhofresh)
219 :
220 : ! fraction of grid cell covered by ponds
221 1351150326 : apondn = apondn * aicen
222 :
223 1351150326 : volpn = hpondn*apondn
224 1351150326 : if (volpn <= c0) then
225 758443097 : volpn = c0
226 758443097 : apondn = c0
227 758443097 : hpondn = c0
228 758443097 : 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 1351150326 : if (ktherm /= 2 .and. hpondn > c0 .and. dpscale > puny) then
238 11405871 : draft = (rhos*hs + rhoi*hi)/rhow + hpondn
239 11405871 : deltah = hpondn + hi - draft
240 11405871 : pressure_head = gravit * rhow * max(deltah, c0)
241 91246968 : Tmlt(:) = -sicen(:) * depressT
242 : call brine_permeability(qicen, &
243 11405871 : sicen, Tmlt, perm)
244 11405871 : if (icepack_warnings_aborted(subname)) return
245 11405871 : drain = perm*pressure_head*dt / (viscosity_dyn*hi) * dpscale
246 11405871 : deltah = min(drain, hpondn)
247 11405871 : dvn = -deltah*apondn
248 11405871 : volpn = volpn + dvn
249 : apondn = max(c0, min(apondn &
250 11405871 : + 0.5*dvn/(pndaspect*apondn), alvl_tmp*aicen))
251 11405871 : hpondn = c0
252 11405871 : if (apondn > puny) hpondn = volpn/apondn
253 : endif
254 :
255 : endif
256 :
257 : !-----------------------------------------------------------
258 : ! Reload tracer array
259 : !-----------------------------------------------------------
260 :
261 1352965794 : hpnd = hpondn
262 1352965794 : apnd = apondn / (aicen*alvl_tmp)
263 1352965794 : 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 11405871 : 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) ! LCOV_EXCL_LINE
279 : salin, & ! salinity (ppt) ! LCOV_EXCL_LINE
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 22811742 : Tin, & ! ice temperature (C) ! LCOV_EXCL_LINE
292 22811742 : 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 91246968 : do k = 1,nilyr
303 91246968 : Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt(k))
304 : enddo
305 :
306 : !-----------------------------------------------------------------
307 : ! brine salinity and liquid fraction
308 : !-----------------------------------------------------------------
309 :
310 91246968 : do k = 1,nilyr
311 79841097 : Sbr = c1/(1.e-3_dbl_kind - depressT/Tin(k)) ! Notz thesis eq 3.6
312 79841097 : phi(k) = salin(k)/Sbr ! liquid fraction
313 91246968 : if (phi(k) < 0.05) phi(k) = c0 ! impermeable
314 : enddo
315 :
316 : !-----------------------------------------------------------------
317 : ! permeability
318 : !-----------------------------------------------------------------
319 :
320 102652839 : perm = 3.0e-8_dbl_kind * (minval(phi))**3
321 :
322 11405871 : end subroutine brine_permeability
323 :
324 : !=======================================================================
325 :
326 : end module icepack_meltpond_lvl
327 :
328 : !=======================================================================
|
