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