Line data Source code
1 : !=======================================================================
2 :
3 : ! Melt pond evolution based on the ice topography as inferred from
4 : ! the ice thickness distribution. This code is based on (but differs
5 : ! from) that described in
6 : !
7 : ! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond
8 : ! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi:
9 : ! 10.1029/2006JC003836.
10 : !
11 : ! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of a
12 : ! physically based melt pond scheme into the sea ice component of a
13 : ! climate model. J. Geophys. Res. 115, C08012, doi: 10.1029/2009JC005568.
14 : !
15 : ! authors Daniela Flocco (UCL)
16 : ! Adrian Turner (UCL)
17 : ! 2010 ECH added module based on original code from Daniela Flocco, UCL
18 : ! 2012 DSCHR modifications
19 :
20 : module icepack_meltpond_topo
21 :
22 : use icepack_kinds
23 : use icepack_parameters, only: c0, c1, c2, p01, p1, p15, p4, p6
24 : use icepack_parameters, only: puny, viscosity_dyn, rhoi, rhos, rhow, Timelt, Lfresh
25 : use icepack_parameters, only: gravit, depressT, kice, ice_ref_salinity
26 : use icepack_tracers, only: ncat, nilyr
27 : use icepack_warnings, only: warnstr, icepack_warnings_add
28 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
29 : use icepack_therm_shared, only: calculate_Tin_from_qin
30 :
31 : implicit none
32 :
33 : private
34 : public :: compute_ponds_topo
35 :
36 : !=======================================================================
37 :
38 : contains
39 :
40 : !=======================================================================
41 :
42 17485104 : subroutine compute_ponds_topo(dt, &
43 : ktherm, & ! LCOV_EXCL_LINE
44 0 : aice, aicen, & ! LCOV_EXCL_LINE
45 0 : vice, vicen, & ! LCOV_EXCL_LINE
46 17485104 : vsno, vsnon, & ! LCOV_EXCL_LINE
47 : meltt, & ! LCOV_EXCL_LINE
48 : fsurf, fpond, & ! LCOV_EXCL_LINE
49 0 : Tsfcn, Tf, & ! LCOV_EXCL_LINE
50 17485104 : qicen, sicen, & ! LCOV_EXCL_LINE
51 34970208 : apnd, hpnd, ipnd )
52 :
53 : integer (kind=int_kind), intent(in) :: &
54 : ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
55 :
56 : real (kind=dbl_kind), intent(in) :: &
57 : dt ! time step (s)
58 :
59 : real (kind=dbl_kind), intent(in) :: &
60 : aice, & ! total ice area fraction ! LCOV_EXCL_LINE
61 : vsno, & ! total snow volume (m) ! LCOV_EXCL_LINE
62 : Tf ! ocean freezing temperature [= ice bottom temperature] (degC)
63 :
64 : real (kind=dbl_kind), intent(inout) :: &
65 : vice, & ! total ice volume (m) ! LCOV_EXCL_LINE
66 : fpond ! fresh water flux to ponds (m)
67 :
68 : real (kind=dbl_kind), dimension (:), intent(in) :: &
69 : aicen, & ! ice area fraction, per category ! LCOV_EXCL_LINE
70 : vsnon ! snow volume, per category (m)
71 :
72 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
73 : vicen ! ice volume, per category (m)
74 :
75 : real (kind=dbl_kind), dimension (:), intent(in) :: &
76 : Tsfcn
77 :
78 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
79 : qicen, & ! LCOV_EXCL_LINE
80 : sicen
81 :
82 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
83 : apnd, & ! LCOV_EXCL_LINE
84 : hpnd, & ! LCOV_EXCL_LINE
85 : ipnd
86 :
87 : real (kind=dbl_kind), intent(in) :: &
88 : meltt, & ! total surface meltwater flux ! LCOV_EXCL_LINE
89 : fsurf ! thermodynamic heat flux at ice/snow surface (W/m^2)
90 :
91 : ! local variables
92 :
93 : real (kind=dbl_kind), dimension (ncat) :: &
94 34970208 : volpn, & ! pond volume per unit area, per category (m) ! LCOV_EXCL_LINE
95 34970208 : vuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
96 :
97 : real (kind=dbl_kind), dimension (ncat) :: &
98 17485104 : apondn,& ! pond area fraction, per category ! LCOV_EXCL_LINE
99 17485104 : hpondn ! pond depth, per category (m)
100 :
101 : real (kind=dbl_kind) :: &
102 : volp ! total volume of pond, per unit area of pond (m)
103 :
104 : real (kind=dbl_kind) :: &
105 : hi, & ! ice thickness (m) ! LCOV_EXCL_LINE
106 : dHui, & ! change in thickness of ice lid (m) ! LCOV_EXCL_LINE
107 : omega, & ! conduction ! LCOV_EXCL_LINE
108 : dTice, & ! temperature difference across ice lid (C) ! LCOV_EXCL_LINE
109 : dvice, & ! change in ice volume (m) ! LCOV_EXCL_LINE
110 : Tavg, & ! mean surface temperature across categories (C) ! LCOV_EXCL_LINE
111 : Tp, & ! pond freezing temperature (C) ! LCOV_EXCL_LINE
112 : rhoi_L,& ! (J/m^3) ! LCOV_EXCL_LINE
113 : dvn ! change in melt pond volume for fresh water budget
114 :
115 : integer (kind=int_kind) :: n ! loop indices
116 :
117 : real (kind=dbl_kind), parameter :: &
118 : hicemin = p1 , & ! minimum ice thickness with ponds (m) ! LCOV_EXCL_LINE
119 : Td = p15 , & ! temperature difference for freeze-up (C) ! LCOV_EXCL_LINE
120 : min_volp = 1.e-4_dbl_kind ! minimum pond volume (m)
121 :
122 : character(len=*),parameter :: subname='(compute_ponds_topo)'
123 :
124 : !---------------------------------------------------------------
125 : ! initialize
126 : !---------------------------------------------------------------
127 :
128 17485104 : volp = c0
129 17485104 : rhoi_L = Lfresh * rhoi ! (J/m^3)
130 :
131 122395728 : do n = 1, ncat
132 : ! load tracers
133 : volp = volp + hpnd(n) &
134 104910624 : * apnd(n) * aicen(n)
135 : vuin (n) = ipnd(n) &
136 104910624 : * apnd(n) * aicen(n)
137 :
138 104910624 : hpondn(n) = c0 ! pond depth, per category
139 122395728 : apondn(n) = c0 ! pond area, per category
140 : enddo
141 :
142 : ! The freezing temperature for meltponds is assumed slightly below 0C,
143 : ! as if meltponds had a little salt in them. The salt budget is not
144 : ! altered for meltponds, but if it were then an actual pond freezing
145 : ! temperature could be computed.
146 :
147 17485104 : Tp = Timelt - Td
148 :
149 : !-----------------------------------------------------------------
150 : ! Identify grid cells with ponds
151 : !-----------------------------------------------------------------
152 :
153 17485104 : hi = c0
154 17485104 : if (aice > puny) hi = vice/aice
155 17485104 : if ( aice > p01 .and. hi > hicemin .and. &
156 : volp > min_volp*aice) then
157 :
158 : !--------------------------------------------------------------
159 : ! calculate pond area and depth
160 : !--------------------------------------------------------------
161 : call pond_area(dt, &
162 : ktherm, & ! LCOV_EXCL_LINE
163 : aice, vice, vsno, & ! LCOV_EXCL_LINE
164 : aicen, vicen, vsnon, & ! LCOV_EXCL_LINE
165 : qicen, sicen, & ! LCOV_EXCL_LINE
166 : volpn, volp, & ! LCOV_EXCL_LINE
167 561128 : apondn, hpondn, dvn )
168 561128 : if (icepack_warnings_aborted(subname)) return
169 :
170 561128 : fpond = fpond - dvn
171 :
172 : ! mean surface temperature
173 561128 : Tavg = c0
174 3927896 : do n = 1, ncat
175 3927896 : Tavg = Tavg + Tsfcn(n)*aicen(n)
176 : enddo
177 561128 : Tavg = Tavg / aice
178 :
179 3366768 : do n = 1, ncat-1
180 :
181 3366768 : if (vuin(n) > puny) then
182 :
183 : !----------------------------------------------------------------
184 : ! melting: floating upper ice layer melts in whole or part
185 : !----------------------------------------------------------------
186 : ! Use Tsfc for each category
187 1227168 : if (Tsfcn(n) > Tp) then
188 :
189 0 : dvice = min(meltt*apondn(n), vuin(n))
190 0 : if (dvice > puny) then
191 0 : vuin (n) = vuin (n) - dvice
192 0 : volpn(n) = volpn(n) + dvice
193 0 : volp = volp + dvice
194 0 : fpond = fpond + dvice
195 :
196 0 : if (vuin(n) < puny .and. volpn(n) > puny) then
197 : ! ice lid melted and category is pond covered
198 0 : volpn(n) = volpn(n) + vuin(n)
199 0 : fpond = fpond + vuin(n)
200 0 : vuin(n) = c0
201 : endif
202 0 : hpondn(n) = volpn(n) / apondn(n)
203 : endif
204 :
205 : !----------------------------------------------------------------
206 : ! freezing: existing upper ice layer grows
207 : !----------------------------------------------------------------
208 :
209 1227168 : else if (volpn(n) > puny) then ! Tsfcn(i,j,n) <= Tp
210 :
211 : ! differential growth of base of surface floating ice layer
212 980519 : dTice = max(-Tsfcn(n)-Td, c0) ! > 0
213 980519 : omega = kice*DTice/rhoi_L
214 : dHui = sqrt(c2*omega*dt + (vuin(n)/aicen(n))**2) &
215 980519 : - vuin(n)/aicen(n)
216 :
217 980519 : dvice = min(dHui*apondn(n), volpn(n))
218 980519 : if (dvice > puny) then
219 980519 : vuin (n) = vuin (n) + dvice
220 980519 : volpn(n) = volpn(n) - dvice
221 980519 : volp = volp - dvice
222 980519 : fpond = fpond - dvice
223 980519 : hpondn(n) = volpn(n) / apondn(n)
224 : endif
225 :
226 : endif ! Tsfcn(i,j,n)
227 :
228 : !----------------------------------------------------------------
229 : ! freezing: upper ice layer begins to form
230 : ! note: albedo does not change
231 : !----------------------------------------------------------------
232 : else ! vuin < puny
233 :
234 : ! thickness of newly formed ice
235 : ! the surface temperature of a meltpond is the same as that
236 : ! of the ice underneath (0C), and the thermodynamic surface
237 : ! flux is the same
238 1578472 : dHui = max(-fsurf*dt/rhoi_L, c0)
239 1578472 : dvice = min(dHui*apondn(n), volpn(n))
240 1578472 : if (dvice > puny) then
241 17294 : vuin (n) = dvice
242 17294 : volpn(n) = volpn(n) - dvice
243 17294 : volp = volp - dvice
244 17294 : fpond = fpond - dvice
245 17294 : hpondn(n)= volpn(n) / apondn(n)
246 : endif
247 :
248 : endif ! vuin
249 :
250 : enddo ! ncat
251 :
252 : else ! remove ponds on thin ice
253 16923976 : fpond = fpond - volp
254 118467832 : volpn(:) = c0
255 118467832 : vuin (:) = c0
256 16923976 : volp = c0
257 : endif
258 :
259 : !---------------------------------------------------------------
260 : ! remove ice lid if there is no liquid pond
261 : ! vuin may be nonzero on category ncat due to dynamics
262 : !---------------------------------------------------------------
263 :
264 122395728 : do n = 1, ncat
265 : if (aicen(n) > puny .and. volpn(n) < puny &
266 104910624 : .and. vuin (n) > puny) then
267 600859 : vuin(n) = c0
268 : endif
269 :
270 : ! reload tracers
271 104910624 : if (apondn(n) > puny) then
272 997813 : ipnd(n) = vuin(n) / apondn(n)
273 : else
274 103912811 : vuin(n) = c0
275 103912811 : ipnd(n) = c0
276 : endif
277 122395728 : if (aicen(n) > puny) then
278 24093747 : apnd(n) = apondn(n) / aicen(n)
279 24093747 : hpnd(n) = hpondn(n)
280 : else
281 80816877 : apnd(n) = c0
282 80816877 : hpnd(n) = c0
283 80816877 : ipnd(n) = c0
284 : endif
285 : enddo ! n
286 :
287 : end subroutine compute_ponds_topo
288 :
289 : !=======================================================================
290 :
291 : ! Computes melt pond area, pond depth and melting rates
292 :
293 561128 : subroutine pond_area(dt, &
294 : ktherm, & ! LCOV_EXCL_LINE
295 : aice, vice, vsno, & ! LCOV_EXCL_LINE
296 561128 : aicen, vicen, vsnon,& ! LCOV_EXCL_LINE
297 561128 : qicen, sicen, & ! LCOV_EXCL_LINE
298 561128 : volpn, volp, & ! LCOV_EXCL_LINE
299 561128 : apondn,hpondn,dvolp )
300 :
301 : integer (kind=int_kind), intent(in) :: &
302 : ktherm ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
303 :
304 : real (kind=dbl_kind), intent(in) :: &
305 : dt, aice, vice, vsno
306 :
307 : real (kind=dbl_kind), dimension(:), intent(in) :: &
308 : aicen, vicen, vsnon
309 :
310 : real (kind=dbl_kind), dimension(:,:), intent(in) :: &
311 : qicen, & ! LCOV_EXCL_LINE
312 : sicen
313 :
314 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
315 : volpn
316 :
317 : real (kind=dbl_kind), intent(inout) :: &
318 : volp, dvolp
319 :
320 : real (kind=dbl_kind), dimension(:), intent(out) :: &
321 : apondn, hpondn
322 :
323 : ! local variables
324 :
325 : integer (kind=int_kind) :: &
326 : n, ns, & ! LCOV_EXCL_LINE
327 : m_index, & ! LCOV_EXCL_LINE
328 : permflag
329 :
330 : real (kind=dbl_kind), dimension(ncat) :: &
331 1122256 : hicen, & ! LCOV_EXCL_LINE
332 1122256 : hsnon, & ! LCOV_EXCL_LINE
333 1122256 : asnon, & ! LCOV_EXCL_LINE
334 561128 : alfan, & ! LCOV_EXCL_LINE
335 1122256 : betan, & ! LCOV_EXCL_LINE
336 1122256 : cum_max_vol, & ! LCOV_EXCL_LINE
337 1122256 : reduced_aicen
338 :
339 : real (kind=dbl_kind), dimension(0:ncat) :: &
340 561128 : cum_max_vol_tmp
341 :
342 : real (kind=dbl_kind) :: &
343 : hpond, & ! LCOV_EXCL_LINE
344 : drain, & ! LCOV_EXCL_LINE
345 : floe_weight, & ! LCOV_EXCL_LINE
346 : pressure_head, & ! LCOV_EXCL_LINE
347 : hsl_rel, & ! LCOV_EXCL_LINE
348 : deltah, & ! LCOV_EXCL_LINE
349 : perm
350 :
351 : character(len=*),parameter :: subname='(pond_area)'
352 :
353 : !-----------|
354 : ! |
355 : ! |-----------|
356 : !___________|___________|______________________________________sea-level
357 : ! | |
358 : ! | |---^--------|
359 : ! | | | |
360 : ! | | | |-----------| |-------
361 : ! | | |alfan(n)| | |
362 : ! | | | | |--------------|
363 : ! | | | | | |
364 : !---------------------------v-------------------------------------------
365 : ! | | ^ | | |
366 : ! | | | | |--------------|
367 : ! | | |betan(n)| | |
368 : ! | | | |-----------| |-------
369 : ! | | | |
370 : ! | |---v------- |
371 : ! | |
372 : ! |-----------|
373 : ! |
374 : !-----------|
375 :
376 : !-------------------------------------------------------------------
377 : ! initialize
378 : !-------------------------------------------------------------------
379 :
380 3927896 : do n = 1, ncat
381 :
382 3366768 : apondn(n) = c0
383 3366768 : hpondn(n) = c0
384 :
385 3366768 : if (aicen(n) < puny) then
386 5892 : hicen(n) = c0
387 5892 : hsnon(n) = c0
388 5892 : reduced_aicen(n) = c0
389 5892 : asnon(n) = c0
390 : else
391 3360876 : hicen(n) = vicen(n) / aicen(n)
392 3360876 : hsnon(n) = vsnon(n) / aicen(n)
393 3360876 : reduced_aicen(n) = c1 ! n=ncat
394 3360876 : if (n < ncat) reduced_aicen(n) = aicen(n) &
395 2801128 : * max(0.2_dbl_kind,(-0.024_dbl_kind*hicen(n) + 0.832_dbl_kind))
396 3360876 : asnon(n) = reduced_aicen(n)
397 : endif
398 :
399 : ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
400 : ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
401 : ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
402 : ! categories. alfa and beta partition the ITD - they are areas not thicknesses!
403 : ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
404 : ! Here, alfa = 60% of the ice area (and since hice is constant in a category,
405 : ! alfan = 60% of the ice volume) in each category lies above the reference line,
406 : ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
407 :
408 3366768 : alfan(n) = p6 * hicen(n)
409 3366768 : betan(n) = p4 * hicen(n)
410 :
411 3366768 : cum_max_vol(n) = c0
412 3927896 : cum_max_vol_tmp(n) = c0
413 :
414 : enddo ! ncat
415 :
416 561128 : cum_max_vol_tmp(0) = c0
417 561128 : drain = c0
418 561128 : dvolp = c0
419 :
420 : !--------------------------------------------------------------------------
421 : ! the maximum amount of water that can be contained up to each ice category
422 : !--------------------------------------------------------------------------
423 :
424 3366768 : do n = 1, ncat-1 ! last category can not hold any volume
425 :
426 2805640 : if (alfan(n+1) >= alfan(n) .and. alfan(n+1) > c0) then
427 :
428 : ! total volume in level including snow
429 : cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
430 11196976 : (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
431 :
432 :
433 : ! subtract snow solid volumes from lower categories in current level
434 11196976 : do ns = 1, n
435 : cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
436 : - rhos/rhow * & ! fraction of snow that is occupied by solid ! LCOV_EXCL_LINE
437 : asnon(ns) * & ! area of snow from that category ! LCOV_EXCL_LINE
438 11196976 : max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), c0)
439 : ! thickness of snow from ns layer in n layer
440 : enddo
441 :
442 : else ! assume higher categories unoccupied
443 5892 : cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
444 : endif
445 3366768 : if (cum_max_vol_tmp(n) < c0) then
446 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
447 0 : call icepack_warnings_add(subname//' topo ponds: negative melt pond volume')
448 0 : return
449 : endif
450 : enddo
451 561128 : cum_max_vol_tmp(ncat) = cum_max_vol_tmp(ncat-1) ! last category holds no volume
452 3927896 : cum_max_vol (1:ncat) = cum_max_vol_tmp(1:ncat)
453 :
454 : !----------------------------------------------------------------
455 : ! is there more meltwater than can be held in the floe?
456 : !----------------------------------------------------------------
457 561128 : if (volp >= cum_max_vol(ncat)) then
458 372 : drain = volp - cum_max_vol(ncat) + puny
459 372 : volp = volp - drain
460 372 : dvolp = drain
461 372 : if (volp < puny) then
462 372 : dvolp = dvolp + volp
463 372 : volp = c0
464 : endif
465 : endif
466 :
467 : ! height and area corresponding to the remaining volume
468 :
469 : call calc_hpond(reduced_aicen, asnon, hsnon, &
470 561128 : alfan, volp, cum_max_vol, hpond, m_index)
471 561128 : if (icepack_warnings_aborted(subname)) return
472 :
473 1558941 : do n=1, m_index
474 997813 : hpondn(n) = max((hpond - alfan(n) + alfan(1)), c0)
475 1558941 : apondn(n) = reduced_aicen(n)
476 : enddo
477 :
478 : !------------------------------------------------------------------------
479 : ! drainage due to ice permeability - Darcy's law
480 : !------------------------------------------------------------------------
481 :
482 : ! sea water level
483 561128 : floe_weight = (vsno*rhos + rhoi*vice + rhow*volp) / aice
484 : hsl_rel = floe_weight / rhow &
485 3927896 : - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
486 :
487 561128 : deltah = hpond - hsl_rel
488 561128 : pressure_head = gravit * rhow * max(deltah, c0)
489 :
490 : ! drain if ice is permeable
491 561128 : permflag = 0
492 561128 : if (ktherm /= 2 .and. pressure_head > c0) then
493 54252 : do n = 1, ncat-1
494 54252 : if (hicen(n) > c0) then
495 40770 : call permeability_phi(qicen(:,n), sicen(:,n), perm)
496 40770 : if (icepack_warnings_aborted(subname)) return
497 40770 : if (perm > c0) permflag = 1
498 40770 : drain = perm*apondn(n)*pressure_head*dt / (viscosity_dyn*hicen(n))
499 40770 : dvolp = dvolp + min(drain, volp)
500 40770 : volp = max(volp - drain, c0)
501 40770 : if (volp < puny) then
502 372 : dvolp = dvolp + volp
503 372 : volp = c0
504 : endif
505 : endif
506 : enddo
507 :
508 : ! adjust melt pond dimensions
509 9042 : if (permflag > 0) then
510 : ! recompute pond depth
511 : call calc_hpond(reduced_aicen, asnon, hsnon, &
512 9042 : alfan, volp, cum_max_vol, hpond, m_index)
513 9042 : if (icepack_warnings_aborted(subname)) return
514 17712 : do n=1, m_index
515 8670 : hpondn(n) = hpond - alfan(n) + alfan(1)
516 17712 : apondn(n) = reduced_aicen(n)
517 : enddo
518 : endif
519 : endif ! pressure_head
520 :
521 : !------------------------------------------------------------------------
522 : ! total melt pond volume in category does not include snow volume
523 : ! snow in melt ponds is not melted
524 : !------------------------------------------------------------------------
525 :
526 : ! Calculate pond volume for lower categories
527 998185 : do n=1,m_index-1
528 : volpn(n) = apondn(n) * hpondn(n) &
529 998185 : - (rhos/rhow) * asnon(n) * min(hsnon(n), hpondn(n))
530 : enddo
531 :
532 : ! Calculate pond volume for highest category = remaining pond volume
533 561128 : if (m_index == 1) volpn(m_index) = volp
534 561128 : if (m_index > 1) then
535 866230 : if (volp > sum(volpn(1:m_index-1))) then
536 866230 : volpn(m_index) = volp - sum(volpn(1:m_index-1))
537 : else
538 0 : volpn(m_index) = c0
539 0 : hpondn(m_index) = c0
540 0 : apondn(m_index) = c0
541 : ! If remaining pond volume is negative reduce pond volume of
542 : ! lower category
543 0 : if (volp+puny < sum(volpn(1:m_index-1))) &
544 : volpn(m_index-1) = volpn(m_index-1) - sum(volpn(1:m_index-1)) + & ! LCOV_EXCL_LINE
545 0 : volp
546 : endif
547 : endif
548 :
549 1558941 : do n=1,m_index
550 1558941 : if (apondn(n) > puny) then
551 997813 : hpondn(n) = volpn(n) / apondn(n)
552 : else
553 0 : dvolp = dvolp + volpn(n)
554 0 : hpondn(n) = c0
555 0 : volpn(n) = c0
556 0 : apondn(n) = c0
557 : end if
558 : enddo
559 2930083 : do n = m_index+1, ncat
560 2368955 : hpondn(n) = c0
561 2368955 : apondn(n) = c0
562 2930083 : volpn (n) = c0
563 : enddo
564 :
565 : end subroutine pond_area
566 :
567 : !=======================================================================
568 :
569 1140340 : subroutine calc_hpond(aicen, asnon, hsnon, &
570 1140340 : alfan, volp, cum_max_vol, hpond, m_index)
571 :
572 : real (kind=dbl_kind), dimension(:), intent(in) :: &
573 : aicen, & ! LCOV_EXCL_LINE
574 : asnon, & ! LCOV_EXCL_LINE
575 : hsnon, & ! LCOV_EXCL_LINE
576 : alfan, & ! LCOV_EXCL_LINE
577 : cum_max_vol
578 :
579 : real (kind=dbl_kind), intent(in) :: &
580 : volp
581 :
582 : real (kind=dbl_kind), intent(out) :: &
583 : hpond
584 :
585 : integer (kind=int_kind), intent(out) :: &
586 : m_index
587 :
588 : integer :: n, ns
589 :
590 : real (kind=dbl_kind), dimension(0:ncat+1) :: &
591 1140340 : hitl, & ! LCOV_EXCL_LINE
592 1140340 : aicetl
593 :
594 : real (kind=dbl_kind) :: &
595 : rem_vol, & ! LCOV_EXCL_LINE
596 : area, & ! LCOV_EXCL_LINE
597 : vol, & ! LCOV_EXCL_LINE
598 : tmp
599 :
600 : character(len=*),parameter :: subname='(calc_hpond)'
601 :
602 : !----------------------------------------------------------------
603 : ! hpond is zero if volp is zero - have we fully drained?
604 : !----------------------------------------------------------------
605 :
606 570170 : if (volp < puny) then
607 744 : hpond = c0
608 744 : m_index = 0
609 : else
610 :
611 : !----------------------------------------------------------------
612 : ! Calculate the category where water fills up to
613 : !----------------------------------------------------------------
614 :
615 : !----------|
616 : ! |
617 : ! |
618 : ! |----------| -- --
619 : !__________|__________|_________________________________________ ^
620 : ! | | rem_vol ^ | Semi-filled
621 : ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer
622 : ! | | | |
623 : ! | | | |hpond
624 : ! | | |----------| | |-------
625 : ! | | | | | |
626 : ! | | | |---v-----|
627 : ! | | m_index | | |
628 : !-------------------------------------------------------------
629 :
630 569426 : m_index = 0 ! 1:m_index categories have water in them
631 1006483 : do n = 1, ncat
632 1006483 : if (volp <= cum_max_vol(n)) then
633 569426 : m_index = n
634 569426 : if (n == 1) then
635 140253 : rem_vol = volp
636 : else
637 429173 : rem_vol = volp - cum_max_vol(n-1)
638 : endif
639 569426 : exit ! to break out of the loop
640 : endif
641 : enddo
642 569426 : m_index = min(ncat-1, m_index)
643 :
644 : !----------------------------------------------------------------
645 : ! semi-filled layer may have m_index different snows in it
646 : !----------------------------------------------------------------
647 :
648 : !----------------------------------------------------------- ^
649 : ! | alfan(m_index+1)
650 : ! |
651 : !hitl(3)--> |----------| |
652 : !hitl(2)--> |------------| * * * * *| |
653 : !hitl(1)--> |----------|* * * * * * |* * * * * | |
654 : !hitl(0)-->------------------------------------------------- | ^
655 : ! various snows from lower categories | |alfa(m_index)
656 :
657 : ! hitl - heights of the snow layers from thinner and current categories
658 : ! aicetl - area of each snow depth in this layer
659 :
660 5124834 : hitl(:) = c0
661 5124834 : aicetl(:) = c0
662 1575909 : do n = 1, m_index
663 : hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
664 1006483 : alfan(m_index+1) - alfan(m_index)), c0)
665 1006483 : aicetl(n) = asnon(n)
666 :
667 1575909 : aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
668 : enddo
669 569426 : hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
670 569426 : aicetl(m_index+1) = c0
671 :
672 : !----------------------------------------------------------------
673 : ! reorder array according to hitl
674 : ! snow heights not necessarily in height order
675 : !----------------------------------------------------------------
676 :
677 2145335 : do ns = 1, m_index+1
678 5172668 : do n = 0, m_index - ns + 1
679 4603242 : if (hitl(n) > hitl(n+1)) then ! swap order
680 0 : tmp = hitl(n)
681 0 : hitl(n) = hitl(n+1)
682 0 : hitl(n+1) = tmp
683 0 : tmp = aicetl(n)
684 0 : aicetl(n) = aicetl(n+1)
685 0 : aicetl(n+1) = tmp
686 : endif
687 : enddo
688 : enddo
689 :
690 : !----------------------------------------------------------------
691 : ! divide semi-filled layer into set of sublayers each vertically homogenous
692 : !----------------------------------------------------------------
693 :
694 : !hitl(3)----------------------------------------------------------------
695 : ! | * * * * * * * *
696 : ! |* * * * * * * * *
697 : !hitl(2)----------------------------------------------------------------
698 : ! | * * * * * * * * | * * * * * * * *
699 : ! |* * * * * * * * * |* * * * * * * * *
700 : !hitl(1)----------------------------------------------------------------
701 : ! | * * * * * * * * | * * * * * * * * | * * * * * * * *
702 : ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
703 : !hitl(0)----------------------------------------------------------------
704 : ! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
705 :
706 : ! move up over layers incrementing volume
707 1160654 : do n = 1, m_index+1
708 :
709 : area = sum(aicetl(:)) - & ! total area of sub-layer
710 17917742 : (rhos/rhow) * sum(aicetl(n:ncat+1)) ! area of sub-layer occupied by snow
711 :
712 1160654 : vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
713 :
714 1160654 : if (vol >= rem_vol) then ! have reached the sub-layer with the depth within
715 569426 : hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
716 569426 : exit
717 : else ! still in sub-layer below the sub-layer with the depth
718 591228 : rem_vol = rem_vol - vol
719 : endif
720 :
721 : enddo
722 :
723 : endif
724 :
725 570170 : end subroutine calc_hpond
726 :
727 : !=======================================================================
728 :
729 : ! determine the liquid fraction of brine in the ice and the permeability
730 :
731 40770 : subroutine permeability_phi(qicen, sicen, perm)
732 :
733 : real (kind=dbl_kind), dimension(:), intent(in) :: &
734 : qicen, & ! energy of melting for each ice layer (J/m2) ! LCOV_EXCL_LINE
735 : sicen ! salinity (ppt)
736 :
737 : real (kind=dbl_kind), intent(out) :: &
738 : perm ! permeability
739 :
740 : ! local variables
741 :
742 : real (kind=dbl_kind) :: &
743 : Tmlt, & ! melting temperature ! LCOV_EXCL_LINE
744 : Sbr ! brine salinity
745 :
746 : real (kind=dbl_kind), dimension(nilyr) :: &
747 81540 : Tin, & ! ice temperature ! LCOV_EXCL_LINE
748 40770 : phi ! liquid fraction
749 :
750 : integer (kind=int_kind) :: k
751 :
752 : character(len=*),parameter :: subname='(permeability_phi)'
753 :
754 : !-----------------------------------------------------------------
755 : ! Compute ice temperatures from enthalpies using quadratic formula
756 : ! NOTE this assumes Tmlt = Si * depressT
757 : !-----------------------------------------------------------------
758 :
759 326160 : do k = 1,nilyr
760 285390 : Tmlt = -sicen(k) * depressT
761 326160 : Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt)
762 : enddo
763 :
764 : !-----------------------------------------------------------------
765 : ! brine salinity and liquid fraction
766 : !-----------------------------------------------------------------
767 :
768 366930 : if (maxval(Tin) <= -c2) then
769 :
770 : ! Assur 1958
771 144912 : do k = 1,nilyr
772 : Sbr = - 1.2_dbl_kind &
773 : -21.8_dbl_kind * Tin(k) & ! LCOV_EXCL_LINE
774 : - 0.919_dbl_kind * Tin(k)**2 & ! LCOV_EXCL_LINE
775 126798 : - 0.01878_dbl_kind * Tin(k)**3
776 144912 : phi(k) = sicen(k)/Sbr ! liquid fraction
777 : enddo ! k
778 :
779 : else
780 :
781 : ! Notz 2005 thesis eq. 3.2
782 181248 : do k = 1,nilyr
783 : Sbr = -17.6_dbl_kind * Tin(k) &
784 : - 0.389_dbl_kind * Tin(k)**2 & ! LCOV_EXCL_LINE
785 158592 : - 0.00362_dbl_kind* Tin(k)**3
786 158592 : if (Sbr == c0) then
787 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
788 0 : call icepack_warnings_add(subname//' topo ponds: zero brine salinity in permeability')
789 0 : return
790 : endif
791 181248 : phi(k) = sicen(k) / Sbr ! liquid fraction
792 : enddo
793 :
794 : endif
795 :
796 : !-----------------------------------------------------------------
797 : ! permeability
798 : !-----------------------------------------------------------------
799 :
800 366930 : perm = 3.0e-08_dbl_kind * (minval(phi))**3
801 :
802 : end subroutine permeability_phi
803 :
804 : !=======================================================================
805 :
806 : end module icepack_meltpond_topo
807 :
808 : !=======================================================================
|
