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