Line data Source code
1 : module icepack_mushy_physics
2 :
3 : use icepack_kinds
4 : use icepack_parameters, only: c0, c1, c2, c4, c1000
5 : use icepack_parameters, only: puny
6 : use icepack_parameters, only: rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh
7 : use icepack_parameters, only: ksno
8 : use icepack_warnings, only: warnstr, icepack_warnings_add
9 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
10 :
11 : implicit none
12 :
13 : private
14 : public :: &
15 : conductivity_mush_array, & ! LCOV_EXCL_LINE
16 : conductivity_snow_array, & ! LCOV_EXCL_LINE
17 : icepack_enthalpy_snow, & ! LCOV_EXCL_LINE
18 : enthalpy_brine, & ! LCOV_EXCL_LINE
19 : icepack_enthalpy_mush, & ! LCOV_EXCL_LINE
20 : enthalpy_mush_liquid_fraction, & ! LCOV_EXCL_LINE
21 : enthalpy_of_melting, & ! LCOV_EXCL_LINE
22 : temperature_snow, & ! LCOV_EXCL_LINE
23 : icepack_mushy_temperature_mush, & ! LCOV_EXCL_LINE
24 : temperature_mush_liquid_fraction, & ! LCOV_EXCL_LINE
25 : liquidus_brine_salinity_mush, & ! LCOV_EXCL_LINE
26 : liquidus_temperature_mush, & ! LCOV_EXCL_LINE
27 : icepack_mushy_liquid_fraction, & ! LCOV_EXCL_LINE
28 : icepack_mushy_density_brine
29 :
30 : !-----------------------------------------------------------------
31 : ! Constants for Liquidus relation from Assur (1958)
32 : !-----------------------------------------------------------------
33 :
34 : ! liquidus relation - higher temperature region
35 : real(kind=dbl_kind), parameter :: &
36 : az1_liq = -18.48_dbl_kind, & ! LCOV_EXCL_LINE
37 : bz1_liq = 0.0_dbl_kind
38 :
39 : ! liquidus relation - lower temperature region
40 : real(kind=dbl_kind), parameter :: &
41 : az2_liq = -10.3085_dbl_kind, & ! LCOV_EXCL_LINE
42 : bz2_liq = 62.4_dbl_kind
43 :
44 : ! liquidus break
45 : real(kind=dbl_kind), parameter :: &
46 : Tb_liq = -7.6362968855167352_dbl_kind, & ! temperature of liquidus break ! LCOV_EXCL_LINE
47 : Sb_liq = 123.66702800276086_dbl_kind ! salinity of liquidus break
48 :
49 : ! basic liquidus relation constants
50 : real(kind=dbl_kind), parameter :: &
51 : az1p_liq = az1_liq / c1000, & ! LCOV_EXCL_LINE
52 : bz1p_liq = bz1_liq / c1000, & ! LCOV_EXCL_LINE
53 : az2p_liq = az2_liq / c1000, & ! LCOV_EXCL_LINE
54 : bz2p_liq = bz2_liq / c1000
55 :
56 : !-----------------------------------------------------------------
57 : ! Other parameters
58 : !-----------------------------------------------------------------
59 :
60 : real(kind=dbl_kind), parameter :: &
61 : ki = 2.3_dbl_kind , & ! fresh ice conductivity (W m-1 K-1) ! LCOV_EXCL_LINE
62 : kb = 0.5375_dbl_kind ! brine conductivity (W m-1 K-1)
63 :
64 : !=======================================================================
65 :
66 : contains
67 :
68 : !=======================================================================
69 : ! Physical Quantities
70 : !=======================================================================
71 : ! Detemine the conductivity of the mush from enthalpy and salinity
72 :
73 33852099 : subroutine conductivity_mush_array(nilyr, zqin, zSin, km)
74 :
75 : integer (kind=int_kind), intent(in) :: &
76 : nilyr ! number of ice layers
77 :
78 : real(kind=dbl_kind), dimension(:), intent(in) :: &
79 : zqin, & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
80 : zSin ! ice layer bulk salinity (ppt)
81 :
82 : real(kind=dbl_kind), dimension(:), intent(out) :: &
83 : km ! ice layer conductivity (W m-1 K-1)
84 :
85 : integer(kind=int_kind) :: &
86 : k ! ice layer index
87 :
88 3086234 : real(kind=dbl_kind) :: Tmush
89 :
90 : character(len=*),parameter :: subname='(conductivity_mush_array)'
91 :
92 270816792 : do k = 1, nilyr
93 :
94 236964693 : Tmush = icepack_mushy_temperature_mush(zqin(k), zSin(k))
95 :
96 270816792 : km(k) = heat_conductivity(Tmush, zSin(k))
97 :
98 : enddo ! k
99 :
100 33852099 : end subroutine conductivity_mush_array
101 :
102 : !=======================================================================
103 : !autodocument_start icepack_mushy_density_brine
104 : ! Compute density of brine from brine salinity
105 :
106 670043277 : function icepack_mushy_density_brine(Sbr) result(rho)
107 :
108 : real(kind=dbl_kind), intent(in) :: &
109 : Sbr ! brine salinity (ppt)
110 :
111 : real(kind=dbl_kind) :: &
112 : rho ! brine density (kg m-3)
113 :
114 : !autodocument_end
115 :
116 : real(kind=dbl_kind), parameter :: &
117 : a = 1000.3_dbl_kind , & ! zeroth empirical coefficient ! LCOV_EXCL_LINE
118 : b = 0.78237_dbl_kind , & ! linear empirical coefficient ! LCOV_EXCL_LINE
119 : c = 2.8008e-4_dbl_kind ! quadratic empirical coefficient
120 :
121 : character(len=*),parameter :: subname='(icepack_mushy_density_brine)'
122 :
123 670043277 : rho = a + b * Sbr + c * Sbr**2
124 :
125 670043277 : end function icepack_mushy_density_brine
126 :
127 : !=======================================================================
128 : ! Snow
129 : !=======================================================================
130 : ! Heat conductivity of the snow
131 :
132 20282724 : subroutine conductivity_snow_array(ks)
133 :
134 : real(kind=dbl_kind), dimension(:), intent(out) :: &
135 : ks ! snow layer conductivity (W m-1 K-1)
136 :
137 : character(len=*),parameter :: subname='(conductivity_snow_array)'
138 :
139 40565448 : ks = ksno
140 :
141 20282724 : end subroutine conductivity_snow_array
142 :
143 : !=======================================================================
144 : !autodocument_start icepack_enthalpy_snow
145 : ! Enthalpy of snow from snow temperature
146 :
147 44178434 : function icepack_enthalpy_snow(zTsn) result(zqsn)
148 :
149 : real(kind=dbl_kind), intent(in) :: &
150 : zTsn ! snow layer temperature (C)
151 :
152 : real(kind=dbl_kind) :: &
153 : zqsn ! snow layer enthalpy (J m-3)
154 :
155 : !autodocument_end
156 :
157 : character(len=*),parameter :: subname='(icepack_enthalpy_snow)'
158 :
159 44178434 : zqsn = -rhos * (-cp_ice * zTsn + Lfresh)
160 :
161 44178434 : end function icepack_enthalpy_snow
162 :
163 : !=======================================================================
164 : ! Temperature of snow from the snow enthalpy
165 :
166 41292873 : function temperature_snow(zqsn) result(zTsn)
167 :
168 : real(kind=dbl_kind), intent(in) :: &
169 : zqsn ! snow layer enthalpy (J m-3)
170 :
171 : real(kind=dbl_kind) :: &
172 : zTsn, & ! snow layer temperature (C) ! LCOV_EXCL_LINE
173 6153078 : A, B
174 :
175 : character(len=*),parameter :: subname='(temperature_snow)'
176 :
177 41292873 : A = c1 / (rhos * cp_ice)
178 41292873 : B = Lfresh / cp_ice
179 41292873 : zTsn = A * zqsn + B
180 :
181 41292873 : end function temperature_snow
182 :
183 : !=======================================================================
184 : ! Mushy Layer Formulation - Assur (1958) liquidus
185 : !=======================================================================
186 : ! Liquidus relation: equilibrium brine salinity as function of temperature
187 : ! based on empirical data from Assur (1958)
188 :
189 2197029654 : function liquidus_brine_salinity_mush(zTin) result(Sbr)
190 :
191 : real(kind=dbl_kind), intent(in) :: &
192 : zTin ! ice layer temperature (C)
193 :
194 : real(kind=dbl_kind) :: &
195 : Sbr ! ice brine salinity (ppt)
196 :
197 : real(kind=dbl_kind) :: &
198 : t_high , & ! mask for high temperature liquidus region ! LCOV_EXCL_LINE
199 211721018 : lsubzero ! mask for sub-zero temperatures
200 :
201 : real(kind=dbl_kind) :: &
202 : J1_liq, K1_liq, L1_liq, & ! temperature to brine salinity ! LCOV_EXCL_LINE
203 211721018 : J2_liq, K2_liq, L2_liq
204 :
205 : character(len=*),parameter :: subname='(liquidus_brine_salinty_mush)'
206 :
207 : ! temperature to brine salinity
208 2197029654 : J1_liq = bz1_liq / az1_liq
209 2197029654 : K1_liq = c1 / c1000
210 2197029654 : L1_liq = (c1 + bz1p_liq) / az1_liq
211 2197029654 : J2_liq = bz2_liq / az2_liq
212 2197029654 : K2_liq = c1 / c1000
213 2197029654 : L2_liq = (c1 + bz2p_liq) / az2_liq
214 :
215 2197029654 : t_high = merge(c1, c0, (zTin > Tb_liq))
216 2197029654 : lsubzero = merge(c1, c0, (zTin <= c0))
217 :
218 : Sbr = ((zTin + J1_liq) / (K1_liq * zTin + L1_liq)) * t_high + &
219 2197029654 : ((zTin + J2_liq) / (K2_liq * zTin + L2_liq)) * (c1 - t_high)
220 :
221 2197029654 : Sbr = Sbr * lsubzero
222 :
223 2197029654 : end function liquidus_brine_salinity_mush
224 :
225 : !=======================================================================
226 : ! Liquidus relation: equilibrium temperature as function of brine salinity
227 : ! based on empirical data from Assur (1958)
228 :
229 532337388 : function liquidus_temperature_mush(Sbr) result(zTin)
230 :
231 : real(kind=dbl_kind), intent(in) :: &
232 : Sbr ! ice brine salinity (ppt)
233 :
234 : real(kind=dbl_kind) :: &
235 : zTin ! ice layer temperature (C)
236 :
237 : real(kind=dbl_kind) :: &
238 49531582 : t_high ! mask for high temperature liquidus region
239 :
240 : real(kind=dbl_kind) :: &
241 : M1_liq, &! brine salinity to temperature ! LCOV_EXCL_LINE
242 : N1_liq, & ! LCOV_EXCL_LINE
243 : O1_liq, & ! LCOV_EXCL_LINE
244 : M2_liq, & ! LCOV_EXCL_LINE
245 : N2_liq, & ! LCOV_EXCL_LINE
246 49531582 : O2_liq
247 :
248 : character(len=*),parameter :: subname='(liquidus_temperature_mush)'
249 :
250 : ! brine salinity to temperature
251 532337388 : M1_liq = az1_liq
252 532337388 : N1_liq = -az1p_liq
253 532337388 : O1_liq = -bz1_liq / az1_liq
254 532337388 : M2_liq = az2_liq
255 532337388 : N2_liq = -az2p_liq
256 532337388 : O2_liq = -bz2_liq / az2_liq
257 :
258 532337388 : t_high = merge(c1, c0, (Sbr <= Sb_liq))
259 :
260 : zTin = ((Sbr / (M1_liq + N1_liq * Sbr)) + O1_liq) * t_high + &
261 532337388 : ((Sbr / (M2_liq + N2_liq * Sbr)) + O2_liq) * (c1 - t_high)
262 :
263 532337388 : end function liquidus_temperature_mush
264 :
265 : !=======================================================================
266 : !autodocument_start icepack_enthalpy_mush
267 : ! Enthalpy of mush from mush temperature and bulk salinity
268 :
269 45684843 : function icepack_enthalpy_mush(zTin, zSin) result(zqin)
270 :
271 : real(kind=dbl_kind), intent(in) :: &
272 : zTin, & ! ice layer temperature (C) ! LCOV_EXCL_LINE
273 : zSin ! ice layer bulk salinity (ppt)
274 :
275 : real(kind=dbl_kind) :: &
276 : zqin ! ice layer enthalpy (J m-3)
277 :
278 : !autodocument_end
279 :
280 : real(kind=dbl_kind) :: &
281 5968394 : phi ! ice liquid fraction
282 :
283 : character(len=*),parameter :: subname='(icepack_enthalpy_mush)'
284 :
285 45684843 : phi = icepack_mushy_liquid_fraction(zTin, zSin)
286 :
287 : zqin = phi * (cp_ocn * rhow - cp_ice * rhoi) * zTin + &
288 45684843 : rhoi * cp_ice * zTin - (c1 - phi) * rhoi * Lfresh
289 :
290 45684843 : end function icepack_enthalpy_mush
291 :
292 : !=======================================================================
293 : ! Enthalpy of mush from mush temperature and liquid fraction
294 :
295 507168144 : function enthalpy_mush_liquid_fraction(zTin, phi) result(zqin)
296 :
297 : real(kind=dbl_kind), intent(in) :: &
298 : zTin, & ! ice layer temperature (C) ! LCOV_EXCL_LINE
299 : phi ! liquid fraction
300 :
301 : real(kind=dbl_kind) :: &
302 : zqin ! ice layer enthalpy (J m-3)
303 :
304 : character(len=*),parameter :: subname='(enthalpy_mush_liquid_fraction)'
305 :
306 : zqin = phi * (cp_ocn * rhow - cp_ice * rhoi) * zTin + &
307 507168144 : rhoi * cp_ice * zTin - (c1 - phi) * rhoi * Lfresh
308 :
309 507168144 : end function enthalpy_mush_liquid_fraction
310 :
311 : !=======================================================================
312 : ! Enthalpy of melting of mush from bulk salinity
313 : ! Energy needed to fully melt mush (T < 0)
314 :
315 270816792 : function enthalpy_of_melting(zSin) result(qm)
316 :
317 : real(kind=dbl_kind), intent(in) :: &
318 : zSin ! ice layer bulk salinity (ppt)
319 :
320 : real(kind=dbl_kind) :: &
321 : qm ! melting ice enthalpy (J m-3)
322 :
323 : character(len=*),parameter :: subname='(enthalpy_of_melting)'
324 :
325 270816792 : qm = cp_ocn * rhow * liquidus_temperature_mush(zSin)
326 :
327 270816792 : end function enthalpy_of_melting
328 :
329 : !=======================================================================
330 : ! Enthalpy of brine (fully liquid) from temperature
331 :
332 1291987279 : function enthalpy_brine(zTin) result(qbr)
333 :
334 : real(kind=dbl_kind), intent(in) :: &
335 : zTin ! ice layer temperature (C)
336 :
337 : real(kind=dbl_kind) :: &
338 : qbr ! brine enthalpy (J m-3)
339 :
340 : character(len=*),parameter :: subname='(enthalpy_brine)'
341 :
342 1291987279 : qbr = cp_ocn * rhow * zTin
343 :
344 1291987279 : end function enthalpy_brine
345 :
346 : !=======================================================================
347 : !autodocument_start icepack_mushy_temperature_mush
348 : ! Temperature of mush from mush enthalpy and bulk salinity
349 :
350 954079630 : function icepack_mushy_temperature_mush(zqin, zSin) result(zTin)
351 :
352 : real(kind=dbl_kind), intent(in) :: &
353 : zqin , & ! ice enthalpy (J m-3) ! LCOV_EXCL_LINE
354 : zSin ! ice layer bulk salinity (ppt)
355 :
356 : real(kind=dbl_kind) :: &
357 : zTin ! ice layer temperature (C)
358 :
359 : !autodocument_end
360 :
361 : real(kind=dbl_kind) :: &
362 : qb , & ! liquidus break enthalpy ! LCOV_EXCL_LINE
363 : q0 , & ! fully melted enthalpy ! LCOV_EXCL_LINE
364 : A , & ! quadratic equation A parameter ! LCOV_EXCL_LINE
365 : B , & ! quadratic equation B parameter ! LCOV_EXCL_LINE
366 : C , & ! quadratic equation C parameter ! LCOV_EXCL_LINE
367 : S_low , & ! mask for salinity less than the liquidus break salinity ! LCOV_EXCL_LINE
368 : t_high , & ! mask for high temperature liquidus region ! LCOV_EXCL_LINE
369 : t_low , & ! mask for low temperature liquidus region ! LCOV_EXCL_LINE
370 90187846 : q_melt ! mask for all mush melted
371 :
372 : ! quadratic constants - higher temperature region
373 : real(kind=dbl_kind) :: &
374 : AS1_liq, AC1_liq, & ! quadratic constants - higher temperature region ! LCOV_EXCL_LINE
375 : BS1_liq, BC1_liq, BQ1_liq, & ! " ! LCOV_EXCL_LINE
376 : CS1_liq, CC1_liq, CQ1_liq, & ! " ! LCOV_EXCL_LINE
377 : AS2_liq, AC2_liq, & ! quadratic constants - lower temperature region ! LCOV_EXCL_LINE
378 : BS2_liq, BC2_liq, BQ2_liq, & ! " ! LCOV_EXCL_LINE
379 : CS2_liq, CC2_liq, CQ2_liq, & ! " ! LCOV_EXCL_LINE
380 : D_liq, E_liq, & ! break enthalpy constants ! LCOV_EXCL_LINE
381 : F1_liq, G1_liq, H1_liq, & ! just fully melted enthapy constants ! LCOV_EXCL_LINE
382 : F2_liq, G2_liq, H2_liq, & ! " ! LCOV_EXCL_LINE
383 90187846 : I_liq ! warmer than fully melted constants
384 :
385 : character(len=*),parameter :: subname='(icepack_mushy_temperature_mush)'
386 :
387 : !--------------------------------------------------------
388 :
389 : ! quadratic constants - higher temperature region
390 954079630 : AS1_liq = az1p_liq * (rhow * cp_ocn - rhoi * cp_ice)
391 954079630 : AC1_liq = rhoi * cp_ice * az1_liq
392 : BS1_liq = (c1 + bz1p_liq) * (rhow * cp_ocn - rhoi * cp_ice) &
393 954079630 : + rhoi * Lfresh * az1p_liq
394 954079630 : BQ1_liq = -az1_liq
395 954079630 : BC1_liq = rhoi * cp_ice * bz1_liq - rhoi * Lfresh * az1_liq
396 954079630 : CS1_liq = rhoi * Lfresh * (c1 + bz1p_liq)
397 954079630 : CQ1_liq = -bz1_liq
398 954079630 : CC1_liq = -rhoi * Lfresh * bz1_liq
399 :
400 : ! quadratic constants - lower temperature region
401 954079630 : AS2_liq = az2p_liq * (rhow * cp_ocn - rhoi * cp_ice)
402 954079630 : AC2_liq = rhoi * cp_ice * az2_liq
403 : BS2_liq = (c1 + bz2p_liq) * (rhow * cp_ocn - rhoi * cp_ice) &
404 954079630 : + rhoi * Lfresh * az2p_liq
405 954079630 : BQ2_liq = -az2_liq
406 954079630 : BC2_liq = rhoi * cp_ice * bz2_liq - rhoi * Lfresh * az2_liq
407 954079630 : CS2_liq = rhoi * Lfresh * (c1 + bz2p_liq)
408 954079630 : CQ2_liq = -bz2_liq
409 954079630 : CC2_liq = -rhoi * Lfresh * bz2_liq
410 :
411 : ! break enthalpy constants
412 : D_liq = ((c1 + az1p_liq*Tb_liq + bz1p_liq) &
413 : / ( az1_liq*Tb_liq + bz1_liq)) & ! LCOV_EXCL_LINE
414 954079630 : * ((cp_ocn*rhow - cp_ice*rhoi)*Tb_liq + Lfresh*rhoi)
415 954079630 : E_liq = cp_ice*rhoi*Tb_liq - Lfresh*rhoi
416 :
417 : ! just fully melted enthapy constants
418 954079630 : F1_liq = ( -c1000 * cp_ocn * rhow) / az1_liq
419 954079630 : G1_liq = -c1000
420 954079630 : H1_liq = (-bz1_liq * cp_ocn * rhow) / az1_liq
421 954079630 : F2_liq = ( -c1000 * cp_ocn * rhow) / az2_liq
422 954079630 : G2_liq = -c1000
423 954079630 : H2_liq = (-bz2_liq * cp_ocn * rhow) / az2_liq
424 :
425 : ! warmer than fully melted constants
426 954079630 : I_liq = c1 / (cp_ocn * rhow)
427 :
428 : ! just melted enthalpy
429 954079630 : S_low = merge(c1, c0, (zSin < Sb_liq))
430 : q0 = ((F1_liq * zSin) / (G1_liq + zSin) + H1_liq) * S_low + &
431 954079630 : ((F2_liq * zSin) / (G2_liq + zSin) + H2_liq) * (c1 - S_low)
432 954079630 : q_melt = merge(c1, c0, (zqin > q0))
433 :
434 : ! break enthalpy
435 954079630 : qb = D_liq * zSin + E_liq
436 954079630 : t_high = merge(c1, c0, (zqin > qb))
437 954079630 : t_low = c1 - t_high
438 :
439 : ! quadratic values
440 : A = (AS1_liq * zSin + AC1_liq) * t_high + &
441 954079630 : (AS2_liq * zSin + AC2_liq) * t_low
442 :
443 : B = (BS1_liq * zSin + BQ1_liq * zqin + BC1_liq) * t_high + &
444 954079630 : (BS2_liq * zSin + BQ2_liq * zqin + BC2_liq) * t_low
445 :
446 : C = (CS1_liq * zSin + CQ1_liq * zqin + CC1_liq) * t_high + &
447 954079630 : (CS2_liq * zSin + CQ2_liq * zqin + CC2_liq) * t_low
448 :
449 954079630 : zTin = (-B + sqrt(max(B**2 - c4 * A * C,puny))) / (c2 * A)
450 :
451 : ! change T if all melted
452 954079630 : zTin = q_melt * zqin * I_liq + (c1 - q_melt) * zTin
453 :
454 954079630 : end function icepack_mushy_temperature_mush
455 :
456 : !=======================================================================
457 : ! Temperature of mush from mush enthalpy and liquid fraction
458 :
459 236964693 : function temperature_mush_liquid_fraction(zqin, phi) result(zTin)
460 :
461 : real(kind=dbl_kind), intent(in) :: &
462 : zqin , & ! ice enthalpy (J m-3) ! LCOV_EXCL_LINE
463 : phi ! liquid fraction
464 :
465 : real(kind=dbl_kind) :: &
466 : zTin ! ice layer temperature (C)
467 :
468 : character(len=*),parameter :: subname='(temperature_mush_liquid_fraction)'
469 :
470 : zTin = (zqin + (c1 - phi) * rhoi * Lfresh) / &
471 236964693 : (phi * (cp_ocn * rhow - cp_ice * rhoi) + rhoi * cp_ice)
472 :
473 236964693 : end function temperature_mush_liquid_fraction
474 :
475 : !=======================================================================
476 : ! Mush heat conductivity from mush temperature and bulk salinity
477 :
478 236964693 : function heat_conductivity(zTin, zSin) result(km)
479 :
480 : real(kind=dbl_kind), intent(in) :: &
481 : zTin , & ! ice layer temperature (C) ! LCOV_EXCL_LINE
482 : zSin ! ice layer bulk salinity (ppt)
483 :
484 : real(kind=dbl_kind) :: &
485 : km ! ice layer conductivity (W m-1 K-1)
486 :
487 : real(kind=dbl_kind) :: &
488 21603638 : phi ! liquid fraction
489 :
490 : character(len=*),parameter :: subname='(heat_conductivity)'
491 :
492 236964693 : phi = icepack_mushy_liquid_fraction(zTin, zSin)
493 :
494 236964693 : km = phi * (kb - ki) + ki
495 :
496 236964693 : end function heat_conductivity
497 :
498 : !=======================================================================
499 : !autodocument_start icepack_mushy_liquid_fraction
500 : ! Liquid fraction of mush from mush temperature and bulk salinity
501 :
502 999764473 : function icepack_mushy_liquid_fraction(zTin, zSin) result(phi)
503 :
504 : real(kind=dbl_kind), intent(in) :: &
505 : zTin, & ! ice layer temperature (C) ! LCOV_EXCL_LINE
506 : zSin ! ice layer bulk salinity (ppt)
507 :
508 : real(kind=dbl_kind) :: &
509 : phi ! liquid fraction
510 :
511 : !autodocument_end
512 :
513 : real(kind=dbl_kind) :: &
514 96156240 : Sbr ! brine salinity (ppt)
515 :
516 : character(len=*),parameter :: subname='(icepack_mushy_liquid_fraction)'
517 :
518 999764473 : Sbr = max(liquidus_brine_salinity_mush(zTin),puny)
519 999764473 : phi = zSin / max(Sbr, zSin)
520 :
521 999764473 : end function icepack_mushy_liquid_fraction
522 :
523 : !=======================================================================
524 :
525 : end module icepack_mushy_physics
|
