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 2621024115 : function calculate_Tin_from_qin (qin, Tmltk) &
66 46425153 : 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 46425153 : aa1,bb1,cc1 ! quadratic solvers
79 :
80 : character(len=*),parameter :: subname='(calculate_Tin_from_qin)'
81 :
82 2621024115 : if (l_brine) then
83 2383335122 : aa1 = cp_ice
84 2383335122 : bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh
85 2383335122 : cc1 = Lfresh * Tmltk
86 : Tin = min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) / &
87 2383335122 : (c2*aa1),Tmltk)
88 :
89 : else ! fresh ice
90 237688993 : Tin = (Lfresh + qin/rhoi) / cp_ice
91 : endif
92 :
93 2621024115 : end function calculate_Tin_from_qin
94 :
95 : !=======================================================================
96 : ! Surface heat flux
97 : !=======================================================================
98 :
99 : ! heat flux into ice
100 :
101 2421875136 : 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 154041222 : TsfK , & ! ice/snow surface temperature (K)
132 154041222 : Qsfc , & ! saturated surface specific humidity (kg/kg)
133 154041222 : qsat , & ! the saturation humidity of air (kg/m^3)
134 154041222 : flwdabs , & ! downward longwave absorbed heat flx (W/m^2)
135 154041222 : tmpvar ! 1/TsfK
136 :
137 : character(len=*),parameter :: subname='(surface_heat_flux)'
138 :
139 : ! ice surface temperature in Kelvin
140 2421875136 : TsfK = Tsf + Tffresh
141 : ! TsfK = max(Tsf + Tffresh, c1)
142 2421875136 : tmpvar = c1/TsfK
143 :
144 : ! saturation humidity
145 2421875136 : qsat = qqqice * exp(-TTTice*tmpvar)
146 2421875136 : Qsfc = qsat / rhoa
147 :
148 : ! longwave radiative flux
149 2421875136 : flwdabs = emissivity * flw
150 2421875136 : flwoutn = -emissivity * stefan_boltzmann * TsfK**4
151 :
152 : ! downward latent and sensible heat fluxes
153 2421875136 : fsensn = shcoef * (potT - TsfK)
154 2421875136 : flatn = lhcoef * (Qa - Qsfc)
155 :
156 : ! combine fluxes
157 2421875136 : fsurfn = fswsfc + flwdabs + flwoutn + fsensn + flatn
158 :
159 2421875136 : end subroutine surface_heat_flux
160 :
161 : !=======================================================================
162 :
163 1756327763 : 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 117234406 : TsfK , & ! ice/snow surface temperature (K)
190 117234406 : dQsfc_dTsf , & ! saturated surface specific humidity (kg/kg)
191 117234406 : qsat , & ! the saturation humidity of air (kg/m^3)
192 117234406 : 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 1756327763 : TsfK = Tsf + Tffresh
199 1756327763 : tmpvar = c1/TsfK
200 :
201 : ! saturation humidity
202 1756327763 : qsat = qqqice * exp(-TTTice*tmpvar)
203 1756327763 : dQsfc_dTsf = TTTice * tmpvar * tmpvar * (qsat / rhoa)
204 :
205 : ! longwave radiative flux
206 1756327763 : dflwoutn_dTsf = -emissivity * stefan_boltzmann * c4*TsfK**3
207 :
208 : ! downward latent and sensible heat fluxes
209 1756327763 : dfsensn_dTsf = -shcoef
210 1756327763 : dflatn_dTsf = -lhcoef * dQsfc_dTsf
211 :
212 : ! combine fluxes
213 1756327763 : dfsurfn_dTsf = dflwoutn_dTsf + dfsensn_dTsf + dflatn_dTsf
214 :
215 1756327763 : 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 2828 : 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 331 : 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 2828 : heat_capacity = .true.
251 2828 : if (ktherm == 0) heat_capacity = .false. ! 0-layer thermodynamics
252 :
253 2828 : l_brine = .false.
254 2828 : 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 2828 : if (l_brine) then
262 20826 : do k = 1, nilyr
263 : zn = (real(k,kind=dbl_kind)-p5) / &
264 18222 : real(nilyr,kind=dbl_kind)
265 18222 : sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
266 20826 : sprofile(k) = max(sprofile(k), min_salin)
267 : enddo ! k
268 2604 : sprofile(nilyr+1) = saltmax
269 :
270 : else ! .not. l_brine
271 672 : do k = 1, nilyr+1
272 672 : sprofile(k) = c0
273 : enddo
274 : endif ! l_brine
275 :
276 2828 : end subroutine icepack_init_thermo
277 :
278 : !=======================================================================
279 : !autodocument_start icepack_init_trcr
280 : !
281 1838640 : subroutine icepack_init_trcr(Tair, Tf, &
282 1838640 : Sprofile, Tprofile, &
283 : Tsfc, &
284 : nilyr, nslyr, &
285 3677280 : 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 372280 : slope, Ti
314 :
315 : character(len=*),parameter :: subname='(icepack_init_trcr)'
316 :
317 : ! surface temperature
318 1838640 : Tsfc = Tf ! default
319 1838640 : if (calc_Tsfc) Tsfc = min(Tsmelt, Tair - Tffresh) ! deg C
320 :
321 1838640 : if (heat_capacity) then
322 :
323 : ! ice enthalpy
324 7388960 : do k = 1, nilyr
325 : ! assume linear temp profile and compute enthalpy
326 6443680 : slope = Tf - Tsfc
327 : Ti = Tsfc + slope*(real(k,kind=dbl_kind)-p5) &
328 6443680 : /real(nilyr,kind=dbl_kind)
329 7388960 : if (ktherm == 2) then
330 4526200 : qin(k) = enthalpy_mush(Ti, Sprofile(k))
331 : else
332 147560 : qin(k) = -(rhoi * (cp_ice*(Tprofile(k)-Ti) &
333 1917480 : + Lfresh*(c1-Tprofile(k)/Ti) - cp_ocn*Tprofile(k)))
334 : endif
335 : enddo ! nilyr
336 :
337 : ! snow enthalpy
338 1890560 : do k = 1, nslyr
339 945280 : Ti = min(c0, Tsfc)
340 1890560 : qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
341 : enddo ! nslyr
342 :
343 : else ! one layer with zero heat capacity
344 :
345 : ! ice energy
346 893360 : qin(1) = -rhoi * Lfresh
347 :
348 : ! snow energy
349 893360 : qsn(1) = -rhos * Lfresh
350 :
351 : endif ! heat_capacity
352 :
353 1838640 : end subroutine icepack_init_trcr
354 :
355 : !=======================================================================
356 : !autodocument_start icepack_liquidus_temperature
357 : ! compute liquidus temperature
358 :
359 5720440 : 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 5720440 : if (ktherm == 2) then
369 :
370 4665264 : Tmlt = liquidus_temperature_mush(Sin)
371 :
372 : else
373 :
374 1055176 : Tmlt = -depressT * Sin
375 :
376 : endif
377 :
378 5720440 : end function icepack_liquidus_temperature
379 :
380 : !=======================================================================
381 : !autodocument_start icepack_sea_freezing_temperature
382 : ! compute ocean freezing temperature
383 :
384 0 : 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 0 : if (trim(tfrz_option) == 'mushy') then
394 :
395 0 : Tf = icepack_liquidus_temperature(sss) ! deg C
396 :
397 0 : elseif (trim(tfrz_option) == 'linear_salt') then
398 :
399 0 : Tf = -depressT * sss ! deg C
400 :
401 : else
402 :
403 0 : Tf = -1.8_dbl_kind
404 :
405 : endif
406 :
407 0 : end function icepack_sea_freezing_temperature
408 :
409 : !=======================================================================
410 : !autodocument_start icepack_ice_temperature
411 : ! compute ice temperature
412 :
413 8477392520 : 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 870809100 : real(kind=dbl_kind) :: Tmlts
421 :
422 : character(len=*),parameter :: subname='(icepack_ice_temperature)'
423 :
424 7606583420 : if (ktherm == 2) then
425 :
426 5564197520 : Tin = icepack_mushy_temperature_mush(qin, Sin)
427 :
428 : else
429 :
430 2042385900 : Tmlts = -depressT * Sin
431 2042385900 : Tin = calculate_Tin_from_qin(qin,Tmlts)
432 :
433 : endif
434 :
435 7606583420 : end function icepack_ice_temperature
436 :
437 : !=======================================================================
438 : !autodocument_start icepack_snow_temperature
439 : ! compute snow temperature
440 :
441 1193284460 : 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 1193284460 : if (ktherm == 2) then
451 :
452 794885360 : Tsn = temperature_snow(qin)
453 :
454 : else
455 :
456 398399100 : Tsn = (Lfresh + qin/rhos)/cp_ice
457 :
458 : endif
459 :
460 1193284460 : 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 : !=======================================================================
|
