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
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, heat_capacity, 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: enthalpy_mush
21 : use icepack_mushy_physics, only: temperature_snow
22 : use icepack_mushy_physics, only: enthalpy_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_thermo, &
33 : icepack_init_trcr, &
34 : icepack_ice_temperature, &
35 : icepack_snow_temperature, &
36 : icepack_liquidus_temperature, &
37 : icepack_sea_freezing_temperature, &
38 : icepack_enthalpy_snow
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 : real (kind=dbl_kind), parameter, public :: &
51 : hfrazilmin = 0.05_dbl_kind ! min thickness of new frazil ice (m)
52 :
53 : real (kind=dbl_kind), public :: &
54 : hi_min ! minimum ice thickness allowed (m)
55 :
56 : !=======================================================================
57 :
58 : contains
59 :
60 : !=======================================================================
61 : !
62 : ! Compute the internal ice temperatures from enthalpy using
63 : ! quadratic formula
64 :
65 9054323 : function calculate_Tin_from_qin (qin, Tmltk) &
66 3232100 : result(Tin)
67 :
68 : real (kind=dbl_kind), intent(in) :: &
69 : qin , & ! enthalpy
70 : Tmltk ! melting temperature at one level
71 :
72 : real (kind=dbl_kind) :: &
73 : Tin ! internal temperature
74 :
75 : ! local variables
76 :
77 : real (kind=dbl_kind) :: &
78 3232100 : aa1,bb1,cc1 ! quadratic solvers
79 :
80 : character(len=*),parameter :: subname='(calculate_Tin_from_qin)'
81 :
82 9054323 : if (l_brine) then
83 8751050 : aa1 = cp_ice
84 8751050 : bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh
85 8751050 : cc1 = Lfresh * Tmltk
86 : Tin = min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) / &
87 8751050 : (c2*aa1),Tmltk)
88 :
89 : else ! fresh ice
90 303273 : Tin = (Lfresh + qin/rhoi) / cp_ice
91 : endif
92 :
93 9054323 : end function calculate_Tin_from_qin
94 :
95 : !=======================================================================
96 : ! Surface heat flux
97 : !=======================================================================
98 :
99 : ! heat flux into ice
100 :
101 11631532 : 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 3942918 : TsfK , & ! ice/snow surface temperature (K)
132 3942918 : Qsfc , & ! saturated surface specific humidity (kg/kg)
133 3942918 : qsat , & ! the saturation humidity of air (kg/m^3)
134 3942918 : flwdabs , & ! downward longwave absorbed heat flx (W/m^2)
135 3942918 : tmpvar ! 1/TsfK
136 :
137 : character(len=*),parameter :: subname='(surface_heat_flux)'
138 :
139 : ! ice surface temperature in Kelvin
140 11631532 : TsfK = Tsf + Tffresh
141 : ! TsfK = max(Tsf + Tffresh, c1)
142 11631532 : tmpvar = c1/TsfK
143 :
144 : ! saturation humidity
145 11631532 : qsat = qqqice * exp(-TTTice*tmpvar)
146 11631532 : Qsfc = qsat / rhoa
147 :
148 : ! longwave radiative flux
149 11631532 : flwdabs = emissivity * flw
150 11631532 : flwoutn = -emissivity * stefan_boltzmann * TsfK**4
151 :
152 : ! downward latent and sensible heat fluxes
153 11631532 : fsensn = shcoef * (potT - TsfK)
154 11631532 : flatn = lhcoef * (Qa - Qsfc)
155 :
156 : ! combine fluxes
157 11631532 : fsurfn = fswsfc + flwdabs + flwoutn + fsensn + flatn
158 :
159 11631532 : end subroutine surface_heat_flux
160 :
161 : !=======================================================================
162 :
163 8637565 : 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 2940290 : TsfK , & ! ice/snow surface temperature (K)
190 2940290 : dQsfc_dTsf , & ! saturated surface specific humidity (kg/kg)
191 2940290 : qsat , & ! the saturation humidity of air (kg/m^3)
192 2940290 : 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 8637565 : TsfK = Tsf + Tffresh
199 8637565 : tmpvar = c1/TsfK
200 :
201 : ! saturation humidity
202 8637565 : qsat = qqqice * exp(-TTTice*tmpvar)
203 8637565 : dQsfc_dTsf = TTTice * tmpvar * tmpvar * (qsat / rhoa)
204 :
205 : ! longwave radiative flux
206 8637565 : dflwoutn_dTsf = -emissivity * stefan_boltzmann * c4*TsfK**3
207 :
208 : ! downward latent and sensible heat fluxes
209 8637565 : dfsensn_dTsf = -shcoef
210 8637565 : dflatn_dTsf = -lhcoef * dQsfc_dTsf
211 :
212 : ! combine fluxes
213 8637565 : dfsurfn_dTsf = dflwoutn_dTsf + dfsensn_dTsf + dflatn_dTsf
214 :
215 8637565 : end subroutine dsurface_heat_flux_dTsf
216 :
217 : !=======================================================================
218 : !autodocument_start icepack_init_thermo
219 : ! Initialize the vertical profile of ice salinity and melting temperature.
220 : !
221 : ! authors: C. M. Bitz, UW
222 : ! William H. Lipscomb, LANL
223 :
224 42 : subroutine icepack_init_thermo(nilyr, sprofile)
225 :
226 : integer (kind=int_kind), intent(in) :: &
227 : nilyr ! number of ice layers
228 :
229 : real (kind=dbl_kind), dimension(:), intent(out) :: &
230 : sprofile ! vertical salinity profile
231 :
232 : !autodocument_end
233 :
234 : real (kind=dbl_kind), parameter :: &
235 : nsal = 0.407_dbl_kind, &
236 : msal = 0.573_dbl_kind
237 :
238 : integer (kind=int_kind) :: k ! ice layer index
239 15 : real (kind=dbl_kind) :: zn ! normalized ice thickness
240 :
241 : character(len=*),parameter :: subname='(icepack_init_thermo)'
242 :
243 : !-----------------------------------------------------------------
244 : ! Determine l_brine based on saltmax.
245 : ! Thermodynamic solver will not converge if l_brine is true and
246 : ! saltmax is close to zero.
247 : ! Set l_brine to false for zero layer thermodynamics
248 : !-----------------------------------------------------------------
249 :
250 42 : heat_capacity = .true.
251 42 : if (ktherm == 0) heat_capacity = .false. ! 0-layer thermodynamics
252 :
253 42 : l_brine = .false.
254 42 : if (saltmax > min_salin .and. heat_capacity) l_brine = .true.
255 :
256 : !-----------------------------------------------------------------
257 : ! Prescibe vertical profile of salinity and melting temperature.
258 : ! Note this profile is only used for BL99 thermodynamics.
259 : !-----------------------------------------------------------------
260 :
261 42 : if (l_brine) then
262 288 : do k = 1, nilyr
263 : zn = (real(k,kind=dbl_kind)-p5) / &
264 252 : real(nilyr,kind=dbl_kind)
265 252 : sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
266 288 : sprofile(k) = max(sprofile(k), min_salin)
267 : enddo ! k
268 36 : sprofile(nilyr+1) = saltmax
269 :
270 : else ! .not. l_brine
271 18 : do k = 1, nilyr+1
272 18 : sprofile(k) = c0
273 : enddo
274 : endif ! l_brine
275 :
276 42 : end subroutine icepack_init_thermo
277 :
278 : !=======================================================================
279 : !autodocument_start icepack_init_trcr
280 : !
281 396 : subroutine icepack_init_trcr(Tair, Tf, &
282 396 : Sprofile, Tprofile, &
283 : Tsfc, &
284 : nilyr, nslyr, &
285 792 : qin, qsn)
286 :
287 : integer (kind=int_kind), intent(in) :: &
288 : nilyr, & ! number of ice layers
289 : nslyr ! number of snow layers
290 :
291 : real (kind=dbl_kind), intent(in) :: &
292 : Tair, & ! air temperature (C)
293 : Tf ! freezing temperature (C)
294 :
295 : real (kind=dbl_kind), dimension(:), intent(in) :: &
296 : Sprofile, & ! vertical salinity profile (ppt)
297 : Tprofile ! vertical temperature profile (C)
298 :
299 : real (kind=dbl_kind), intent(out) :: &
300 : Tsfc ! surface temperature (C)
301 :
302 : real (kind=dbl_kind), dimension(:), intent(out) :: &
303 : qin, & ! ice enthalpy profile (J/m3)
304 : qsn ! snow enthalpy profile (J/m3)
305 :
306 : !autodocument_end
307 :
308 : ! local variables
309 :
310 : integer (kind=int_kind) :: k
311 :
312 : real (kind=dbl_kind) :: &
313 142 : slope, Ti
314 :
315 : character(len=*),parameter :: subname='(icepack_init_trcr)'
316 :
317 : ! surface temperature
318 396 : Tsfc = Tf ! default
319 396 : if (calc_Tsfc) Tsfc = min(Tsmelt, Tair - Tffresh) ! deg C
320 :
321 396 : if (heat_capacity) then
322 :
323 : ! ice enthalpy
324 2880 : do k = 1, nilyr
325 : ! assume linear temp profile and compute enthalpy
326 2520 : slope = Tf - Tsfc
327 : Ti = Tsfc + slope*(real(k,kind=dbl_kind)-p5) &
328 2520 : /real(nilyr,kind=dbl_kind)
329 2880 : if (ktherm == 2) then
330 1890 : qin(k) = enthalpy_mush(Ti, Sprofile(k))
331 : else
332 210 : qin(k) = -(rhoi * (cp_ice*(Tprofile(k)-Ti) &
333 630 : + Lfresh*(c1-Tprofile(k)/Ti) - cp_ocn*Tprofile(k)))
334 : endif
335 : enddo ! nilyr
336 :
337 : ! snow enthalpy
338 840 : do k = 1, nslyr
339 480 : Ti = min(c0, Tsfc)
340 840 : qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
341 : enddo ! nslyr
342 :
343 : else ! one layer with zero heat capacity
344 :
345 : ! ice energy
346 36 : qin(1) = -rhoi * Lfresh
347 :
348 : ! snow energy
349 36 : qsn(1) = -rhos * Lfresh
350 :
351 : endif ! heat_capacity
352 :
353 396 : end subroutine icepack_init_trcr
354 :
355 : !=======================================================================
356 : !autodocument_start icepack_liquidus_temperature
357 : ! compute liquidus temperature
358 :
359 791160 : function icepack_liquidus_temperature(Sin) result(Tmlt)
360 :
361 : real(dbl_kind), intent(in) :: Sin
362 : real(dbl_kind) :: Tmlt
363 :
364 : !autodocument_end
365 :
366 : character(len=*),parameter :: subname='(icepack_liquidus_temperature)'
367 :
368 791160 : if (ktherm == 2) then
369 :
370 791100 : Tmlt = liquidus_temperature_mush(Sin)
371 :
372 : else
373 :
374 60 : Tmlt = -depressT * Sin
375 :
376 : endif
377 :
378 791160 : end function icepack_liquidus_temperature
379 :
380 : !=======================================================================
381 : !autodocument_start icepack_sea_freezing_temperature
382 : ! compute ocean freezing temperature
383 :
384 1273632 : function icepack_sea_freezing_temperature(sss) result(Tf)
385 :
386 : real(dbl_kind), intent(in) :: sss
387 : real(dbl_kind) :: Tf
388 :
389 : !autodocument_end
390 :
391 : character(len=*),parameter :: subname='(icepack_sea_freezing_temperature)'
392 :
393 1273632 : if (trim(tfrz_option) == 'mushy') then
394 :
395 790992 : Tf = icepack_liquidus_temperature(sss) ! deg C
396 :
397 482640 : elseif (trim(tfrz_option) == 'linear_salt') then
398 :
399 386112 : Tf = -depressT * sss ! deg C
400 :
401 : else
402 :
403 96528 : Tf = -1.8_dbl_kind
404 :
405 : endif
406 :
407 1273632 : end function icepack_sea_freezing_temperature
408 :
409 : !=======================================================================
410 : !autodocument_start icepack_ice_temperature
411 : ! compute ice temperature
412 :
413 0 : function icepack_ice_temperature(qin, Sin) result(Tin)
414 :
415 : real(kind=dbl_kind), intent(in) :: qin, Sin
416 : real(kind=dbl_kind) :: Tin
417 :
418 : !autodocument_end
419 :
420 0 : real(kind=dbl_kind) :: Tmlts
421 :
422 : character(len=*),parameter :: subname='(icepack_ice_temperature)'
423 :
424 0 : if (ktherm == 2) then
425 :
426 0 : Tin = icepack_mushy_temperature_mush(qin, Sin)
427 :
428 : else
429 :
430 0 : Tmlts = -depressT * Sin
431 0 : Tin = calculate_Tin_from_qin(qin,Tmlts)
432 :
433 : endif
434 :
435 0 : end function icepack_ice_temperature
436 :
437 : !=======================================================================
438 : !autodocument_start icepack_snow_temperature
439 : ! compute snow temperature
440 :
441 0 : function icepack_snow_temperature(qin) result(Tsn)
442 :
443 : real(kind=dbl_kind), intent(in) :: qin
444 : real(kind=dbl_kind) :: Tsn
445 :
446 : !autodocument_end
447 :
448 : character(len=*),parameter :: subname='(icepack_snow_temperature)'
449 :
450 0 : if (ktherm == 2) then
451 :
452 0 : Tsn = temperature_snow(qin)
453 :
454 : else
455 :
456 0 : Tsn = (Lfresh + qin/rhos)/cp_ice
457 :
458 : endif
459 :
460 0 : end function icepack_snow_temperature
461 :
462 : !=======================================================================
463 : !autodocument_start icepack_enthalpy_snow
464 : ! compute snow enthalpy
465 :
466 0 : function icepack_enthalpy_snow(zTsn) result(qsn)
467 :
468 : real(kind=dbl_kind), intent(in) :: zTsn
469 : real(kind=dbl_kind) :: qsn
470 :
471 : !autodocument_end
472 :
473 : character(len=*),parameter :: subname='(icepack_enthalpy_snow)'
474 :
475 0 : qsn = enthalpy_snow(zTsn)
476 :
477 0 : end function icepack_enthalpy_snow
478 :
479 : !=======================================================================
480 :
481 : end module icepack_therm_shared
482 :
483 : !=======================================================================
|
