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 158016 : subroutine compute_ponds_topo(dt, &
43 : ktherm, &
44 0 : aice, aicen, &
45 0 : vice, vicen, &
46 158016 : vsno, vsnon, &
47 : meltt, &
48 : fsurf, fpond, &
49 0 : Tsfcn, Tf, &
50 158016 : qicen, sicen, &
51 316032 : 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
61 : vsno, & ! total snow volume (m)
62 : Tf ! ocean freezing temperature [= ice bottom temperature] (degC)
63 :
64 : real (kind=dbl_kind), intent(inout) :: &
65 : vice, & ! total ice volume (m)
66 : fpond ! fresh water flux to ponds (m)
67 :
68 : real (kind=dbl_kind), dimension (:), intent(in) :: &
69 : aicen, & ! ice area fraction, per category
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, &
80 : sicen
81 :
82 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
83 : apnd, &
84 : hpnd, &
85 : ipnd
86 :
87 : real (kind=dbl_kind), intent(in) :: &
88 : meltt, & ! total surface meltwater flux
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 316032 : volpn, & ! pond volume per unit area, per category (m)
95 316032 : vuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
96 :
97 : real (kind=dbl_kind), dimension (ncat) :: &
98 158016 : apondn,& ! pond area fraction, per category
99 158016 : 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)
106 : dHui, & ! change in thickness of ice lid (m)
107 : omega, & ! conduction
108 : dTice, & ! temperature difference across ice lid (C)
109 : dvice, & ! change in ice volume (m)
110 : Tavg, & ! mean surface temperature across categories (C)
111 : Tp, & ! pond freezing temperature (C)
112 : rhoi_L,& ! (J/m^3)
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)
119 : Td = p15 , & ! temperature difference for freeze-up (C)
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 158016 : volp = c0
129 158016 : rhoi_L = Lfresh * rhoi ! (J/m^3)
130 :
131 948096 : do n = 1, ncat
132 : ! load tracers
133 : volp = volp + hpnd(n) &
134 790080 : * apnd(n) * aicen(n)
135 : vuin (n) = ipnd(n) &
136 790080 : * apnd(n) * aicen(n)
137 :
138 790080 : hpondn(n) = c0 ! pond depth, per category
139 948096 : 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 158016 : Tp = Timelt - Td
148 :
149 : !-----------------------------------------------------------------
150 : ! Identify grid cells with ponds
151 : !-----------------------------------------------------------------
152 :
153 158016 : hi = c0
154 158016 : if (aice > puny) hi = vice/aice
155 158016 : 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, &
163 : aice, vice, vsno, &
164 : aicen, vicen, vsnon, &
165 : qicen, sicen, &
166 : volpn, volp, &
167 26635 : apondn, hpondn, dvn )
168 26635 : if (icepack_warnings_aborted(subname)) return
169 :
170 26635 : fpond = fpond - dvn
171 :
172 : ! mean surface temperature
173 26635 : Tavg = c0
174 159810 : do n = 1, ncat
175 159810 : Tavg = Tavg + Tsfcn(n)*aicen(n)
176 : enddo
177 26635 : Tavg = Tavg / aice
178 :
179 133175 : do n = 1, ncat-1
180 :
181 133175 : 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 4547 : if (Tsfcn(n) > Tp) then
188 :
189 1555 : dvice = min(meltt*apondn(n), vuin(n))
190 1555 : if (dvice > puny) then
191 1499 : vuin (n) = vuin (n) - dvice
192 1499 : volpn(n) = volpn(n) + dvice
193 1499 : volp = volp + dvice
194 1499 : fpond = fpond + dvice
195 :
196 1499 : if (vuin(n) < puny .and. volpn(n) > puny) then
197 : ! ice lid melted and category is pond covered
198 610 : volpn(n) = volpn(n) + vuin(n)
199 610 : fpond = fpond + vuin(n)
200 610 : vuin(n) = c0
201 : endif
202 1499 : hpondn(n) = volpn(n) / apondn(n)
203 : endif
204 :
205 : !----------------------------------------------------------------
206 : ! freezing: existing upper ice layer grows
207 : !----------------------------------------------------------------
208 :
209 2992 : else if (volpn(n) > puny) then ! Tsfcn(i,j,n) <= Tp
210 :
211 : ! differential growth of base of surface floating ice layer
212 2992 : dTice = max(-Tsfcn(n)-Td, c0) ! > 0
213 2992 : omega = kice*DTice/rhoi_L
214 : dHui = sqrt(c2*omega*dt + (vuin(n)/aicen(n))**2) &
215 2992 : - vuin(n)/aicen(n)
216 :
217 2992 : dvice = min(dHui*apondn(n), volpn(n))
218 2992 : if (dvice > puny) then
219 2990 : vuin (n) = vuin (n) + dvice
220 2990 : volpn(n) = volpn(n) - dvice
221 2990 : volp = volp - dvice
222 2990 : fpond = fpond - dvice
223 2990 : 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 101993 : dHui = max(-fsurf*dt/rhoi_L, c0)
239 101993 : dvice = min(dHui*apondn(n), volpn(n))
240 101993 : if (dvice > puny) then
241 619 : vuin (n) = dvice
242 619 : volpn(n) = volpn(n) - dvice
243 619 : volp = volp - dvice
244 619 : fpond = fpond - dvice
245 619 : 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 131381 : fpond = fpond - volp
254 788286 : volpn(:) = c0
255 788286 : vuin (:) = c0
256 131381 : 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 948096 : do n = 1, ncat
265 : if (aicen(n) > puny .and. volpn(n) < puny &
266 790080 : .and. vuin (n) > puny) then
267 3 : vuin(n) = c0
268 : endif
269 :
270 : ! reload tracers
271 790080 : if (apondn(n) > puny) then
272 44104 : ipnd(n) = vuin(n) / apondn(n)
273 : else
274 745976 : vuin(n) = c0
275 745976 : ipnd(n) = c0
276 : endif
277 948096 : if (aicen(n) > puny) then
278 452058 : apnd(n) = apondn(n) / aicen(n)
279 452058 : hpnd(n) = hpondn(n)
280 : else
281 338022 : apnd(n) = c0
282 338022 : hpnd(n) = c0
283 338022 : 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 26635 : subroutine pond_area(dt, &
294 : ktherm, &
295 : aice, vice, vsno, &
296 26635 : aicen, vicen, vsnon,&
297 26635 : qicen, sicen, &
298 26635 : volpn, volp, &
299 26635 : 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, &
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, &
327 : m_index, &
328 : permflag
329 :
330 : real (kind=dbl_kind), dimension(ncat) :: &
331 53270 : hicen, &
332 53270 : hsnon, &
333 53270 : asnon, &
334 26635 : alfan, &
335 53270 : betan, &
336 53270 : cum_max_vol, &
337 53270 : reduced_aicen
338 :
339 : real (kind=dbl_kind), dimension(0:ncat) :: &
340 26635 : cum_max_vol_tmp
341 :
342 : real (kind=dbl_kind) :: &
343 : hpond, &
344 : drain, &
345 : floe_weight, &
346 : pressure_head, &
347 : hsl_rel, &
348 : deltah, &
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 159810 : do n = 1, ncat
381 :
382 133175 : apondn(n) = c0
383 133175 : hpondn(n) = c0
384 :
385 133175 : if (aicen(n) < puny) then
386 9505 : hicen(n) = c0
387 9505 : hsnon(n) = c0
388 9505 : reduced_aicen(n) = c0
389 9505 : asnon(n) = c0
390 : else
391 123670 : hicen(n) = vicen(n) / aicen(n)
392 123670 : hsnon(n) = vsnon(n) / aicen(n)
393 123670 : reduced_aicen(n) = c1 ! n=ncat
394 123670 : if (n < ncat) reduced_aicen(n) = aicen(n) &
395 97035 : * max(0.2_dbl_kind,(-0.024_dbl_kind*hicen(n) + 0.832_dbl_kind))
396 123670 : 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 133175 : alfan(n) = p6 * hicen(n)
409 133175 : betan(n) = p4 * hicen(n)
410 :
411 133175 : cum_max_vol(n) = c0
412 159810 : cum_max_vol_tmp(n) = c0
413 :
414 : enddo ! ncat
415 :
416 26635 : cum_max_vol_tmp(0) = c0
417 26635 : drain = c0
418 26635 : dvolp = c0
419 :
420 : !--------------------------------------------------------------------------
421 : ! the maximum amount of water that can be contained up to each ice category
422 : !--------------------------------------------------------------------------
423 :
424 133175 : do n = 1, ncat-1 ! last category can not hold any volume
425 :
426 106540 : 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 360402 : (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))
431 :
432 :
433 : ! subtract snow solid volumes from lower categories in current level
434 360402 : 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
437 : asnon(ns) * & ! area of snow from that category
438 360402 : 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 5260 : cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
444 : endif
445 133175 : 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 26635 : cum_max_vol_tmp(ncat) = cum_max_vol_tmp(ncat-1) ! last category holds no volume
452 159810 : 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 26635 : if (volp >= cum_max_vol(ncat)) then
458 0 : drain = volp - cum_max_vol(ncat) + puny
459 0 : volp = volp - drain
460 0 : dvolp = drain
461 0 : if (volp < puny) then
462 0 : dvolp = dvolp + volp
463 0 : 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 26635 : alfan, volp, cum_max_vol, hpond, m_index)
471 26635 : if (icepack_warnings_aborted(subname)) return
472 :
473 82372 : do n=1, m_index
474 55737 : hpondn(n) = max((hpond - alfan(n) + alfan(1)), c0)
475 82372 : 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 26635 : floe_weight = (vsno*rhos + rhoi*vice + rhow*volp) / aice
484 : hsl_rel = floe_weight / rhow &
485 159810 : - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
486 :
487 26635 : deltah = hpond - hsl_rel
488 26635 : pressure_head = gravit * rhow * max(deltah, c0)
489 :
490 : ! drain if ice is permeable
491 26635 : permflag = 0
492 26635 : if (ktherm /= 2 .and. pressure_head > c0) then
493 17970 : do n = 1, ncat-1
494 17970 : if (hicen(n) > c0) then
495 5668 : call permeability_phi(qicen(:,n), sicen(:,n), perm)
496 5668 : if (icepack_warnings_aborted(subname)) return
497 5668 : if (perm > c0) permflag = 1
498 5668 : drain = perm*apondn(n)*pressure_head*dt / (viscosity_dyn*hicen(n))
499 5668 : dvolp = dvolp + min(drain, volp)
500 5668 : volp = max(volp - drain, c0)
501 5668 : if (volp < puny) then
502 2048 : dvolp = dvolp + volp
503 2048 : volp = c0
504 : endif
505 : endif
506 : enddo
507 :
508 : ! adjust melt pond dimensions
509 3594 : if (permflag > 0) then
510 : ! recompute pond depth
511 : call calc_hpond(reduced_aicen, asnon, hsnon, &
512 3594 : alfan, volp, cum_max_vol, hpond, m_index)
513 3594 : if (icepack_warnings_aborted(subname)) return
514 8232 : do n=1, m_index
515 4638 : hpondn(n) = hpond - alfan(n) + alfan(1)
516 8232 : 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 49593 : do n=1,m_index-1
528 : volpn(n) = apondn(n) * hpondn(n) &
529 49593 : - (rhos/rhow) * asnon(n) * min(hsnon(n), hpondn(n))
530 : enddo
531 :
532 : ! Calculate pond volume for highest category = remaining pond volume
533 26635 : if (m_index == 1) volpn(m_index) = volp
534 26635 : if (m_index > 1) then
535 40247 : if (volp > sum(volpn(1:m_index-1))) then
536 40247 : 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)) + &
545 0 : volp
546 : endif
547 : endif
548 :
549 74180 : do n=1,m_index
550 74180 : if (apondn(n) > puny) then
551 44104 : hpondn(n) = volpn(n) / apondn(n)
552 : else
553 3441 : dvolp = dvolp + volpn(n)
554 3441 : hpondn(n) = c0
555 3441 : volpn(n) = c0
556 3441 : apondn(n) = c0
557 : end if
558 : enddo
559 112265 : do n = m_index+1, ncat
560 85630 : hpondn(n) = c0
561 85630 : apondn(n) = c0
562 112265 : volpn (n) = c0
563 : enddo
564 :
565 : end subroutine pond_area
566 :
567 : !=======================================================================
568 :
569 60458 : subroutine calc_hpond(aicen, asnon, hsnon, &
570 60458 : alfan, volp, cum_max_vol, hpond, m_index)
571 :
572 : real (kind=dbl_kind), dimension(:), intent(in) :: &
573 : aicen, &
574 : asnon, &
575 : hsnon, &
576 : alfan, &
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 60458 : hitl, &
592 60458 : aicetl
593 :
594 : real (kind=dbl_kind) :: &
595 : rem_vol, &
596 : area, &
597 : vol, &
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 30229 : if (volp < puny) then
607 2048 : hpond = c0
608 2048 : 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 28181 : m_index = 0 ! 1:m_index categories have water in them
631 60375 : do n = 1, ncat
632 60375 : if (volp <= cum_max_vol(n)) then
633 28181 : m_index = n
634 28181 : if (n == 1) then
635 7298 : rem_vol = volp
636 : else
637 20883 : rem_vol = volp - cum_max_vol(n-1)
638 : endif
639 28181 : exit ! to break out of the loop
640 : endif
641 : enddo
642 28181 : 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 225448 : hitl(:) = c0
661 225448 : aicetl(:) = c0
662 88556 : do n = 1, m_index
663 : hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
664 60375 : alfan(m_index+1) - alfan(m_index)), c0)
665 60375 : aicetl(n) = asnon(n)
666 :
667 88556 : aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
668 : enddo
669 28181 : hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
670 28181 : aicetl(m_index+1) = c0
671 :
672 : !----------------------------------------------------------------
673 : ! reorder array according to hitl
674 : ! snow heights not necessarily in height order
675 : !----------------------------------------------------------------
676 :
677 116737 : do ns = 1, m_index+1
678 311262 : do n = 0, m_index - ns + 1
679 283081 : 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 88534 : do n = 1, m_index+1
708 :
709 : area = sum(aicetl(:)) - & ! total area of sub-layer
710 1133529 : (rhos/rhow) * sum(aicetl(n:ncat+1)) ! area of sub-layer occupied by snow
711 :
712 88534 : vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
713 :
714 88534 : if (vol >= rem_vol) then ! have reached the sub-layer with the depth within
715 28181 : hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
716 28181 : exit
717 : else ! still in sub-layer below the sub-layer with the depth
718 60353 : rem_vol = rem_vol - vol
719 : endif
720 :
721 : enddo
722 :
723 : endif
724 :
725 30229 : end subroutine calc_hpond
726 :
727 : !=======================================================================
728 :
729 : ! determine the liquid fraction of brine in the ice and the permeability
730 :
731 5668 : 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)
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
744 : Sbr ! brine salinity
745 :
746 : real (kind=dbl_kind), dimension(nilyr) :: &
747 11336 : Tin, & ! ice temperature
748 5668 : 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 45344 : do k = 1,nilyr
760 39676 : Tmlt = -sicen(k) * depressT
761 45344 : Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt)
762 : enddo
763 :
764 : !-----------------------------------------------------------------
765 : ! brine salinity and liquid fraction
766 : !-----------------------------------------------------------------
767 :
768 51012 : if (maxval(Tin) <= -c2) then
769 :
770 : ! Assur 1958
771 0 : do k = 1,nilyr
772 : Sbr = - 1.2_dbl_kind &
773 : -21.8_dbl_kind * Tin(k) &
774 : - 0.919_dbl_kind * Tin(k)**2 &
775 0 : - 0.01878_dbl_kind * Tin(k)**3
776 0 : phi(k) = sicen(k)/Sbr ! liquid fraction
777 : enddo ! k
778 :
779 : else
780 :
781 : ! Notz 2005 thesis eq. 3.2
782 45344 : do k = 1,nilyr
783 : Sbr = -17.6_dbl_kind * Tin(k) &
784 : - 0.389_dbl_kind * Tin(k)**2 &
785 39676 : - 0.00362_dbl_kind* Tin(k)**3
786 39676 : 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 45344 : phi(k) = sicen(k) / Sbr ! liquid fraction
792 : enddo
793 :
794 : endif
795 :
796 : !-----------------------------------------------------------------
797 : ! permeability
798 : !-----------------------------------------------------------------
799 :
800 51012 : 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 : !=======================================================================
|
