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