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