Line data Source code
1 : !=========================================================================
2 : !
3 : ! Shared thermo variables, subroutines
4 : !
5 : ! authors: Elizabeth C. Hunke, LANL
6 :
7 : module icepack_therm_shared
8 :
9 : use icepack_kinds
10 :
11 : use icepack_parameters, only: c0, c1, c2, c4, p5, pi, puny, Tocnfrz
12 : use icepack_parameters, only: cp_ocn, cp_ice, rhoi, rhos, Tffresh, TTTice, qqqice
13 : use icepack_parameters, only: stefan_boltzmann, emissivity, Lfresh, Tsmelt
14 : use icepack_parameters, only: saltmax, min_salin, depressT
15 : use icepack_parameters, only: ktherm, tfrz_option
16 : use icepack_parameters, only: calc_Tsfc
17 : use icepack_tracers, only: nilyr, nslyr
18 : use icepack_warnings, only: warnstr, icepack_warnings_add
19 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
20 :
21 : use icepack_mushy_physics, only: icepack_enthalpy_mush
22 : use icepack_mushy_physics, only: temperature_snow
23 : use icepack_mushy_physics, only: icepack_mushy_temperature_mush
24 : use icepack_mushy_physics, only: liquidus_temperature_mush
25 :
26 : implicit none
27 :
28 : private
29 : public :: calculate_Tin_from_qin, &
30 : surface_heat_flux, &
31 : dsurface_heat_flux_dTsf, &
32 : icepack_init_salinity, &
33 : icepack_salinity_profile, &
34 : icepack_init_enthalpy, &
35 : icepack_ice_temperature, &
36 : icepack_snow_temperature, &
37 : icepack_liquidus_temperature, &
38 : icepack_sea_freezing_temperature, &
39 : adjust_enthalpy
40 :
41 : real (kind=dbl_kind), parameter, public :: &
42 : ferrmax = 1.0e-3_dbl_kind ! max allowed energy flux error (W m-2)
43 : ! recommend ferrmax < 0.01 W m-2
44 :
45 : real (kind=dbl_kind), parameter, public :: &
46 : Tmin = -100.0_dbl_kind ! min allowed internal temperature (deg C)
47 :
48 : logical (kind=log_kind), public :: &
49 : l_brine ! if true, treat brine pocket effects
50 :
51 : !=======================================================================
52 :
53 : contains
54 :
55 : !=======================================================================
56 : !
57 : ! Compute the internal ice temperatures from enthalpy using
58 : ! quadratic formula
59 :
60 11550422 : function calculate_Tin_from_qin (qin, Tmltk) &
61 : result(Tin)
62 :
63 : real (kind=dbl_kind), intent(in) :: &
64 : qin , & ! enthalpy
65 : Tmltk ! melting temperature at one level
66 :
67 : real (kind=dbl_kind) :: &
68 : Tin ! internal temperature
69 :
70 : ! local variables
71 :
72 : real (kind=dbl_kind) :: &
73 : aa1,bb1,cc1,csqrt ! quadratic solvers
74 :
75 : character(len=*),parameter :: subname='(calculate_Tin_from_qin)'
76 :
77 11550422 : if (l_brine) then
78 11550422 : aa1 = cp_ice
79 11550422 : bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh
80 11550422 : cc1 = Lfresh * Tmltk
81 11550422 : csqrt = bb1*bb1 - c4*aa1*cc1
82 11550422 : if (csqrt < c0) then
83 0 : call icepack_warnings_add(subname//' sqrt error: ')
84 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
85 0 : return
86 : endif
87 11550422 : Tin = min((-bb1 - sqrt(csqrt)) / (c2*aa1),Tmltk)
88 :
89 : else ! fresh ice
90 0 : Tin = (Lfresh + qin/rhoi) / cp_ice
91 : endif
92 :
93 11550422 : end function calculate_Tin_from_qin
94 :
95 : !=======================================================================
96 : ! Surface heat flux
97 : !=======================================================================
98 :
99 : ! heat flux into ice
100 :
101 26660276 : subroutine surface_heat_flux(Tsf, fswsfc, &
102 : rhoa, flw, &
103 : potT, Qa, &
104 : shcoef, lhcoef, &
105 : flwoutn, fsensn, &
106 : flatn, fsurfn)
107 :
108 : ! input surface temperature
109 : real(kind=dbl_kind), intent(in) :: &
110 : Tsf ! ice/snow surface temperature (C)
111 :
112 : ! input variables
113 : real(kind=dbl_kind), intent(in) :: &
114 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2)
115 : rhoa , & ! air density (kg/m^3)
116 : flw , & ! incoming longwave radiation (W/m^2)
117 : potT , & ! air potential temperature (K)
118 : Qa , & ! specific humidity (kg/kg)
119 : shcoef , & ! transfer coefficient for sensible heat
120 : lhcoef ! transfer coefficient for latent heat
121 :
122 : ! output
123 : real(kind=dbl_kind), intent(out) :: &
124 : fsensn , & ! surface downward sensible heat (W m-2)
125 : flatn , & ! surface downward latent heat (W m-2)
126 : flwoutn , & ! upward LW at surface (W m-2)
127 : fsurfn ! net flux to top surface, excluding fcondtopn
128 :
129 : ! local variables
130 : real(kind=dbl_kind) :: &
131 : TsfK , & ! ice/snow surface temperature (K)
132 : Qsfc , & ! saturated surface specific humidity (kg/kg)
133 : qsat , & ! the saturation humidity of air (kg/m^3)
134 : flwdabs , & ! downward longwave absorbed heat flx (W/m^2)
135 : tmpvar ! 1/TsfK
136 :
137 : character(len=*),parameter :: subname='(surface_heat_flux)'
138 :
139 : ! ice surface temperature in Kelvin
140 26660276 : TsfK = Tsf + Tffresh
141 : ! TsfK = max(Tsf + Tffresh, c1)
142 26660276 : tmpvar = c1/TsfK
143 :
144 : ! saturation humidity
145 26660276 : qsat = qqqice * exp(-TTTice*tmpvar)
146 26660276 : Qsfc = qsat / rhoa
147 :
148 : ! longwave radiative flux
149 26660276 : flwdabs = emissivity * flw
150 26660276 : flwoutn = -emissivity * stefan_boltzmann * TsfK**4
151 :
152 : ! downward latent and sensible heat fluxes
153 26660276 : fsensn = shcoef * (potT - TsfK)
154 26660276 : flatn = lhcoef * (Qa - Qsfc)
155 :
156 : ! combine fluxes
157 26660276 : fsurfn = fswsfc + flwdabs + flwoutn + fsensn + flatn
158 :
159 26660276 : end subroutine surface_heat_flux
160 :
161 : !=======================================================================
162 :
163 18969472 : subroutine dsurface_heat_flux_dTsf(Tsf, rhoa, &
164 : shcoef, lhcoef, &
165 : dfsurfn_dTsf, dflwoutn_dTsf, &
166 : dfsensn_dTsf, dflatn_dTsf)
167 :
168 : ! input surface temperature
169 : real(kind=dbl_kind), intent(in) :: &
170 : Tsf ! ice/snow surface temperature (C)
171 :
172 : ! input variables
173 : real(kind=dbl_kind), intent(in) :: &
174 : rhoa , & ! air density (kg/m^3)
175 : shcoef , & ! transfer coefficient for sensible heat
176 : lhcoef ! transfer coefficient for latent heat
177 :
178 : ! output
179 : real(kind=dbl_kind), intent(out) :: &
180 : dfsurfn_dTsf ! derivative of net flux to top surface, excluding fcondtopn
181 :
182 : real(kind=dbl_kind), intent(out) :: &
183 : dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature
184 : dfsensn_dTsf , & ! derivative of sensible heat flux wrt surface temperature
185 : dflatn_dTsf ! derivative of latent heat flux wrt surface temperature
186 :
187 : ! local variables
188 : real(kind=dbl_kind) :: &
189 : TsfK , & ! ice/snow surface temperature (K)
190 : dQsfc_dTsf , & ! saturated surface specific humidity (kg/kg)
191 : qsat , & ! the saturation humidity of air (kg/m^3)
192 : tmpvar ! 1/TsfK
193 :
194 : character(len=*),parameter :: subname='(dsurface_heat_flux_dTsf)'
195 :
196 : ! ice surface temperature in Kelvin
197 : ! TsfK = max(Tsf + Tffresh, c1)
198 18969472 : TsfK = Tsf + Tffresh
199 18969472 : tmpvar = c1/TsfK
200 :
201 : ! saturation humidity
202 18969472 : qsat = qqqice * exp(-TTTice*tmpvar)
203 18969472 : dQsfc_dTsf = TTTice * tmpvar * tmpvar * (qsat / rhoa)
204 :
205 : ! longwave radiative flux
206 18969472 : dflwoutn_dTsf = -emissivity * stefan_boltzmann * c4*TsfK**3
207 :
208 : ! downward latent and sensible heat fluxes
209 18969472 : dfsensn_dTsf = -shcoef
210 18969472 : dflatn_dTsf = -lhcoef * dQsfc_dTsf
211 :
212 : ! combine fluxes
213 18969472 : dfsurfn_dTsf = dflwoutn_dTsf + dfsensn_dTsf + dflatn_dTsf
214 :
215 18969472 : end subroutine dsurface_heat_flux_dTsf
216 :
217 : !=======================================================================
218 : !autodocument_start icepack_init_salinity
219 : ! Initialize the vertical profile of ice salinity.
220 : ! This subroutine was renamed from icepack_init_thermo in Oct 2024
221 : !
222 : ! authors: C. M. Bitz, UW
223 : ! William H. Lipscomb, LANL
224 :
225 83 : subroutine icepack_init_salinity(sprofile)
226 :
227 : real (kind=dbl_kind), dimension(:), intent(out) :: &
228 : sprofile ! vertical salinity profile
229 :
230 : !autodocument_end
231 :
232 : integer (kind=int_kind) :: k ! ice layer index
233 : real (kind=dbl_kind) :: zn ! normalized ice thickness
234 :
235 : character(len=*),parameter :: subname='(icepack_init_salinity)'
236 :
237 : !-----------------------------------------------------------------
238 : ! Determine l_brine based on saltmax.
239 : ! Thermodynamic solver will not converge if l_brine is true and
240 : ! saltmax is close to zero.
241 : ! Set l_brine to false for zero layer thermodynamics
242 : !-----------------------------------------------------------------
243 :
244 83 : l_brine = .false.
245 83 : if (saltmax > min_salin) l_brine = .true.
246 :
247 : !-----------------------------------------------------------------
248 : ! Prescibe vertical profile of salinity.
249 : ! Note this profile is only used for BL99 thermodynamics.
250 : !-----------------------------------------------------------------
251 :
252 83 : if (l_brine) then
253 628 : do k = 1, nilyr
254 : zn = (real(k,kind=dbl_kind)-p5) / &
255 545 : real(nilyr,kind=dbl_kind)
256 : ! sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
257 545 : sprofile(k)=icepack_salinity_profile(zn)
258 628 : sprofile(k) = max(sprofile(k), min_salin)
259 : enddo ! k
260 83 : sprofile(nilyr+1) = saltmax
261 :
262 : else ! .not. l_brine
263 0 : do k = 1, nilyr+1
264 0 : sprofile(k) = c0
265 : enddo
266 : endif ! l_brine
267 :
268 83 : end subroutine icepack_init_salinity
269 :
270 : !=======================================================================
271 : !autodocument_start icepack_salinity_profile
272 : ! Initial salinity profile
273 : !
274 : ! authors: C. M. Bitz, UW
275 : ! William H. Lipscomb, LANL
276 :
277 545 : function icepack_salinity_profile(zn) result(salinity)
278 :
279 : real(kind=dbl_kind), intent(in) :: &
280 : zn ! depth
281 :
282 : real(kind=dbl_kind) :: &
283 : salinity ! initial salinity profile
284 :
285 : !autodocument_end
286 :
287 : real (kind=dbl_kind), parameter :: &
288 : nsal = 0.407_dbl_kind, &
289 : msal = 0.573_dbl_kind
290 :
291 : character(len=*),parameter :: subname='(icepack_salinity_profile)'
292 :
293 545 : salinity = (saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
294 :
295 545 : end function icepack_salinity_profile
296 :
297 : !=======================================================================
298 : !autodocument_start icepack_init_enthalpy
299 : ! This subroutine was renamed from icepack_init_trcr in Oct 2024
300 : !
301 806 : subroutine icepack_init_enthalpy(Tair, Tf, &
302 806 : Sprofile, Tprofile, &
303 : Tsfc, &
304 1612 : qin, qsn)
305 :
306 : real (kind=dbl_kind), intent(in) :: &
307 : Tair, & ! air temperature (K)
308 : Tf ! freezing temperature (C)
309 :
310 : real (kind=dbl_kind), dimension(:), intent(in) :: &
311 : Sprofile, & ! vertical salinity profile (ppt)
312 : Tprofile ! vertical temperature profile (C)
313 :
314 : real (kind=dbl_kind), intent(out) :: &
315 : Tsfc ! surface temperature (C)
316 :
317 : real (kind=dbl_kind), dimension(:), intent(out) :: &
318 : qin, & ! ice enthalpy profile (J/m3)
319 : qsn ! snow enthalpy profile (J/m3)
320 :
321 : !autodocument_end
322 :
323 : ! local variables
324 :
325 : integer (kind=int_kind) :: k
326 :
327 : real (kind=dbl_kind) :: &
328 : slope, Ti
329 :
330 : character(len=*),parameter :: subname='(icepack_init_enthalpy)'
331 :
332 : ! surface temperature
333 806 : Tsfc = Tf ! default
334 806 : if (calc_Tsfc) Tsfc = min(Tsmelt, Tair - Tffresh) ! deg C
335 :
336 : ! ice enthalpy
337 6232 : do k = 1, nilyr
338 : ! assume linear temp profile and compute enthalpy
339 5426 : slope = Tf - Tsfc
340 : Ti = Tsfc + slope*(real(k,kind=dbl_kind)-p5) &
341 5426 : /real(nilyr,kind=dbl_kind)
342 6232 : if (ktherm == 2) then
343 4550 : qin(k) = icepack_enthalpy_mush(Ti, Sprofile(k))
344 : else
345 : qin(k) = -(rhoi * (cp_ice*(Tprofile(k)-Ti) &
346 876 : + Lfresh*(c1-Tprofile(k)/Ti) - cp_ocn*Tprofile(k)))
347 : endif
348 : enddo ! nilyr
349 :
350 : ! snow enthalpy
351 2212 : do k = 1, nslyr
352 1406 : Ti = min(c0, Tsfc)
353 2212 : qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
354 : enddo ! nslyr
355 :
356 806 : end subroutine icepack_init_enthalpy
357 :
358 : !=======================================================================
359 : !autodocument_start icepack_liquidus_temperature
360 : ! compute liquidus temperature
361 :
362 2147180 : function icepack_liquidus_temperature(Sin) result(Tmlt)
363 :
364 : real(dbl_kind), intent(in) :: Sin
365 : real(dbl_kind) :: Tmlt
366 :
367 : !autodocument_end
368 :
369 : character(len=*),parameter :: subname='(icepack_liquidus_temperature)'
370 :
371 2147180 : if (ktherm == 2) then
372 :
373 2050580 : Tmlt = liquidus_temperature_mush(Sin)
374 :
375 : else
376 :
377 96600 : Tmlt = -depressT * Sin
378 :
379 : endif
380 :
381 2147180 : end function icepack_liquidus_temperature
382 :
383 : !=======================================================================
384 : !autodocument_start icepack_sea_freezing_temperature
385 : ! compute ocean freezing temperature
386 :
387 2629488 : function icepack_sea_freezing_temperature(sss) result(Tf)
388 :
389 : real(dbl_kind), intent(in) :: sss
390 : real(dbl_kind) :: Tf
391 :
392 : !autodocument_end
393 :
394 : character(len=*),parameter :: subname='(icepack_sea_freezing_temperature)'
395 :
396 2629488 : if (trim(tfrz_option) == 'mushy') then
397 :
398 2146848 : Tf = icepack_liquidus_temperature(sss) ! deg C
399 :
400 482640 : elseif (trim(tfrz_option) == 'linear_salt') then
401 :
402 386112 : Tf = -depressT * sss ! deg C
403 :
404 96528 : elseif (trim(tfrz_option) == 'constant') then
405 :
406 0 : Tf = Tocnfrz
407 :
408 96528 : elseif (trim(tfrz_option) == 'minus1p8') then
409 :
410 96528 : Tf = -1.8_dbl_kind
411 :
412 : else
413 :
414 0 : call icepack_warnings_add(subname//' tfrz_option unsupported: '//trim(tfrz_option))
415 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
416 0 : return
417 :
418 : endif
419 :
420 2629488 : end function icepack_sea_freezing_temperature
421 :
422 : !=======================================================================
423 : !autodocument_start icepack_ice_temperature
424 : ! compute ice temperature
425 :
426 885331 : function icepack_ice_temperature(qin, Sin) result(Tin)
427 :
428 : real(kind=dbl_kind), intent(in) :: qin, Sin
429 : real(kind=dbl_kind) :: Tin
430 :
431 : !autodocument_end
432 :
433 : real(kind=dbl_kind) :: Tmlts
434 :
435 : character(len=*),parameter :: subname='(icepack_ice_temperature)'
436 :
437 885331 : if (ktherm == 2) then
438 :
439 885331 : Tin = icepack_mushy_temperature_mush(qin, Sin)
440 :
441 : else
442 :
443 0 : Tmlts = -depressT * Sin
444 0 : Tin = calculate_Tin_from_qin(qin,Tmlts)
445 :
446 : endif
447 :
448 885331 : end function icepack_ice_temperature
449 :
450 : !=======================================================================
451 : !autodocument_start icepack_snow_temperature
452 : ! compute snow temperature
453 :
454 0 : function icepack_snow_temperature(qin) result(Tsn)
455 :
456 : real(kind=dbl_kind), intent(in) :: qin
457 : real(kind=dbl_kind) :: Tsn
458 :
459 : !autodocument_end
460 :
461 : character(len=*),parameter :: subname='(icepack_snow_temperature)'
462 :
463 0 : if (ktherm == 2) then
464 :
465 0 : Tsn = temperature_snow(qin)
466 :
467 : else
468 :
469 0 : Tsn = (Lfresh + qin/rhos)/cp_ice
470 :
471 : endif
472 :
473 0 : end function icepack_snow_temperature
474 :
475 : !=======================================================================
476 : !
477 : ! Conserving energy, compute the new enthalpy of equal-thickness ice
478 : ! or snow layers.
479 : !
480 : ! authors William H. Lipscomb, LANL
481 : ! C. M. Bitz, UW
482 :
483 20764175 : subroutine adjust_enthalpy (nlyr, &
484 20764175 : z1, z2, &
485 : hlyr, hn, &
486 20764175 : qn)
487 :
488 : integer (kind=int_kind), intent(in) :: &
489 : nlyr ! number of layers (nilyr or nslyr)
490 :
491 : real (kind=dbl_kind), dimension (:), intent(in) :: &
492 : z1 , & ! interface depth for old, unequal layers (m)
493 : z2 ! interface depth for new, equal layers (m)
494 :
495 : real (kind=dbl_kind), intent(in) :: &
496 : hlyr ! new layer thickness (m)
497 :
498 : real (kind=dbl_kind), intent(in) :: &
499 : hn ! total thickness (m)
500 :
501 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
502 : qn ! layer quantity (enthalpy, salinity...)
503 :
504 : ! local variables
505 :
506 : integer (kind=int_kind) :: &
507 : k, k1, k2 ! vertical indices
508 :
509 : real (kind=dbl_kind) :: &
510 : hovlp ! overlap between old and new layers (m)
511 :
512 : real (kind=dbl_kind) :: &
513 : rhlyr, & ! 1./hlyr
514 : qtot ! total h*q in the column
515 :
516 : real (kind=dbl_kind), dimension (nlyr) :: &
517 41528350 : hq ! h * q for a layer
518 :
519 : character(len=*),parameter :: subname='(adjust_enthalpy)'
520 :
521 : !-----------------------------------------------------------------
522 : ! Compute reciprocal layer thickness.
523 : !-----------------------------------------------------------------
524 :
525 20764175 : rhlyr = c0
526 20764175 : if (hn > puny) then
527 20462188 : rhlyr = c1 / hlyr
528 :
529 : !-----------------------------------------------------------------
530 : ! Compute h*q for new layers (k2) given overlap with old layers (k1)
531 : !-----------------------------------------------------------------
532 :
533 153112310 : do k2 = 1, nlyr
534 153112310 : hq(k2) = c0
535 : enddo ! k
536 20462188 : k1 = 1
537 20462188 : k2 = 1
538 264056825 : do while (k1 <= nlyr .and. k2 <= nlyr)
539 : hovlp = min (z1(k1+1), z2(k2+1)) &
540 243594637 : - max (z1(k1), z2(k2))
541 243594637 : hovlp = max (hovlp, c0)
542 243594637 : hq(k2) = hq(k2) + hovlp*qn(k1)
543 243594637 : if (z1(k1+1) > z2(k2+1)) then
544 112747620 : k2 = k2 + 1
545 : else
546 130847017 : k1 = k1 + 1
547 : endif
548 : enddo ! while
549 :
550 : !-----------------------------------------------------------------
551 : ! Compute new enthalpies.
552 : !-----------------------------------------------------------------
553 :
554 153112310 : do k = 1, nlyr
555 153112310 : qn(k) = hq(k) * rhlyr
556 : enddo ! k
557 :
558 : else
559 :
560 301987 : qtot = c0
561 1811938 : do k = 1, nlyr
562 1811938 : qtot = qtot + qn(k) * (z1(k+1)-z1(k))
563 : enddo
564 301987 : if (hn > c0) then
565 1626 : do k = 1, nlyr
566 1626 : qn(k) = qtot/hn
567 : enddo
568 : else
569 1810312 : do k = 1, nlyr
570 1810312 : qn(k) = c0
571 : enddo
572 : endif
573 :
574 : endif
575 :
576 20764175 : end subroutine adjust_enthalpy
577 :
578 : !=======================================================================
579 :
580 : end module icepack_therm_shared
581 :
582 : !=======================================================================
|
