CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
flake.F90
1
2!=======================================================================
3! Current Code Owner: DWD, Ulrich Schaettler
4! phone: +49 69 8062 2739
5! fax: +49 69 8236 1493
6! email: uschaettler@dwd.d400.de
7!
8! History:
9! Version Date Name
10! ---------- ---------- ----
11! 1.1 1998/03/11 Ulrich Schaettler
12! Initial release
13! 1.8 1998/08/03 Ulrich Schaettler
14! Eliminated intgribf, intgribc, irealgrib, iwlength and put it to data_io.
15! 1.10 1998/09/29 Ulrich Schaettler
16! Eliminated parameters for grid point and diagnostic calculations.
17! !VERSION! !DATE! <Your name>
18! <Modification comments>
19!
20! Code Description:
21! Language: Fortran 90.
22! Software Standards: "European Standards for Writing and
23! Documenting Exchangeable Fortran 90 Code".
24!
25! reorganize the FLake to module_FLake.F90 by Shaobo Zhang in 2016-7-13
26! added a new layer for deep lakes by Shaobo Zhang in 2016-11-15
27!
28!=======================================================================
29
30!------------------------------------------------------------------------------
31
33
34!------------------------------------------------------------------------------
35!
36! Description:
37! Global parameters for the program are defined.
38! Actually, scratch that. We'll import them from machine.F instead.
39!
40 use machine, only: ireals=>kind_phys, iintegers=>kind_integer
41
42!=======================================================================
43
44END MODULE data_parameters
45
46!------------------------------------------------------------------------------
47
49
50!------------------------------------------------------------------------------
51!
52! Description:
53!
54! This module contains "reference" values of albedo
55! for the lake water, lake ice and snow.
56! As in "flake_paramoptic_ref", two ice categories, viz. white ice and blue ice,
57! and two snow categories, viz. dry snow and melting snow, are used.
58!
59
60USE data_parameters, ONLY : &
61 ireals , & ! KIND-type parameter for real variables
62 iintegers ! KIND-type parameter for "normal" integer variables
63
64use machine, only: kind_phys
65!==============================================================================
66
67IMPLICIT NONE
68
69!==============================================================================
70!
71! Declarations
72
73! Albedo for water, ice and snow.
74!REAL (KIND = ireals), PARAMETER :: &
75! albedo_water_ref = 0.070 , & ! Water
76! albedo_whiteice_ref = 0.600 , & ! White ice
77! albedo_blueice_ref = 0.100 , & ! Blue ice
78! albedo_drysnow_ref = 0.600 , & ! Dry snow
79! albedo_meltingsnow_ref = 0.100 ! Melting snow
80
81! Empirical parameters.
82!REAL (KIND = ireals), PARAMETER :: &
83! c_albice_MR = 95.60 ! Constant in the interpolation formula for
84 ! the ice albedo (Mironov and Ritter 2004)
85! Albedo for water, ice and snow.
86real(kind = kind_phys), PARAMETER :: &
87 albedo_water_ref = 0.07 , & ! Water
88 albedo_whiteice_ref = 0.60 , & ! White ice
89 albedo_blueice_ref = 0.10 , & ! Blue ice
90! albedo_drysnow_ref = 0.60 , & ! Dry snow
91 albedo_drysnow_ref = 0.90 , & ! Dry snow
92 albedo_meltingsnow_ref = 0.10 ! Melting snow
93
94! Empirical parameters.
95real(kind = kind_phys), PARAMETER :: &
96 c_albice_mr = 95.6 ! Constant in the interpolation formula for
97 ! the ice albedo (Mironov and Ritter 2004)
98
99
100!==============================================================================
101
102END MODULE flake_albedo_ref
103
104!------------------------------------------------------------------------------
105
107
108!------------------------------------------------------------------------------
109!
110! Description:
111!
112! Switches and reference values of parameters
113! that configure the lake model FLake are set.
114!
115
116USE data_parameters , ONLY : &
117 ireals , & ! KIND-type parameter for real variables
118 iintegers ! KIND-type parameter for "normal" integer variables
119
120use machine, only: kind_phys
121!==============================================================================
122
123IMPLICIT NONE
124
125!==============================================================================
126!
127! Declarations
128! changed by Shaobo Zhang
129LOGICAL lflk_botsed_use
130!LOGICAL, PARAMETER :: &
131! lflk_botsed_use = .TRUE. ! .TRUE. indicates that the bottom-sediment scheme is used
132 ! to compute the depth penetrated by the thermal wave,
133 ! the temperature at this depth and the bottom heat flux.
134 ! Otherwise, the heat flux at the water-bottom sediment interface
135 ! is set to zero, the depth penetrated by the thermal wave
136 ! is set to a reference value defined below,
137 ! and the temperature at this depth is set to
138 ! the temperature of maximum density of the fresh water.
139
140!REAL (KIND = ireals), PARAMETER :: &
141! rflk_depth_bs_ref = 10.00 ! Reference value of the depth of the thermally active
142 ! layer of bottom sediments [m].
143 ! This value is used to (formally) define
144 ! the depth penetrated by the thermal wave
145 ! in case the bottom-sediment scheme is not used.
146
147real(kind = kind_phys), PARAMETER :: &
148 rflk_depth_bs_ref = 10.0
149
150!==============================================================================
151
152END MODULE flake_configure
153
154!------------------------------------------------------------------------------
155
157
158!------------------------------------------------------------------------------
159!
160! Description:
161!
162! Derived type(s) is(are) defined.
163!
164
165USE data_parameters , ONLY : &
166 ireals , & ! KIND-type parameter for real variables
167 iintegers ! KIND-type parameter for "normal" integer variables
168
169!==============================================================================
170
171IMPLICIT NONE
172
173!==============================================================================
174!
175! Declarations
176
177! Maximum value of the wave-length bands
178! in the exponential decay law for the radiation flux.
179! A storage for a ten-band approximation is allocated,
180! although a smaller number of bands is actually used.
181INTEGER (KIND = iintegers), PARAMETER :: &
182 nband_optic_max = 10_iintegers
183
184! Define TYPE "opticpar_medium"
186 INTEGER (KIND = iintegers) :: &
187 nband_optic ! Number of wave-length bands
188 REAL (kind = ireals), DIMENSION (nband_optic_max) :: &
189 frac_optic , & ! Fractions of total radiation flux
190 extincoef_optic ! Extinction coefficients
191END TYPE opticpar_medium
192
193!==============================================================================
194
195END MODULE flake_derivedtypes
196
197!------------------------------------------------------------------------------
198
200
201!------------------------------------------------------------------------------
202!
203! Description:
204!
205! This module contains "reference" values of the optical characteristics
206! of the lake water, lake ice and snow. These reference values may be used
207! if no information about the optical characteristics of the lake in question
208! is available. An exponential decay law for the solar radiation flux is assumed.
209! In the simplest one-band approximation,
210! the extinction coefficient for water is set to a large value,
211! leading to the absorption of 95% of the incoming radiation
212! within the uppermost 1 m of the lake water.
213! The extinction coefficients for ice and snow are taken from
214! Launiainen and Cheng (1998). The estimates for the ice correspond
215! to the uppermost 0.1 m of the ice layer and to the clear sky conditions
216! (see Table 2 in op. cit.).
217! Very large values of the extinction coefficients for ice and snow ("opaque")
218! can be used to prevent penetration of the solar radiation
219! through the snow-ice cover.
220!
221
222USE data_parameters, ONLY : &
223 ireals , & ! KIND-type parameter for real variables
224 iintegers ! KIND-type parameter for "normal" integer variables
225
226USE flake_derivedtypes, ONLY : &
227 nband_optic_max , & ! Maximum value of the wave-length bands
228 opticpar_medium ! Derived TYPE
229
230!==============================================================================
231
232IMPLICIT NONE
233
234!==============================================================================
235!
236! Declarations
237
238INTEGER (KIND = iintegers), PRIVATE :: & ! Help variable(s)
239 i ! DO loop index
240
241! Optical characteristics for water, ice and snow.
242! The simplest one-band approximation is used as a reference.
243!TYPE (opticpar_medium), PARAMETER :: &
244! opticpar_water_ref = opticpar_medium(1, & ! Water (reference)
245! (/1.0, (0.0,i=2,nband_optic_max)/), &
246! (/3.0, (1.E+100,i=2,nband_optic_max)/)) , &
247! opticpar_water_trans = opticpar_medium(2, & ! Transparent Water (two-band)
248! (/0.100, 0.900, (0.0,i=3,nband_optic_max)/), &
249! (/2.00, 0.200, (1.E+100,i=3,nband_optic_max)/)) , &
250!!_nu opticpar_water_trans = opticpar_medium(1, & ! Transparent Water (one-band)
251!!_nu (/1.0, (0.0,i=2,nband_optic_max)/), &
252!!_nu (/0.300, (1.E+100,i=2,nband_optic_max)/)) , &
253! opticpar_whiteice_ref = opticpar_medium(1, & ! White ice
254! (/1.0, (0.0,i=2,nband_optic_max)/), &
255! (/17.10, (1.E+100,i=2,nband_optic_max)/)) , &
256! opticpar_blueice_ref = opticpar_medium(1, & ! Blue ice
257! (/1.0, (0.0,i=2,nband_optic_max)/), &
258! (/8.40, (1.E+100,i=2,nband_optic_max)/)) , &
259! opticpar_drysnow_ref = opticpar_medium(1, & ! Dry snow
260! (/1.0, (0.0,i=2,nband_optic_max)/), &
261! (/25.00, (1.E+100,i=2,nband_optic_max)/)) , &
262! opticpar_meltingsnow_ref = opticpar_medium(1, & ! Melting snow
263! (/1.0, (0.0,i=2,nband_optic_max)/), &
264! (/15.00, (1.E+100,i=2,nband_optic_max)/)) , &
265! opticpar_ice_opaque = opticpar_medium(1, & ! Opaque ice
266! (/1.0, (0.0,i=2,nband_optic_max)/), &
267! (/1.0E+070, (1.E+100,i=2,nband_optic_max)/)) , &
268! opticpar_snow_opaque = opticpar_medium(1, & ! Opaque snow
269! (/1.0, (0.0,i=2,nband_optic_max)/), &
270! (/1.0E+070, (1.E+100,i=2,nband_optic_max)/))
271
272TYPE (opticpar_medium), PARAMETER :: &
273 opticpar_water_ref = opticpar_medium(1, & ! Water (reference)
274 (/1., (0.,i=2,nband_optic_max)/), &
275 (/3., (1.e+10,i=2,nband_optic_max)/)) , &
276 opticpar_water_trans = opticpar_medium(2, & ! Transparent Water (two-band)
277 (/0.10, 0.90, (0.,i=3,nband_optic_max)/), &
278 (/2.0, 0.20, (1.e+10,i=3,nband_optic_max)/)) , &
279 opticpar_whiteice_ref = opticpar_medium(1, & ! White ice
280 (/1., (0.,i=2,nband_optic_max)/), &
281 (/17.1, (1.e+10,i=2,nband_optic_max)/)) , &
282 opticpar_blueice_ref = opticpar_medium(1, & ! Blue ice
283 (/1., (0.,i=2,nband_optic_max)/), &
284 (/8.4, (1.e+10,i=2,nband_optic_max)/)) , &
285 opticpar_drysnow_ref = opticpar_medium(1, & ! Dry snow
286 (/1., (0.,i=2,nband_optic_max)/), &
287 (/25.0, (1.e+10,i=2,nband_optic_max)/)) , &
288 opticpar_meltingsnow_ref = opticpar_medium(1, & ! Melting snow
289 (/1., (0.,i=2,nband_optic_max)/), &
290 (/15.0, (1.e+10,i=2,nband_optic_max)/)) , &
291 opticpar_ice_opaque = opticpar_medium(1, & ! Opaque ice
292 (/1., (0.,i=2,nband_optic_max)/), &
293 (/1.0e+07, (1.e+10,i=2,nband_optic_max)/)) , &
294 opticpar_snow_opaque = opticpar_medium(1, & ! Opaque snow
295 (/1., (0.,i=2,nband_optic_max)/), &
296 (/1.0e+07, (1.e+10,i=2,nband_optic_max)/))
297
298
299!==============================================================================
300
301END MODULE flake_paramoptic_ref
302
303!------------------------------------------------------------------------------
304
306
307!------------------------------------------------------------------------------
308!
309! Description:
310!
311! Values of empirical constants of the lake model FLake
312! and of several thermodynamic parameters are set.
313!
314
315USE data_parameters , ONLY : &
316 ireals , & ! KIND-type parameter for real variables
317 iintegers ! KIND-type parameter for "normal" integer variables
318
319use machine, only: kind_phys
320
321!==============================================================================
322
323IMPLICIT NONE
324
325!==============================================================================
326!
327! Declarations
328
329! Dimensionless constants
330! in the equations for the mixed-layer depth
331! and for the shape factor with respect to the temperature profile in the thermocline
332real(kind = kind_phys), PARAMETER :: &
333! c_cbl_1 = 0.170 , & ! Constant in the CBL entrainment equation
334! c_cbl_2 = 1.0 , & ! Constant in the CBL entrainment equation
335! c_sbl_ZM_n = 0.50 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
336! c_sbl_ZM_s = 10.0 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
337! c_sbl_ZM_i = 20.0 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
338! c_relax_h = 0.0300 , & ! Constant in the relaxation equation for the SBL depth
339! c_relax_C = 0.00300 ! Constant in the relaxation equation for the shape factor
340 ! with respect to the temperature profile in the thermocline
341 c_cbl_1 = 0.17 , & ! Constant in the CBL entrainment equation
342 c_cbl_2 = 1. , & ! Constant in the CBL entrainment equation
343 c_sbl_zm_n = 0.5 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
344 c_sbl_zm_s = 10. , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
345 c_sbl_zm_i = 20. , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
346 c_relax_h = 0.030 , & ! Constant in the relaxation equation for the SBL depth
347 c_relax_c = 0.0030
348
349! Parameters of the shape functions
350! Indices refer to T - thermocline, S - snow, I - ice,
351! B1 - upper layer of the bottom sediments, B2 - lower layer of the bottom sediments.
352! "pr0" and "pr1" denote zeta derivatives of the corresponding shape function
353! at "zeta=0" ad "zeta=1", respectively.
354real(kind = kind_phys), PARAMETER :: &
355 c_t_min = 0.5 , & ! Minimum value of the shape factor C_T (thermocline)
356 c_t_max = 0.8 , & ! Maximum value of the shape factor C_T (thermocline)
357 phi_t_pr0_1 = 40.0/3.0 , & ! Constant in the expression for the T shape-function derivative
358 phi_t_pr0_2 = 20.0/3.0 , & ! Constant in the expression for the T shape-function derivative
359 c_tt_1 = 11.0/18.0 , & ! Constant in the expression for C_TT (thermocline)
360 c_tt_2 = 7.0/45.0 , & ! Constant in the expression for C_TT (thermocline)
361 c_b1 = 2.0/3.0 , & ! Shape factor (upper layer of bottom sediments)
362 c_b2 = 3.0/5.0 , & ! Shape factor (lower layer of bottom sediments)
363 phi_b1_pr0 = 2.0 , & ! B1 shape-function derivative
364 c_s_lin = 0.5 , & ! Shape factor (linear temperature profile in the snow layer)
365 phi_s_pr0_lin = 1.0 , & ! S shape-function derivative (linear profile)
366 c_i_lin = 0.5 , & ! Shape factor (linear temperature profile in the ice layer)
367 phi_i_pr0_lin = 1.0 , & ! I shape-function derivative (linear profile)
368 phi_i_pr1_lin = 1.0 , & ! I shape-function derivative (linear profile)
369 phi_i_ast_mr = 2.0 , & ! Constant in the MR2004 expression for I shape factor
370 c_i_mr = 1.0/12.0 , & ! Constant in the MR2004 expression for I shape factor
371 h_ice_max = 3.0 ! Maximum ice tickness in
372 ! the Mironov and Ritter (2004, MR2004) ice model [m]
373
374! Security constants
375real(kind = kind_phys), PARAMETER :: &
376 h_snow_min_flk = 1.0e-5 , & ! Minimum snow thickness [m]
377 h_ice_min_flk = 1.0e-9 , & ! Minimum ice thickness [m]
378 h_ml_min_flk = 1.0e-2 , & ! Minimum mixed-layer depth [m]
379 h_ml_max_flk = 1.0e+3 , & ! Maximum mixed-layer depth [m]
380 h_b1_min_flk = 1.0e-3 , & ! Minimum thickness of the upper layer of bottom sediments [m]
381 u_star_min_flk = 1.0e-6 ! Minimum value of the surface friction velocity [m s^{-1}]
382
383! Security constant(s)
384real(kind = kind_phys), PARAMETER :: &
385 c_small_flk = 1.0e-10 ! A small number
386
387! Thermodynamic parameters
388real(kind = kind_phys), PARAMETER :: &
389 tpl_grav = 9.81 , & ! Acceleration due to gravity [m s^{-2}]
390 tpl_t_r = 277.13 , & ! Temperature of maximum density of fresh water [K]
391 tpl_t_f = 273.15 , & ! Fresh water freezing point [K]
392 tpl_a_t = 1.6509e-05 , & ! Constant in the fresh-water equation of state [K^{-2}]
393 tpl_rho_w_r = 1.0e+03 , & ! Maximum density of fresh water [kg m^{-3}]
394 tpl_rho_i = 9.1e+02 , & ! Density of ice [kg m^{-3}]
395 tpl_rho_s_min = 1.0e+02 , & ! Minimum snow density [kg m^{-3}]
396 tpl_rho_s_max = 4.0e+02 , & ! Maximum snow density [kg m^{-3}]
397 tpl_gamma_rho_s = 2.0e+02 , & ! Empirical parameter [kg m^{-4}]
398 ! in the expression for the snow density
399 tpl_l_f = 3.3e+05 , & ! Latent heat of fusion [J kg^{-1}]
400 tpl_c_w = 4.2e+03 , & ! Specific heat of water [J kg^{-1} K^{-1}]
401 tpl_c_i = 2.1e+03 , & ! Specific heat of ice [J kg^{-1} K^{-1}]
402 tpl_c_s = 2.1e+03 , & ! Specific heat of snow [J kg^{-1} K^{-1}]
403 tpl_kappa_w = 5.46e-01 , & ! Molecular heat conductivity of water [J m^{-1} s^{-1} K^{-1}]
404 tpl_kappa_i = 2.29 , & ! Molecular heat conductivity of ice [J m^{-1} s^{-1} K^{-1}]
405 tpl_kappa_s_min = 0.2 , & ! Minimum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
406 tpl_kappa_s_max = 1.5 , & ! Maximum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
407 tpl_gamma_kappa_s = 1.3 ! Empirical parameter [J m^{-2} s^{-1} K^{-1}]
408 ! in the expression for the snow heat conductivity
409
410!==============================================================================
411
412END MODULE flake_parameters
413
414!------------------------------------------------------------------------------
415
416MODULE flake
417
418!------------------------------------------------------------------------------
419!
420! Description:
421!
422! The main program unit of the lake model FLake,
423! containing most of the FLake procedures.
424! Most FLake variables and local parameters are declared.
425!
426! FLake (Fresh-water Lake) is a lake model capable of predicting the surface temperature
427! in lakes of various depth on the time scales from a few hours to a year.
428! The model is based on a two-layer parametric representation of
429! the evolving temperature profile, where the structure of the stratified layer between the
430! upper mixed layer and the basin bottom, the lake thermocline,
431! is described using the concept of self-similarity of the temperature-depth curve.
432! The concept was put forward by Kitaigorodskii and Miropolsky (1970)
433! to describe the vertical temperature structure of the oceanic seasonal thermocline.
434! It has since been successfully used in geophysical applications.
435! The concept of self-similarity of the evolving temperature profile
436! is also used to describe the vertical structure of the thermally active upper layer
437! of bottom sediments and of the ice and snow cover.
438!
439! The lake model incorporates the heat budget equations
440! for the four layers in question, viz., snow, ice, water and bottom sediments,
441! developed with due regard for the vertically distributed character
442! of solar radiation heating.
443! The entrainment equation that incorporates the Zilitinkevich (1975) spin-up term
444! is used to compute the depth of a convectively-mixed layer.
445! A relaxation-type equation is used
446! to compute the wind-mixed layer depth in stable and neutral stratification,
447! where a multi-limit formulation for the equilibrium mixed-layer depth
448! proposed by Zilitinkevich and Mironov (1996)
449! accounts for the effects of the earth's rotation, of the surface buoyancy flux
450! and of the static stability in the thermocline.
451! The equations for the mixed-layer depth are developed with due regard for
452! the volumetric character of the radiation heating.
453! Simple thermodynamic arguments are invoked to develop
454! the evolution equations for the ice thickness and for the snow thickness.
455! The heat flux through the water-bottom sediment interface is computed,
456! using a parameterization proposed by Golosov et al. (1998).
457! The heat flux trough the air-water interface
458! (or through the air-ice or air-snow interface)
459! is provided by the driving atmospheric model.
460!
461! Empirical constants and parameters of the lake model
462! are estimated, using independent empirical and numerical data.
463! They should not be re-evaluated when the model is applied to a particular lake.
464! The only lake-specific parameters are the lake depth,
465! the optical characteristics of lake water,
466! the temperature at the bottom of the thermally active layer
467! of bottom sediments and the depth of that layer.
468!
469! A detailed description of the lake model is given in
470! Mironov, D. V., 2005:
471! Parameterization of Lakes in Numerical Weather Prediction.
472! Part 1: Description of a Lake Model.
473! Manuscript is available from the author.
474! Dmitrii Mironov
475! German Weather Service, Kaiserleistr. 29/35, D-63067 Offenbach am Main, Germany.
476! dmitrii.mironov@dwd.de
477!
478! Lines embraced with "!_tmp" contain temporary parts of the code.
479! Lines embraced/marked with "!_dev" may be replaced
480! as improved parameterizations are developed and tested.
481! Lines embraced/marked with "!_dm" are DM's comments
482! that may be helpful to a user.
483! Lines embraced/marked with "!_dbg" are used
484! for debugging purposes only.
485!
486
487USE data_parameters , ONLY : &
488 ireals , & ! KIND-type parameter for real variables
489 iintegers ! KIND-type parameter for "normal" integer variables
490
491use machine, only: kind_phys
492!==============================================================================
493
494IMPLICIT NONE
495
496!==============================================================================
497!
498! Declarations
499!
500! The variables declared below
501! are accessible to all program units of the MODULE flake.
502! Some of them should be USEd by the driving routines that call flake routines.
503! These are basically the quantities computed by FLake.
504! All variables declared below have a suffix "flk".
505
506! FLake variables of type REAL
507
508! Temperatures at the previous time step ("p") and the updated temperatures ("n")
509real(kind = kind_phys) :: &
510 t_mnw_p_flk, t_mnw_n_flk , & ! Mean temperature of the water column [K]
511 t_snow_p_flk, t_snow_n_flk , & ! Temperature at the air-snow interface [K]
512 t_ice_p_flk, t_ice_n_flk , & ! Temperature at the snow-ice or air-ice interface [K]
513 t_wml_p_flk, t_wml_n_flk , & ! Mixed-layer temperature [K]
514 t_bot_p_flk, t_bot_n_flk , & ! Temperature at the water-bottom sediment interface [K]
515 t_b1_p_flk, t_b1_n_flk ! Temperature at the bottom of the upper layer of the sediments [K]
516
517! Thickness of various layers at the previous time step ("p") and the updated values ("n")
518real(kind = kind_phys) :: &
519 h_snow_p_flk, h_snow_n_flk , & ! Snow thickness [m]
520 h_ice_p_flk, h_ice_n_flk , & ! Ice thickness [m]
521 h_ml_p_flk, h_ml_n_flk , & ! Thickness of the mixed-layer [m]
522 h_b1_p_flk, h_b1_n_flk ! Thickness of the upper layer of bottom sediments [m]
523
524! The shape factor(s) at the previous time step ("p") and the updated value(s) ("n")
525real(kind = kind_phys) :: &
526 c_t_p_flk, c_t_n_flk , & ! Shape factor (thermocline)
527 c_tt_flk , & ! Dimensionless parameter (thermocline)
528 c_q_flk , & ! Shape factor with respect to the heat flux (thermocline)
529 c_i_flk , & ! Shape factor (ice)
530 c_s_flk ! Shape factor (snow)
531
532! Derivatives of the shape functions
533real(kind = kind_phys) :: &
534 phi_t_pr0_flk , & ! d\Phi_T(0)/d\zeta (thermocline)
535 phi_i_pr0_flk , & ! d\Phi_I(0)/d\zeta_I (ice)
536 phi_i_pr1_flk , & ! d\Phi_I(1)/d\zeta_I (ice)
537 phi_s_pr0_flk ! d\Phi_S(0)/d\zeta_S (snow)
538
539! Heat and radiation fluxes
540real(kind = kind_phys) :: &
541 q_snow_flk , & ! Heat flux through the air-snow interface [W m^{-2}]
542 q_ice_flk , & ! Heat flux through the snow-ice or air-ice interface [W m^{-2}]
543 q_w_flk , & ! Heat flux through the ice-water or air-water interface [W m^{-2}]
544 q_bot_flk , & ! Heat flux through the water-bottom sediment interface [W m^{-2}]
545 i_atm_flk , & ! Radiation flux at the lower boundary of the atmosphere [W m^{-2}],
546 ! i.e. the incident radiation flux with no regard for the surface albedo.
547 i_snow_flk , & ! Radiation flux through the air-snow interface [W m^{-2}]
548 i_ice_flk , & ! Radiation flux through the snow-ice or air-ice interface [W m^{-2}]
549 i_w_flk , & ! Radiation flux through the ice-water or air-water interface [W m^{-2}]
550 i_h_flk , & ! Radiation flux through the mixed-layer-thermocline interface [W m^{-2}]
551 i_bot_flk , & ! Radiation flux through the water-bottom sediment interface [W m^{-2}]
552 i_intm_0_h_flk , & ! Mean radiation flux over the mixed layer [W m^{-1}]
553 i_intm_h_d_flk , & ! Mean radiation flux over the thermocline [W m^{-1}]
554 i_intm_d_h_flk , & ! Mean radiation flux over the deeper layer defined by Shaobo Zhang [W m^{-1}]
555 i_hh_flk , & ! Radiation flux through the bottom of the deeper layer defined by Shaobo Zhang [W m^{-2}]
556 q_star_flk ! A generalized heat flux scale [W m^{-2}]
557
558! Velocity scales
559real(kind = kind_phys) :: &
560 u_star_w_flk , & ! Friction velocity in the surface layer of lake water [m s^{-1}]
561 w_star_sfc_flk ! Convective velocity scale,
562 ! using a generalized heat flux scale [m s^{-1}]
563
564! The rate of snow accumulation
565real(kind = kind_phys) :: &
566 dmsnowdt_flk ! The rate of snow accumulation [kg m^{-2} s^{-1}]
567! The secondary layer temp
568real(kind = kind_phys) :: &
569 t_bot_2_in_flk
570
571!==============================================================================
572! Procedures
573!==============================================================================
574
575CONTAINS
576
577!==============================================================================
578! The codes of the FLake procedures are stored in separate "*.incf" files
579! and are included below.
580!------------------------------------------------------------------------------
581
582!==============================================================================
583! include 'flake_radflux.incf'
584!------------------------------------------------------------------------------
585! changed by Shaobo Zhang
586
587SUBROUTINE flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
588 opticpar_water, opticpar_ice, opticpar_snow, &
589 depth_bs )
590
591!------------------------------------------------------------------------------
592!
593! Description:
594!
595! Computes the radiation fluxes
596! at the snow-ice, ice-water, air-water,
597! mixed layer-thermocline and water column-bottom sediment interfaces,
598! the mean radiation flux over the mixed layer,
599! and the mean radiation flux over the thermocline.
600!
601!
602! Declarations:
603!
604! Modules used:
605
606!_dm Parameters are USEd in module "flake".
607!_nu USE data_parameters , ONLY : &
608!_nu ireals, & ! KIND-type parameter for real variables
609!_nu iintegers ! KIND-type parameter for "normal" integer variables
610
611USE flake_derivedtypes ! Definitions of derived TYPEs
612
613USE flake_parameters , ONLY : &
614 h_snow_min_flk , & ! Minimum snow thickness [m]
615 h_ice_min_flk , & ! Minimum ice thickness [m]
616 h_ml_min_flk ! Minimum mixed-layer depth [m]
617
618use machine, only: kind_phys
619!==============================================================================
620
621IMPLICIT NONE
622
623!==============================================================================
624!
625! Declarations
626
627! Input (procedure arguments)
628
629real(kind = kind_phys), INTENT(IN) :: &
630 depth_w , & ! The lake depth [m]
631 depth_bs , & ! The depth_bs added by Shaobo Zhang
632 albedo_water , & ! Albedo of the water surface
633 albedo_ice , & ! Albedo of the ice surface
634 albedo_snow ! Albedo of the snow surface
635
636TYPE (opticpar_medium), INTENT(IN) :: &
637 opticpar_water , & ! Optical characteristics of water
638 opticpar_ice , & ! Optical characteristics of ice
639 opticpar_snow ! Optical characteristics of snow
640
641
642! Local variables of type INTEGER
643INTEGER (KIND = iintegers) :: & ! Help variable(s)
644 i ! DO loop index
645
646!==============================================================================
647! Start calculations
648!------------------------------------------------------------------------------
649
650 IF(h_ice_p_flk.GE.h_ice_min_flk) THEN ! Ice exists
651 IF(h_snow_p_flk.GE.h_snow_min_flk) THEN ! There is snow above the ice
652 i_snow_flk = i_atm_flk*(1.0-albedo_snow)
653 i_bot_flk = 0.0
654 DO i=1, opticpar_snow%nband_optic
655 i_bot_flk = i_bot_flk + &
656 opticpar_snow%frac_optic(i)*exp(-opticpar_snow%extincoef_optic(i)*h_snow_p_flk)
657 END DO
658 i_ice_flk = i_snow_flk*i_bot_flk
659 ELSE ! No snow above the ice
660 i_snow_flk = i_atm_flk
661 i_ice_flk = i_atm_flk*(1.0-albedo_ice)
662 END IF
663 i_bot_flk = 0.0
664 DO i=1, opticpar_ice%nband_optic
665 i_bot_flk = i_bot_flk + &
666 opticpar_ice%frac_optic(i)*exp(-opticpar_ice%extincoef_optic(i)*h_ice_p_flk)
667 END DO
668 i_w_flk = i_ice_flk*i_bot_flk
669 ELSE ! No ice-snow cover
670 i_snow_flk = i_atm_flk
671 i_ice_flk = i_atm_flk
672 i_w_flk = i_atm_flk*(1.0-albedo_water)
673 END IF
674
675 IF(h_ml_p_flk.GE.h_ml_min_flk) THEN ! Radiation flux at the bottom of the mixed layer
676 i_bot_flk = 0.0
677 DO i=1, opticpar_water%nband_optic
678 i_bot_flk = i_bot_flk + &
679 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk)
680! print*,'nband_optic=',opticpar_water%nband_optic
681! print*,'Extinction=',opticpar_water%extincoef_optic(i)
682 END DO
683 i_h_flk = i_w_flk*i_bot_flk
684 ELSE ! Mixed-layer depth is less then a minimum value
685 i_h_flk = i_w_flk
686 END IF
687
688 i_bot_flk = 0.0 ! Radiation flux at the lake bottom
689 DO i=1, opticpar_water%nband_optic
690 i_bot_flk = i_bot_flk + &
691 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*depth_w)
692 END DO
693 i_bot_flk = i_w_flk*i_bot_flk
694
695 IF(h_ml_p_flk.GE.h_ml_min_flk) THEN ! Integral-mean radiation flux over the mixed layer
696 i_intm_0_h_flk = 0.0
697 DO i=1, opticpar_water%nband_optic
698 i_intm_0_h_flk = i_intm_0_h_flk + &
699 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
700 (1.0 - exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk))
701 END DO
702 i_intm_0_h_flk = i_w_flk*i_intm_0_h_flk/h_ml_p_flk
703 ELSE
704 i_intm_0_h_flk = i_h_flk
705 END IF
706
707 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk) THEN ! Integral-mean radiation flux over the thermocline
708 i_intm_h_d_flk = 0.0
709 DO i=1, opticpar_water%nband_optic
710 i_intm_h_d_flk = i_intm_h_d_flk + &
711 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
712 ( exp(-opticpar_water%extincoef_optic(i)*h_ml_p_flk) &
713 - exp(-opticpar_water%extincoef_optic(i)*depth_w) )
714 END DO
715 i_intm_h_d_flk = i_w_flk*i_intm_h_d_flk/(depth_w-h_ml_p_flk)
716 ELSE
717 i_intm_h_d_flk = i_h_flk
718 END IF
719
720! Added by Shaobo Zhang
721
722 IF(depth_bs.GE.h_ml_min_flk) then! Integral-mean radiation flux over the deeper layer defined by Shaobo Zhang
723 i_intm_d_h_flk = 0.0
724 DO i=1, opticpar_water%nband_optic
725 i_intm_d_h_flk = i_intm_d_h_flk + &
726 opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
727 ( exp(-opticpar_water%extincoef_optic(i)*depth_w) &
728 - exp(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs)) )
729 END DO
730 i_intm_d_h_flk = i_w_flk*i_intm_d_h_flk/depth_bs
731 ELSE
732 i_intm_d_h_flk = i_bot_flk
733 END IF
734
735! Radiation flux at the bottom of the deeper layer defined by Shaobo Zhang
736 i_hh_flk = 0.0
737 DO i=1, opticpar_water%nband_optic
738 i_hh_flk = i_hh_flk + &
739 opticpar_water%frac_optic(i)*exp(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs))
740 END DO
741 i_hh_flk = i_w_flk*i_hh_flk
742
743!------------------------------------------------------------------------------
744! End calculations
745!==============================================================================
746
747END SUBROUTINE flake_radflux
748
749!==============================================================================
750
751!==============================================================================
752! include 'flake_main.incf'
753!------------------------------------------------------------------------------
754
755SUBROUTINE flake_main ( depthw, depthbs, T_bs, par_Coriolis, &
756 extincoef_water_typ, &
757 del_time, T_sfc_p, T_sfc_n, T_bot_2_in, &
758 T_bot_2_out )
759
760!------------------------------------------------------------------------------
761!
762! Description:
763!
764! The main driving routine of the lake model FLake
765! where computations are performed.
766! Advances the surface temperature
767! and other FLake variables one time step.
768! At the moment, the Euler explicit scheme is used.
769!
770! Lines embraced with "!_tmp" contain temporary parts of the code.
771! Lines embraced/marked with "!_dev" may be replaced
772! as improved parameterizations are developed and tested.
773! Lines embraced/marked with "!_dm" are DM's comments
774! that may be helpful to a user.
775! Lines embraced/marked with "!_dbg" are used
776! for debugging purposes only.
777!
778! Declarations:
779!
780! Modules used:
781
782!_dm Parameters are USEd in module "flake".
783!_nu USE data_parameters , ONLY : &
784!_nu ireals, & ! KIND-type parameter for real variables
785!_nu iintegers ! KIND-type parameter for "normal" integer variables
786
787USE flake_parameters ! Thermodynamic parameters and dimensionless constants of FLake
788
789USE flake_configure ! Switches and parameters that configure FLake
790
791use machine, only: kind_phys
792! ADDED by Shaobo Zhang
793! USE mod_dynparam, only : lake_depth_max
794
795!==============================================================================
796
797IMPLICIT NONE
798
799!==============================================================================
800!
801! Declarations
802
803! Input (procedure arguments)
804
805! changed by Shaobo Zhang
806real(kind = kind_phys), INTENT(IN) :: &
807 depthw , & ! The lake depth [m]
808 depthbs , & ! Depth of the thermally active layer of bottom sediments [m]
809 t_bs , & ! Temperature at the outer edge of
810 ! the thermally active layer of bottom sediments [K]
811 par_coriolis , & ! The Coriolis parameter [s^{-1}]
812 extincoef_water_typ , & ! "Typical" extinction coefficient of the lake water [m^{-1}],
813 ! used to compute the equilibrium CBL depth
814 del_time , & ! The model time step [s]
815 t_sfc_p , & ! Surface temperature at the previous time step [K]
816 t_bot_2_in
817
818real(kind = kind_phys) :: &
819 depth_w , & ! The lake depth [m]
820 depth_bs ! Depth of the thermally active layer of bottom sediments [m]
821
822! Output (procedure arguments)
823
824real(kind = kind_phys), INTENT(OUT) :: &
825 t_sfc_n , & ! Updated surface temperature [K]
826 ! (equal to the updated value of either T_ice, T_snow or T_wML)
827 t_bot_2_out
828
829
830! Local variables of type LOGICAL
831LOGICAL :: &
832 l_ice_create , & ! Switch, .TRUE. = ice does not exist but should be created
833 l_snow_exists , & ! Switch, .TRUE. = there is snow above the ice
834 l_ice_meltabove ! Switch, .TRUE. = snow/ice melting from above takes place
835
836! Local variables of type INTEGER
837INTEGER (KIND = iintegers) :: &
838 i ! Loop index
839
840! Local variables of type REAL
841real(kind = kind_phys) :: &
842 d_t_mnw_dt , & ! Time derivative of T_mnw [K s^{-1}]
843 d_t_ice_dt , & ! Time derivative of T_ice [K s^{-1}]
844 d_t_bot_dt , & ! Time derivative of T_bot [K s^{-1}]
845 d_t_b1_dt , & ! Time derivative of T_B1 [K s^{-1}]
846 d_h_snow_dt , & ! Time derivative of h_snow [m s^{-1}]
847 d_h_ice_dt , & ! Time derivative of h_ice [m s^{-1}]
848 d_h_ml_dt , & ! Time derivative of h_ML [m s^{-1}]
849 d_h_b1_dt , & ! Time derivative of H_B1 [m s^{-1}]
850 d_h_d_dt , & ! Time derivative of h_D, new defined by Shaobo Zhang
851 d_t_h_dt , & ! Time derivative of T_H, new defined by Shaobo Zhang
852 d_c_t_dt ! Time derivative of C_T [s^{-1}]
853
854! Local variables of type REAL
855real(kind = kind_phys) :: &
856 n_t_mean , & ! The mean buoyancy frequency in the thermocline [s^{-1}]
857 tmp , & ! temperary variable
858 zm_h_scale , & ! The ZM96 equilibrium SBL depth scale [m]
859 conv_equil_h_scale ! The equilibrium CBL depth scale [m]
860
861! Local variables of type REAL
862real(kind = kind_phys) :: &
863 h_ice_threshold , & ! If h_ice<h_ice_threshold, use quasi-equilibrium ice model
864 flk_str_1 , & ! Help storage variable
865 flk_str_2 , & ! Help storage variable
866 r_h_icesnow , & ! Dimensionless ratio, used to store intermediate results
867 r_rho_c_icesnow , & ! Dimensionless ratio, used to store intermediate results
868 r_ti_icesnow , & ! Dimensionless ratio, used to store intermediate results
869 r_tstar_icesnow ! Dimensionless ratio, used to store intermediate results
870
871! ADDED by Shaobo Zhang
872!REAL (KIND = kind_phys) :: T_bot_2_in, T_bot_2_out
873real(kind = kind_phys) :: ct, ctt, cq
874real(kind = kind_phys) :: lake_depth_max
875
876!==============================================================================
877! Start calculations
878!------------------------------------------------------------------------------
879
880!_dm
881! Security. Set time-rate-of-change of prognostic variables to zero.
882! Set prognostic variables to their values at the previous time step.
883! (This is to avoid spurious changes of prognostic variables
884! when FLake is used within a 3D model, e.g. to avoid spurious generation of ice
885! at the neighbouring lake points as noticed by Burkhardt Rockel.)
886!_dm
887
888! print*,'T_sfc_p=',T_sfc_p,' T_bot_2_in=',T_bot_2_in
889! print*,'h_snow_p_flk=',h_snow_p_flk,' h_ice_p_flk=',h_ice_p_flk
890! print*,'T_snow_p_flk=',T_snow_p_flk, ' T_ice_p_flk=',T_ice_p_flk
891! print*,'T_wML_p_flk=',T_wML_p_flk, ' T_mnw_p_flk=',T_mnw_p_flk
892! print*,'T_bot_p_flk=',T_bot_p_flk, ' T_B1_p_flk=',T_B1_p_flk
893! print*,'h_ML_p_flk=',h_ML_p_flk, ' H_B1_p_flk=',H_B1_p_flk
894! print*,'C_T_p_flk=',C_T_p_flk
895! print*,'depthw= ', depthw
896! print*,'depthbs=',depthbs
897
898lake_depth_max = 60.0
899depth_w =depthw
900depth_bs = depthbs
901d_t_mnw_dt = 0.0
902d_t_ice_dt = 0.0
903d_t_bot_dt = 0.0
904d_t_b1_dt = 0.0
905d_h_snow_dt = 0.0
906d_h_ice_dt = 0.0
907d_h_ml_dt = 0.0
908d_h_b1_dt = 0.0
909d_c_t_dt = 0.0
910t_snow_n_flk = t_snow_p_flk
911t_ice_n_flk = t_ice_p_flk
912t_wml_n_flk = t_wml_p_flk
913t_mnw_n_flk = t_mnw_p_flk
914t_bot_n_flk = t_bot_p_flk
915t_b1_n_flk = t_b1_p_flk
916h_snow_n_flk = h_snow_p_flk
917h_ice_n_flk = h_ice_p_flk
918h_ml_n_flk = h_ml_p_flk
919h_b1_n_flk = h_b1_p_flk
920c_t_n_flk = c_t_p_flk
921
922!------------------------------------------------------------------------------
923! Compute fluxes, using variables from the previous time step.
924!------------------------------------------------------------------------------
925
926!_dm
927! At this point, the heat and radiation fluxes, namely,
928! Q_snow_flk, Q_ice_flk, Q_w_flk,
929! I_atm_flk, I_snow_flk, I_ice_flk, I_w_flk, I_h_flk, I_bot_flk,
930! the mean radiation flux over the mixed layer, I_intm_0_h_flk,
931! and the mean radiation flux over the thermocline, I_intm_h_D_flk,
932! should be known.
933! They are computed within "flake_interface" (or within the driving model)
934! and are available to "flake_main"
935! through the above variables declared in the MODULE "flake".
936! In case a lake is ice-covered, Q_w_flk is re-computed below.
937!_dm
938
939! Heat flux through the ice-water interface
940IF(h_ice_p_flk.GE.h_ice_min_flk) THEN ! Ice exists
941 IF(h_ml_p_flk.LE.h_ml_min_flk) THEN ! Mixed-layer depth is zero, compute flux
942 q_w_flk = -tpl_kappa_w*(t_bot_p_flk-t_wml_p_flk)/depth_w ! Flux with linear T(z)
943 phi_t_pr0_flk = phi_t_pr0_1*c_t_p_flk-phi_t_pr0_2 ! d\Phi(0)/d\zeta (thermocline)
944 q_w_flk = q_w_flk*max(phi_t_pr0_flk, 1.0) ! Account for an increased d\Phi(0)/d\zeta
945 ELSE
946 q_w_flk = 0.0 ! Mixed-layer depth is greater than zero, set flux to zero
947 END IF
948ELSE ! If Ice doesn't exist, set flux to zero, YWu 2019
949 q_w_flk = 0.0
950END IF
951
952! A generalized heat flux scale
953q_star_flk = q_w_flk + i_w_flk + i_h_flk - 2.0*i_intm_0_h_flk
954
955! changed by Shaobo Zhang
956! Heat flux through the water-bottom sediment interface
957!IF(lflk_botsed_use) THEN
958! Q_bot_flk = -tpl_kappa_w*(T_B1_p_flk-T_bot_p_flk)/MAX(H_B1_p_flk, H_B1_min_flk)*Phi_B1_pr0
959!ELSE
960! Q_bot_flk = 0.0 ! The bottom-sediment scheme is not used
961!END IF
962
963IF(lflk_botsed_use) THEN ! The bottom-sediment scheme is used
964 q_bot_flk = -tpl_kappa_w*(t_b1_p_flk-t_bot_p_flk)/max(h_b1_p_flk, h_b1_min_flk)*phi_b1_pr0
965ELSE ! The scheme written by Shaobo Zhang for the deeper layer of a deep lake is used
966 q_bot_flk = -tpl_kappa_w*(t_bot_2_in-t_bot_p_flk)/max(depth_bs, h_b1_min_flk)*phi_b1_pr0
967END IF
968
969!------------------------------------------------------------------------------
970! Check if ice exists or should be created.
971! If so, compute the thickness and the temperature of ice and snow.
972!------------------------------------------------------------------------------
973
974!_dm
975! Notice that a quasi-equilibrium ice-snow model is used
976! to avoid numerical instability when the ice is thin.
977! This is always the case when new ice is created.
978!_dm
979
980!_dev
981! The dependence of snow density and of snow heat conductivity
982! on the snow thickness is accounted for parametrically.
983! That is, the time derivatives of \rho_S and \kappa_S are neglected.
984! The exception is the equation for the snow thickness
985! in case of snow accumulation and no melting,
986! where d\rho_S/dt is incorporated.
987! Furthermore, some (presumably small) correction terms incorporating
988! the snow density and the snow heat conductivity are dropped out.
989! Those terms may be included as better formulations
990! for \rho_S and \kappa_S are available.
991!_dev
992
993! Default values
994l_ice_create = .false.
995l_ice_meltabove = .false.
996
997ice_exist: IF(h_ice_p_flk.LT.h_ice_min_flk) THEN ! Ice does not exist
998
999 l_ice_create = t_wml_p_flk.LE.(tpl_t_f+c_small_flk).AND.q_w_flk.LT.0.0
1000 IF(l_ice_create) THEN ! Ice does not exist but should be created
1001 d_h_ice_dt = -q_w_flk/tpl_rho_i/tpl_l_f
1002 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time ! Advance h_ice
1003 t_ice_n_flk = tpl_t_f + h_ice_n_flk*q_w_flk/tpl_kappa_i/phi_i_pr0_lin ! Ice temperature
1004 d_h_snow_dt = dmsnowdt_flk/tpl_rho_s_min
1005 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
1006 phi_i_pr1_flk = phi_i_pr1_lin &
1007 + phi_i_ast_mr*min(1.0, h_ice_n_flk/h_ice_max) ! d\Phi_I(1)/d\zeta_I (ice)
1008 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_n_flk) &
1009 * h_snow_n_flk/max(h_ice_n_flk, h_ice_min_flk)
1010 t_snow_n_flk = t_ice_n_flk + r_h_icesnow*(t_ice_n_flk-tpl_t_f) ! Snow temperature
1011 END IF
1012
1013ELSE ice_exist ! Ice exists
1014
1015 l_snow_exists = h_snow_p_flk.GE.h_snow_min_flk ! Check if there is snow above the ice
1016
1017 melting: IF(t_snow_p_flk.GE.(tpl_t_f-c_small_flk)) THEN ! T_sfc = T_f, check for melting from above
1018 ! T_snow = T_ice if snow is absent
1019 IF(l_snow_exists) THEN ! There is snow above the ice
1020 flk_str_1 = q_snow_flk + i_snow_flk - i_ice_flk ! Atmospheric forcing
1021 IF(flk_str_1.GE.0.0) THEN ! Melting of snow and ice from above
1022 l_ice_meltabove = .true.
1023 d_h_snow_dt = (-flk_str_1/tpl_l_f+dmsnowdt_flk)/flake_snowdensity(h_snow_p_flk)
1024 d_h_ice_dt = -(i_ice_flk - i_w_flk - q_w_flk)/tpl_l_f/tpl_rho_i
1025 END IF
1026 ELSE ! No snow above the ice
1027 flk_str_1 = q_ice_flk + i_ice_flk - i_w_flk - q_w_flk ! Atmospheric forcing + heating from the water
1028 IF(flk_str_1.GE.0.0) THEN ! Melting of ice from above, snow accumulation may occur
1029 l_ice_meltabove = .true.
1030 d_h_ice_dt = -flk_str_1/tpl_l_f/tpl_rho_i
1031 d_h_snow_dt = dmsnowdt_flk/tpl_rho_s_min
1032 END IF
1033 END IF
1034 IF(l_ice_meltabove) THEN ! Melting from above takes place
1035 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time ! Advance h_ice
1036 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
1037 t_ice_n_flk = tpl_t_f ! Set T_ice to the freezing point
1038 t_snow_n_flk = tpl_t_f ! Set T_snow to the freezing point
1039 END IF
1040
1041 END IF melting
1042
1043 no_melting: IF(.NOT.l_ice_meltabove) THEN ! No melting from above
1044
1045 d_h_snow_dt = flake_snowdensity(h_snow_p_flk)
1046 IF(d_h_snow_dt.LT.tpl_rho_s_max) THEN ! Account for d\rho_S/dt
1047 flk_str_1 = h_snow_p_flk*tpl_gamma_rho_s/tpl_rho_w_r
1048 flk_str_1 = flk_str_1/(1.0-flk_str_1)
1049 ELSE ! Snow density is equal to its maximum value, d\rho_S/dt=0
1050 flk_str_1 = 0.0
1051 END IF
1052 d_h_snow_dt = dmsnowdt_flk/d_h_snow_dt/(1.0+flk_str_1) ! Snow accumulation
1053 h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
1054
1055 phi_i_pr0_flk = h_ice_p_flk/h_ice_max ! h_ice relative to its maximum value
1056 c_i_flk = c_i_lin - c_i_mr*(1.0+phi_i_ast_mr)*phi_i_pr0_flk ! Shape factor (ice)
1057 phi_i_pr1_flk = phi_i_pr1_lin + phi_i_ast_mr*phi_i_pr0_flk ! d\Phi_I(1)/d\zeta_I (ice)
1058 phi_i_pr0_flk = phi_i_pr0_lin - phi_i_pr0_flk ! d\Phi_I(0)/d\zeta_I (ice)
1059
1060 h_ice_threshold = max(1.0, 2.0*c_i_flk*tpl_c_i*(tpl_t_f-t_ice_p_flk)/tpl_l_f)
1061 h_ice_threshold = phi_i_pr0_flk/c_i_flk*tpl_kappa_i/tpl_rho_i/tpl_c_i*h_ice_threshold
1062 h_ice_threshold = sqrt(h_ice_threshold*del_time) ! Threshold value of h_ice
1063 h_ice_threshold = min(0.9*h_ice_max, max(h_ice_threshold, h_ice_min_flk))
1064 ! h_ice(threshold) < 0.9*H_Ice_max
1065
1066 IF(h_ice_p_flk.LT.h_ice_threshold) THEN ! Use a quasi-equilibrium ice model
1067
1068 IF(l_snow_exists) THEN ! Use fluxes at the air-snow interface
1069 flk_str_1 = q_snow_flk + i_snow_flk - i_w_flk
1070 ELSE ! Use fluxes at the air-ice interface
1071 flk_str_1 = q_ice_flk + i_ice_flk - i_w_flk
1072 END IF
1073 d_h_ice_dt = -(flk_str_1-q_w_flk)/tpl_l_f/tpl_rho_i
1074 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time ! Advance h_ice
1075 t_ice_n_flk = tpl_t_f + h_ice_n_flk*flk_str_1/tpl_kappa_i/phi_i_pr0_flk ! Ice temperature
1076
1077 ELSE ! Use a complete ice model
1078
1079 d_h_ice_dt = tpl_kappa_i*(tpl_t_f-t_ice_p_flk)/h_ice_p_flk*phi_i_pr0_flk
1080 d_h_ice_dt = (q_w_flk+d_h_ice_dt)/tpl_l_f/tpl_rho_i
1081 h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time ! Advance h_ice
1082
1083 r_ti_icesnow = tpl_c_i*(tpl_t_f-t_ice_p_flk)/tpl_l_f ! Dimensionless parameter
1084 r_tstar_icesnow = 1.0 - c_i_flk ! Dimensionless parameter
1085 IF(l_snow_exists) THEN ! There is snow above the ice
1086 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_p_flk) &
1087 * h_snow_p_flk/h_ice_p_flk
1088 r_rho_c_icesnow = flake_snowdensity(h_snow_p_flk)*tpl_c_s/tpl_rho_i/tpl_c_i
1089!_dev
1090!_dm
1091! These terms should be included as an improved understanding of the snow scheme is gained,
1092! of the effect of snow density in particular.
1093!_dm
1094!_nu R_Tstar_icesnow = R_Tstar_icesnow &
1095!_nu + (1.0+C_S_lin*h_snow_p_flk/h_ice_p_flk)*R_H_icesnow*R_rho_c_icesnow
1096!_dev
1097
1098 r_tstar_icesnow = r_tstar_icesnow*r_ti_icesnow ! Dimensionless parameter
1099
1100!_dev
1101!_nu R_Tstar_icesnow = R_Tstar_icesnow &
1102!_nu + (1.0-R_rho_c_icesnow)*tpl_c_I*T_ice_p_flk/tpl_L_f
1103!_dev
1104 flk_str_2 = q_snow_flk+i_snow_flk-i_w_flk ! Atmospheric fluxes
1105 flk_str_1 = c_i_flk*h_ice_p_flk + (1.0+c_s_lin*r_h_icesnow)*r_rho_c_icesnow*h_snow_p_flk
1106 d_t_ice_dt = -(1.0-2.0*c_s_lin)*r_h_icesnow*(tpl_t_f-t_ice_p_flk) &
1107 * tpl_c_s*dmsnowdt_flk ! Effect of snow accumulation
1108 ELSE ! No snow above the ice
1109 r_tstar_icesnow = r_tstar_icesnow*r_ti_icesnow ! Dimensionless parameter
1110 flk_str_2 = q_ice_flk+i_ice_flk-i_w_flk ! Atmospheric fluxes
1111 flk_str_1 = c_i_flk*h_ice_p_flk
1112 d_t_ice_dt = 0.0
1113 END IF
1114 d_t_ice_dt = d_t_ice_dt + tpl_kappa_i*(tpl_t_f-t_ice_p_flk)/h_ice_p_flk*phi_i_pr0_flk &
1115 * (1.0-r_tstar_icesnow) ! Add flux due to heat conduction
1116 d_t_ice_dt = d_t_ice_dt - r_tstar_icesnow*q_w_flk ! Add flux from water to ice
1117 d_t_ice_dt = d_t_ice_dt + flk_str_2 ! Add atmospheric fluxes
1118 d_t_ice_dt = d_t_ice_dt/tpl_rho_i/tpl_c_i ! Total forcing
1119 d_t_ice_dt = d_t_ice_dt/flk_str_1 ! dT_ice/dt
1120 t_ice_n_flk = t_ice_p_flk + d_t_ice_dt*del_time ! Advance T_ice
1121 END IF
1122
1123 phi_i_pr1_flk = min(1.0, h_ice_n_flk/h_ice_max) ! h_ice relative to its maximum value
1124 phi_i_pr1_flk = phi_i_pr1_lin + phi_i_ast_mr*phi_i_pr1_flk ! d\Phi_I(1)/d\zeta_I (ice)
1125 r_h_icesnow = phi_i_pr1_flk/phi_s_pr0_lin*tpl_kappa_i/flake_snowheatconduct(h_snow_n_flk) &
1126 *h_snow_n_flk/max(h_ice_n_flk, h_ice_min_flk)
1127 t_snow_n_flk = t_ice_n_flk + r_h_icesnow*(t_ice_n_flk-tpl_t_f) ! Snow temperature
1128
1129 END IF no_melting
1130
1131END IF ice_exist
1132
1133! Security, limit h_ice by its maximum value
1134h_ice_n_flk = min(h_ice_n_flk, h_ice_max)
1135
1136! Security, limit the ice and snow temperatures by the freezing point
1137t_snow_n_flk = min(t_snow_n_flk, tpl_t_f)
1138t_ice_n_flk = min(t_ice_n_flk, tpl_t_f)
1139
1140!_tmp
1141! Security, avoid too low values (these constraints are used for debugging purposes)
1142! T_snow_n_flk = MAX(T_snow_n_flk, 73.15)
1143! T_ice_n_flk = MAX(T_ice_n_flk, 73.15)
1144! The lowest natural temperature ever
1145! directly recorded at ground level on Earth is -89.2 C (-128.6 F; 184.0 K)
1146! at the Soviet Vostok Station in Antarctica on 21 July, 1983 by ground
1147! measurements.
1148! https://en.wikipedia.org/wiki/Lowest_temperature_recorded_on_Earth
1149
1150 t_snow_n_flk = max(t_snow_n_flk, 184.0)
1151 t_ice_n_flk = max(t_ice_n_flk, 184.0)
1152!_tmp
1153
1154! Remove too thin ice and/or snow
1155IF(h_ice_n_flk.LT.h_ice_min_flk) THEN ! Check ice
1156 h_ice_n_flk = 0.0 ! Ice is too thin, remove it, and
1157 t_ice_n_flk = tpl_t_f ! set T_ice to the freezing point.
1158 h_snow_n_flk = 0.0 ! Remove snow when there is no ice, and
1159 t_snow_n_flk = tpl_t_f ! set T_snow to the freezing point.
1160 l_ice_create = .false. ! "Exotic" case, ice has been created but proved to be too thin
1161ELSE IF(h_snow_n_flk.LT.h_snow_min_flk) THEN ! Ice exists, check snow
1162 h_snow_n_flk = 0.0 ! Snow is too thin, remove it,
1163 t_snow_n_flk = t_ice_n_flk ! and set the snow temperature equal to the ice temperature.
1164END IF
1165
1166
1167!------------------------------------------------------------------------------
1168! Compute the mean temperature of the water column.
1169!------------------------------------------------------------------------------
1170
1171IF(l_ice_create) q_w_flk = 0.0 ! Ice has just been created, set Q_w to zero
1172d_t_mnw_dt = (q_w_flk - q_bot_flk + i_w_flk - i_bot_flk)/tpl_rho_w_r/tpl_c_w/depth_w
1173
1174!print*,'d_T_mnw_dt= ',d_T_mnw_dt
1175!print*,'Q_w_flk=',Q_w_flk
1176!print*,'Q_bot_flk=',Q_bot_flk
1177!print*,'I_w_flk=',I_w_flk
1178!print*,'I_bot_flk=',I_bot_flk
1179!print*,'tpl_rho_w_r=',tpl_rho_w_r
1180!print*,'tpl_c_w=',tpl_c_w
1181!print*,'depth_w=',depth_w
1182
1183t_mnw_n_flk = t_mnw_p_flk + d_t_mnw_dt*del_time ! Advance T_mnw
1184
1185t_mnw_n_flk = max(t_mnw_n_flk, tpl_t_f) ! Limit T_mnw by the freezing point
1186
1187
1188!------------------------------------------------------------------------------
1189! Compute the mixed-layer depth, the mixed-layer temperature,
1190! the bottom temperature and the shape factor
1191! with respect to the temperature profile in the thermocline.
1192! Different formulations are used, depending on the regime of mixing.
1193!------------------------------------------------------------------------------
1194
1195htc_water: IF(h_ice_n_flk.GE.h_ice_min_flk) THEN ! Ice exists
1196
1197 t_mnw_n_flk = min(t_mnw_n_flk, tpl_t_r) ! Limit the mean temperature under the ice by T_r
1198
1199 t_wml_n_flk = tpl_t_f ! The mixed-layer temperature is equal to the freezing point
1200
1201 IF(l_ice_create) THEN ! Ice has just been created
1202 IF(h_ml_p_flk.GE.depth_w-h_ml_min_flk) THEN ! h_ML=D when ice is created
1203 h_ml_n_flk = 0.0 ! Set h_ML to zero
1204 c_t_n_flk = c_t_min ! Set C_T to its minimum value
1205 ELSE ! h_ML<D when ice is created
1206 h_ml_n_flk = h_ml_p_flk ! h_ML remains unchanged
1207 c_t_n_flk = c_t_p_flk ! C_T (thermocline) remains unchanged
1208 END IF
1209 t_bot_n_flk = t_wml_n_flk - (t_wml_n_flk-t_mnw_n_flk)/c_t_n_flk/(1.0-h_ml_n_flk/depth_w)
1210 ! Update the bottom temperature
1211
1212 ELSE IF(t_bot_p_flk.LT.tpl_t_r) THEN ! Ice exists and T_bot < T_r, molecular heat transfer
1213 h_ml_n_flk = h_ml_p_flk ! h_ML remains unchanged
1214 c_t_n_flk = c_t_p_flk ! C_T (thermocline) remains unchanged
1215 t_bot_n_flk = t_wml_n_flk - (t_wml_n_flk-t_mnw_n_flk)/c_t_n_flk/(1.0-h_ml_n_flk/depth_w)
1216 ! Update the bottom temperature
1217
1218 ELSE ! Ice exists and T_bot = T_r, convection due to bottom heating
1219 t_bot_n_flk = tpl_t_r ! T_bot is equal to the temperature of maximum density
1220 IF(h_ml_p_flk.GE.c_small_flk) THEN ! h_ML > 0
1221 c_t_n_flk = c_t_p_flk ! C_T (thermocline) remains unchanged
1222 h_ml_n_flk = depth_w*(1.0-(t_wml_n_flk-t_mnw_n_flk)/(t_wml_n_flk-t_bot_n_flk)/c_t_n_flk)
1223 h_ml_n_flk = max(h_ml_n_flk, 0.0) ! Update the mixed-layer depth
1224 ELSE ! h_ML = 0
1225 h_ml_n_flk = h_ml_p_flk ! h_ML remains unchanged
1226 c_t_n_flk = (t_wml_n_flk-t_mnw_n_flk)/(t_wml_n_flk-t_bot_n_flk)
1227 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min)) ! Update the shape factor (thermocline)
1228 END IF
1229 END IF
1230
1231 t_bot_n_flk = min(t_bot_n_flk, tpl_t_r) ! Security, limit the bottom temperature by T_r
1232
1233ELSE htc_water ! Open water
1234
1235! Generalised buoyancy flux scale and convective velocity scale
1236 flk_str_1 = flake_buoypar(t_wml_p_flk)*q_star_flk/tpl_rho_w_r/tpl_c_w
1237 IF(flk_str_1.LT.0.0) THEN
1238 w_star_sfc_flk = (-flk_str_1*h_ml_p_flk)**(1.0/3.0) ! Convection
1239 ELSE
1240 w_star_sfc_flk = 0.0 ! Neutral or stable stratification
1241 END IF
1242
1243!_dm
1244! The equilibrium depth of the CBL due to surface cooling with the volumetric heating
1245! is not computed as a solution to the transcendental equation.
1246! Instead, an algebraic formula is used
1247! that interpolates between the two asymptotic limits.
1248!_dm
1249 conv_equil_h_scale = -q_w_flk/max(i_w_flk, c_small_flk)
1250 IF(conv_equil_h_scale.GT.0.0 .AND. conv_equil_h_scale.LT.1.0 &
1251 .AND. t_wml_p_flk.GT.tpl_t_r) THEN ! The equilibrium CBL depth scale is only used above T_r
1252 conv_equil_h_scale = sqrt(6.0*conv_equil_h_scale) &
1253 + 2.0*conv_equil_h_scale/(1.0-conv_equil_h_scale)
1254 conv_equil_h_scale = min(depth_w, conv_equil_h_scale/extincoef_water_typ)
1255! print*,'extincoef_water_typ=',extincoef_water_typ
1256 ELSE
1257 conv_equil_h_scale = 0.0 ! Set the equilibrium CBL depth to zero
1258 END IF
1259
1260! Mean buoyancy frequency in the thermocline
1261 n_t_mean = flake_buoypar(0.5*(t_wml_p_flk+t_bot_p_flk))*(t_wml_p_flk-t_bot_p_flk)
1262 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk) THEN
1263 tmp = n_t_mean/(depth_w-h_ml_p_flk)
1264 n_t_mean = sqrt(abs(tmp)) ! Compute N
1265 ELSE
1266 n_t_mean = 0.0 ! h_ML=D, set N to zero
1267 END IF
1268
1269
1270! The rate of change of C_T
1271 d_c_t_dt = max(w_star_sfc_flk, u_star_w_flk, u_star_min_flk)**2_iintegers
1272 d_c_t_dt = n_t_mean*(depth_w-h_ml_p_flk)**2_iintegers &
1273 / c_relax_c/d_c_t_dt ! Relaxation time scale for C_T
1274 d_c_t_dt = (c_t_max-c_t_min)/max(d_c_t_dt, c_small_flk) ! Rate-of-change of C_T
1275
1276! Compute the shape factor and the mixed-layer depth,
1277! using different formulations for convection and wind mixing
1278
1279 c_tt_flk = c_tt_1*c_t_p_flk-c_tt_2 ! C_TT, using C_T at the previous time step
1280 c_q_flk = 2.0*c_tt_flk/c_t_p_flk ! C_Q using C_T at the previous time step
1281
1282 mixing_regime: IF(flk_str_1.LT.0.0) THEN ! Convective mixing
1283
1284 c_t_n_flk = c_t_p_flk + d_c_t_dt*del_time ! Update C_T, assuming dh_ML/dt>0
1285 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min)) ! Limit C_T
1286 d_c_t_dt = (c_t_n_flk-c_t_p_flk)/del_time ! Re-compute dC_T/dt
1287
1288 IF(h_ml_p_flk.LE.depth_w-h_ml_min_flk) THEN ! Compute dh_ML/dt
1289 IF(h_ml_p_flk.LE.h_ml_min_flk) THEN ! Use a reduced entrainment equation (spin-up)
1290 d_h_ml_dt = c_cbl_1/c_cbl_2*max(w_star_sfc_flk, c_small_flk)
1291
1292!_dbg
1293! print*, ' FLake: reduced entrainment eq. D_time*d_h_ML_dt = ', d_h_ML_dt*del_time
1294! print*, ' w_* = ', w_star_sfc_flk
1295! print*, ' \beta*Q_* = ', flk_str_1
1296!_dbg
1297
1298 ELSE ! Use a complete entrainment equation
1299 r_h_icesnow = depth_w/h_ml_p_flk
1300 r_rho_c_icesnow = r_h_icesnow-1.0
1301 r_ti_icesnow = c_t_p_flk/c_tt_flk
1302 r_tstar_icesnow = (r_ti_icesnow/2.0-1.0)*r_rho_c_icesnow + 1.0
1303 d_h_ml_dt = -q_star_flk*(r_tstar_icesnow*(1.0+c_cbl_1)-1.0) - q_bot_flk
1304 d_h_ml_dt = d_h_ml_dt/tpl_rho_w_r/tpl_c_w ! Q_* and Q_b flux terms
1305 flk_str_2 = (depth_w-h_ml_p_flk)*(t_wml_p_flk-t_bot_p_flk)*c_tt_2/c_tt_flk*d_c_t_dt
1306 d_h_ml_dt = d_h_ml_dt + flk_str_2 ! Add dC_T/dt term
1307 flk_str_2 = i_bot_flk + (r_ti_icesnow-1.0)*i_h_flk - r_ti_icesnow*i_intm_h_d_flk
1308 flk_str_2 = flk_str_2 + (r_ti_icesnow-2.0)*r_rho_c_icesnow*(i_h_flk-i_intm_0_h_flk)
1309 flk_str_2 = flk_str_2/tpl_rho_w_r/tpl_c_w
1310 d_h_ml_dt = d_h_ml_dt + flk_str_2 ! Add radiation terms
1311 flk_str_2 = -c_cbl_2*r_tstar_icesnow*q_star_flk/tpl_rho_w_r/tpl_c_w/max(w_star_sfc_flk, c_small_flk)
1312 flk_str_2 = flk_str_2 + c_t_p_flk*(t_wml_p_flk-t_bot_p_flk)
1313 d_h_ml_dt = d_h_ml_dt/flk_str_2 ! dh_ML/dt = r.h.s.
1314 END IF
1315!_dm
1316! Notice that dh_ML/dt may appear to be negative
1317! (e.g. due to buoyancy loss to bottom sediments and/or
1318! the effect of volumetric radiation heating),
1319! although a negative generalized buoyancy flux scale indicates
1320! that the equilibrium CBL depth has not yet been reached
1321! and convective deepening of the mixed layer should take place.
1322! Physically, this situation reflects an approximate character of the lake model.
1323! Using the self-similar temperature profile in the thermocline,
1324! there is always communication between the mixed layer, the thermocline
1325! and the lake bottom. As a result, the rate of change of the CBL depth
1326! is always dependent on the bottom heat flux and the radiation heating of the thermocline.
1327! In reality, convective mixed-layer deepening may be completely decoupled
1328! from the processes underneath. In order to account for this fact,
1329! the rate of CBL deepening is set to a small value
1330! if dh_ML/dt proves to be negative.
1331! This is "double insurance" however,
1332! as a negative dh_ML/dt is encountered very rarely.
1333!_dm
1334
1335!_dbg
1336! IF(d_h_ML_dt.LT.0.0) THEN
1337! print*, 'FLake: negative d_h_ML_dt during convection, = ', d_h_ML_dt
1338! print*, ' d_h_ML_dt*del_time = ', MAX(d_h_ML_dt, c_small_flk)*del_time
1339! print*, ' u_* = ', u_star_w_flk
1340! print*, ' w_* = ', w_star_sfc_flk
1341! print*, ' h_CBL_eqi = ', conv_equil_h_scale
1342! print*, ' ZM scale = ', ZM_h_scale
1343! print*, ' h_ML_p_flk = ', h_ML_p_flk
1344! END IF
1345! print*, 'FLake: Convection, = ', d_h_ML_dt
1346! print*, ' Q_* = ', Q_star_flk
1347! print*, ' \beta*Q_* = ', flk_str_1
1348!_dbg
1349
1350 d_h_ml_dt = max(d_h_ml_dt, c_small_flk)
1351 h_ml_n_flk = h_ml_p_flk + d_h_ml_dt*del_time ! Update h_ML
1352 h_ml_n_flk = max(h_ml_min_flk, min(h_ml_n_flk, depth_w)) ! Security, limit h_ML
1353 ELSE ! Mixing down to the lake bottom
1354 h_ml_n_flk = depth_w
1355 END IF
1356
1357 ELSE mixing_regime ! Wind mixing
1358
1359 d_h_ml_dt = max(u_star_w_flk, u_star_min_flk) ! The surface friction velocity
1360 zm_h_scale = (abs(par_coriolis)/c_sbl_zm_n + n_t_mean/c_sbl_zm_i)*d_h_ml_dt**2_iintegers
1361 zm_h_scale = zm_h_scale + flk_str_1/c_sbl_zm_s
1362 zm_h_scale = max(zm_h_scale, c_small_flk)
1363 zm_h_scale = d_h_ml_dt**3_iintegers/zm_h_scale
1364 zm_h_scale = max(h_ml_min_flk, min(zm_h_scale, h_ml_max_flk)) ! The ZM96 SBL depth scale
1365 zm_h_scale = max(zm_h_scale, conv_equil_h_scale) ! Equilibrium mixed-layer depth
1366
1367!_dm
1368! In order to avoid numerical discretization problems,
1369! an analytical solution to the evolution equation
1370! for the wind-mixed layer depth is used.
1371! That is, an exponential relaxation formula is applied
1372! over the time interval equal to the model time step.
1373!_dm
1374
1375 d_h_ml_dt = c_relax_h*d_h_ml_dt/zm_h_scale*del_time
1376 h_ml_n_flk = zm_h_scale - (zm_h_scale-h_ml_p_flk)*exp(-d_h_ml_dt) ! Update h_ML
1377 h_ml_n_flk = max(h_ml_min_flk, min(h_ml_n_flk, depth_w)) ! Limit h_ML
1378 d_h_ml_dt = (h_ml_n_flk-h_ml_p_flk)/del_time ! Re-compute dh_ML/dt
1379
1380 IF(h_ml_n_flk.LE.h_ml_p_flk) &
1381 d_c_t_dt = -d_c_t_dt ! Mixed-layer retreat or stationary state, dC_T/dt<0
1382 c_t_n_flk = c_t_p_flk + d_c_t_dt*del_time ! Update C_T
1383 c_t_n_flk = min(c_t_max, max(c_t_n_flk, c_t_min)) ! Limit C_T
1384 d_c_t_dt = (c_t_n_flk-c_t_p_flk)/del_time ! Re-compute dC_T/dt
1385
1386!_dbg
1387! print*, 'FLake: wind mixing: d_h_ML_dt*del_time = ', d_h_ML_dt*del_time
1388! print*, ' h_CBL_eqi = ', conv_equil_h_scale
1389! print*, ' ZM scale = ', ZM_h_scale
1390! print*, ' w_* = ', w_star_sfc_flk
1391! print*, ' u_* = ', u_star_w_flk
1392! print*, ' h_ML_p_flk = ', h_ML_p_flk
1393!_dbg
1394
1395 END IF mixing_regime
1396
1397! Compute the time-rate-of-change of the the bottom temperature,
1398! depending on the sign of dh_ML/dt
1399! Update the bottom temperature and the mixed-layer temperature
1400
1401 IF(h_ml_n_flk.LE.depth_w-h_ml_min_flk) THEN ! Mixing did not reach the bottom
1402
1403 IF(h_ml_n_flk.GT.h_ml_p_flk) THEN ! Mixed-layer deepening
1404 r_h_icesnow = h_ml_p_flk/depth_w
1405 r_rho_c_icesnow = 1.0-r_h_icesnow
1406 r_ti_icesnow = 0.5*c_t_p_flk*r_rho_c_icesnow+c_tt_flk*(2.0*r_h_icesnow-1.0)
1407 r_tstar_icesnow = (0.5+c_tt_flk-c_q_flk)/r_ti_icesnow
1408 r_ti_icesnow = (1.0-c_t_p_flk*r_rho_c_icesnow)/r_ti_icesnow
1409
1410 d_t_bot_dt = (q_w_flk-q_bot_flk+i_w_flk-i_bot_flk)/tpl_rho_w_r/tpl_c_w
1411 d_t_bot_dt = d_t_bot_dt - c_t_p_flk*(t_wml_p_flk-t_bot_p_flk)*d_h_ml_dt
1412 d_t_bot_dt = d_t_bot_dt*r_tstar_icesnow/depth_w ! Q+I fluxes and dh_ML/dt term
1413
1414 flk_str_2 = i_intm_h_d_flk - (1.0-c_q_flk)*i_h_flk - c_q_flk*i_bot_flk
1415 flk_str_2 = flk_str_2*r_ti_icesnow/(depth_w-h_ml_p_flk)/tpl_rho_w_r/tpl_c_w
1416 d_t_bot_dt = d_t_bot_dt + flk_str_2 ! Add radiation-flux term
1417
1418 flk_str_2 = (1.0-c_tt_2*r_ti_icesnow)/c_t_p_flk
1419 flk_str_2 = flk_str_2*(t_wml_p_flk-t_bot_p_flk)*d_c_t_dt
1420 d_t_bot_dt = d_t_bot_dt + flk_str_2 ! Add dC_T/dt term
1421
1422 ELSE ! Mixed-layer retreat or stationary state
1423 d_t_bot_dt = 0.0 ! dT_bot/dt=0
1424 END IF
1425
1426 t_bot_n_flk = t_bot_p_flk + d_t_bot_dt*del_time ! Update T_bot
1427 t_bot_n_flk = max(t_bot_n_flk, tpl_t_f) ! Security, limit T_bot by the freezing point
1428 flk_str_2 = (t_bot_n_flk-tpl_t_r)*flake_buoypar(t_mnw_n_flk)
1429 IF(flk_str_2.LT.0.0) t_bot_n_flk = tpl_t_r ! Security, avoid T_r crossover
1430 t_wml_n_flk = c_t_n_flk*(1.0-h_ml_n_flk/depth_w)
1431 t_wml_n_flk = (t_mnw_n_flk-t_bot_n_flk*t_wml_n_flk)/(1.0-t_wml_n_flk)
1432 t_wml_n_flk = max(t_wml_n_flk, tpl_t_f) ! Security, limit T_wML by the freezing point
1433! print*,'mark 2 T_wML_n_flk=',T_wML_n_flk
1434
1435 ELSE ! Mixing down to the lake bottom
1436
1437 h_ml_n_flk = depth_w
1438 t_wml_n_flk = t_mnw_n_flk
1439 t_bot_n_flk = t_mnw_n_flk
1440 c_t_n_flk = c_t_min
1441! print*,'mark 3 T_wML_n_flk=',T_wML_n_flk
1442
1443 END IF
1444
1445END IF htc_water
1446
1447
1448!------------------------------------------------------------------------------
1449! Compute the depth of the upper layer of bottom sediments
1450! and the temperature at that depth.
1451!------------------------------------------------------------------------------
1452
1453use_sediment: IF(lflk_botsed_use) THEN ! The bottom-sediment scheme is used
1454
1455 IF(h_b1_p_flk.GE.depth_bs-h_b1_min_flk) THEN ! No T(z) maximum (no thermal wave)
1456 h_b1_p_flk = 0.0 ! Set H_B1_p to zero
1457 t_b1_p_flk = t_bot_p_flk ! Set T_B1_p to the bottom temperature
1458 END IF
1459
1460 flk_str_1 = 2.0*phi_b1_pr0/(1.0-c_b1)*tpl_kappa_w/tpl_rho_w_r/tpl_c_w*del_time
1461 h_ice_threshold = sqrt(flk_str_1) ! Threshold value of H_B1
1462 h_ice_threshold = min(0.9*depth_bs, h_ice_threshold) ! Limit H_B1
1463 flk_str_2 = c_b2/(1.0-c_b2)*(t_bs-t_b1_p_flk)/(depth_bs-h_b1_p_flk)
1464
1465 IF(h_b1_p_flk.LT.h_ice_threshold) THEN ! Use a truncated equation for H_B1(t)
1466 h_b1_n_flk = sqrt(h_b1_p_flk**2_iintegers+flk_str_1) ! Advance H_B1
1467 d_h_b1_dt = (h_b1_n_flk-h_b1_p_flk)/del_time ! Re-compute dH_B1/dt
1468 ELSE ! Use a full equation for H_B1(t)
1469 flk_str_1 = (q_bot_flk+i_bot_flk)/h_b1_p_flk/tpl_rho_w_r/tpl_c_w
1470 flk_str_1 = flk_str_1 - (1.0-c_b1)*(t_bot_n_flk-t_bot_p_flk)/del_time
1471 d_h_b1_dt = (1.0-c_b1)*(t_bot_p_flk-t_b1_p_flk)/h_b1_p_flk + c_b1*flk_str_2
1472 d_h_b1_dt = flk_str_1/d_h_b1_dt
1473 h_b1_n_flk = h_b1_p_flk + d_h_b1_dt*del_time ! Advance H_B1
1474 END IF
1475 d_t_b1_dt = flk_str_2*d_h_b1_dt
1476 t_b1_n_flk = t_b1_p_flk + d_t_b1_dt*del_time ! Advance T_B1
1477
1478!_dbg
1479! print*, 'BS module: '
1480! print*, ' Q_bot = ', Q_bot_flk
1481! print*, ' d_H_B1_dt = ', d_H_B1_dt
1482! print*, ' d_T_B1_dt = ', d_T_B1_dt
1483! print*, ' H_B1 = ', H_B1_n_flk
1484! print*, ' T_bot = ', T_bot_n_flk
1485! print*, ' T_B1 = ', T_B1_n_flk
1486! print*, ' T_bs = ', T_bs
1487!_dbg
1488
1489!_nu
1490! Use a very simplistic procedure, where only the upper layer profile is used,
1491! H_B1 is always set to depth_bs, and T_B1 is always set to T_bs.
1492! Then, the time derivatives are zero, and the sign of the bottom heat flux depends on
1493! whether T_bot is smaller or greater than T_bs.
1494! This is, of course, an oversimplified scheme.
1495!_nu d_H_B1_dt = 0.0
1496!_nu d_T_B1_dt = 0.0
1497!_nu H_B1_n_flk = H_B1_p_flk + d_H_B1_dt*del_time ! Advance H_B1
1498!_nu T_B1_n_flk = T_B1_p_flk + d_T_B1_dt*del_time ! Advance T_B1
1499!_nu
1500
1501 l_snow_exists = h_b1_n_flk.GE.depth_bs-h_b1_min_flk & ! H_B1 reached depth_bs, or
1502 .OR. h_b1_n_flk.LT.h_b1_min_flk & ! H_B1 decreased to zero, or
1503 .OR.(t_bot_n_flk-t_b1_n_flk)*(t_bs-t_b1_n_flk).LE.0.0 ! there is no T(z) maximum
1504 IF(l_snow_exists) THEN
1505 h_b1_n_flk = depth_bs ! Set H_B1 to the depth of the thermally active layer
1506 t_b1_n_flk = t_bs ! Set T_B1 to the climatological temperature
1507 END IF
1508
1509! changed by Shaobo Zhang
1510
1511!ELSE Use_sediment ! The bottom-sediment scheme is not used
1512!
1513! H_B1_n_flk = rflk_depth_bs_ref ! H_B1 is set to a reference value
1514! T_B1_n_flk = tpl_T_r ! T_B1 is set to the temperature of maximum density
1515
1516ELSE use_sediment ! The scheme written by Shaobo Zhang for the deeper layer of a deep lake is used
1517
1518 if ( abs(t_bot_p_flk-t_bot_2_in) .lt. 0.01 ) then
1519 depth_bs = depthbs
1520 depth_w = depthw
1521 t_bot_2_out = t_bot_2_in
1522 else
1523
1524 ct = 0.65
1525 ctt = 11.0/18.0 * ct - 7.0/45.0
1526 cq = 2.0 * ctt / ct
1527
1528! compute d_h_D_dt
1529 flk_str_1 = ( q_bot_flk * cq + i_bot_flk - i_intm_d_h_flk/depth_bs ) /tpl_rho_w_r/tpl_c_w
1530 flk_str_1 = flk_str_1 - depth_bs * (0.5-ctt) * (t_bot_n_flk-t_bot_p_flk)/del_time
1531 flk_str_1 = flk_str_1 - ctt/ct*( (q_bot_flk+i_bot_flk-i_hh_flk)/tpl_rho_w_r/tpl_c_w - &
1532 depth_bs * ( 1.0 - ct ) * (t_bot_n_flk-t_bot_p_flk)/del_time )
1533 flk_str_2 = ctt * (t_bot_p_flk-t_bot_2_in)
1534 if(abs(flk_str_2)<0.01) then
1535 d_h_d_dt = 0.0
1536 else
1537 d_h_d_dt = flk_str_1/flk_str_2
1538 endif
1539
1540! compute d_T_H_dt
1541 flk_str_1 = (q_bot_flk+i_bot_flk-i_hh_flk)/tpl_rho_w_r/tpl_c_w
1542 flk_str_1 = flk_str_1 - depth_bs*(1.0-ct)*(t_bot_n_flk-t_bot_p_flk)/del_time
1543 flk_str_1 = flk_str_1 - ct * (t_bot_p_flk-t_bot_2_in) * d_h_d_dt
1544 flk_str_2 = ct * depth_bs
1545 d_t_h_dt = flk_str_1/flk_str_2
1546
1547! update T_bot_2_out
1548 t_bot_2_out = t_bot_2_in + d_t_h_dt * del_time
1549
1550 if ( (t_bot_2_out-tpl_t_r)*(t_bot_n_flk-tpl_t_r) .le. 0.0 ) then
1551 t_bot_2_out = tpl_t_r
1552 else if ( (t_bot_n_flk-tpl_t_r)/(t_bot_2_out-tpl_t_r) .le. 1.0 ) then
1553 t_bot_2_out = t_bot_n_flk
1554 end if
1555
1556! update depth_w
1557 flk_str_2 = depth_w + depth_bs ! The total depth of the deep lake
1558
1559 depth_w = depth_w - d_h_d_dt * del_time ! Advance depth_bs
1560! depth_w = max (lake_depth_max, min(depth_w, 0.4*flk_str_2)) ! Limit depth_bs
1561 depth_w = max(max(lake_depth_max,0.3*flk_str_2), &
1562 min(depth_w, min(0.4*flk_str_2, 80.0))) ! Limit depth_bs
1563 depth_bs = flk_str_2 - depth_w ! Update depth_w
1564 end if
1565
1566END IF use_sediment
1567
1568
1569!------------------------------------------------------------------------------
1570! Impose additional constraints.
1571!------------------------------------------------------------------------------
1572
1573! In case of unstable stratification, force mixing down to the bottom
1574flk_str_2 = (t_wml_n_flk-t_bot_n_flk)*flake_buoypar(t_mnw_n_flk)
1575IF(flk_str_2.LT.0.0) THEN
1576
1577!_dbg
1578! print*, 'FLake: inverse (unstable) stratification !!! '
1579! print*, ' Mixing down to the bottom is forced.'
1580! print*, ' T_wML_p, T_wML_n ', T_wML_p_flk-tpl_T_f, T_wML_n_flk-tpl_T_f
1581! print*, ' T_mnw_p, T_mnw_n ', T_mnw_p_flk-tpl_T_f, T_mnw_n_flk-tpl_T_f
1582! print*, ' T_bot_p, T_bot_n ', T_bot_p_flk-tpl_T_f, T_bot_n_flk-tpl_T_f
1583! print*, ' h_ML_p, h_ML_n ', h_ML_p_flk, h_ML_n_flk
1584! print*, ' C_T_p, C_T_n ', C_T_p_flk, C_T_n_flk
1585!_dbg
1586
1587 h_ml_n_flk = depth_w
1588 t_wml_n_flk = t_mnw_n_flk
1589 t_bot_n_flk = t_mnw_n_flk
1590 c_t_n_flk = c_t_min
1591! print*,'mark 4 T_wML_n_flk=',T_wML_n_flk
1592
1593END IF
1594
1595
1596!------------------------------------------------------------------------------
1597! Update the surface temperature.
1598!------------------------------------------------------------------------------
1599
1600IF(h_snow_n_flk.GE.h_snow_min_flk) THEN
1601 t_sfc_n = t_snow_n_flk ! Snow exists, use the snow temperature
1602ELSE IF(h_ice_n_flk.GE.h_ice_min_flk) THEN
1603 t_sfc_n = t_ice_n_flk ! Ice exists but there is no snow, use the ice temperature
1604ELSE
1605 t_sfc_n = t_wml_n_flk ! No ice-snow cover, use the mixed-layer temperature
1606END IF
1607if(t_sfc_n.lt.200.0 .or. t_sfc_n.gt.350.0) then
1608 print*,'h_snow_n_flk=',h_snow_n_flk,' h_Snow_min_flk=',h_snow_min_flk
1609 print*,'h_ice_n_flk=',h_ice_n_flk, ' h_Ice_min_flk= ',h_ice_min_flk
1610 print*,'T_wML_n_flk=',t_wml_n_flk
1611 print*,'T_sfc_n=',t_sfc_n
1612endif
1613!------------------------------------------------------------------------------
1614! End calculations
1615!==============================================================================
1616
1617END SUBROUTINE flake_main
1618
1619!==============================================================================
1620
1621!==============================================================================
1622! include 'flake_buoypar.incf'
1623!------------------------------------------------------------------------------
1624
1625REAL (KIND = ireals) FUNCTION flake_buoypar (T_water)
1626
1627!------------------------------------------------------------------------------
1628!
1629! Description:
1630!
1631! Computes the buoyancy parameter,
1632! using a quadratic equation of state for the fresh-water.
1633!
1634! Declarations:
1635!
1636! Modules used:
1637
1638!_dm Parameters are USEd in module "flake".
1639!_nu USE data_parameters , ONLY : &
1640!_nu ireals, & ! KIND-type parameter for real variables
1641!_nu iintegers ! KIND-type parameter for "normal" integer variables
1642
1643USE flake_parameters , ONLY : &
1644 tpl_grav , & ! Acceleration due to gravity [m s^{-2}]
1645 tpl_t_r , & ! Temperature of maximum density of fresh water [K]
1646 tpl_a_t ! Constant in the fresh-water equation of state [K^{-2}]
1647
1648use machine, only: kind_phys
1649!==============================================================================
1650
1651IMPLICIT NONE
1652
1653!==============================================================================
1654!
1655! Declarations
1656
1657! Input (function argument)
1658real(kind = kind_phys), INTENT(IN) :: &
1659 t_water ! Water temperature [K]
1660
1661!==============================================================================
1662! Start calculations
1663!------------------------------------------------------------------------------
1664
1665! Buoyancy parameter [m s^{-2} K^{-1}]
1666
1667 flake_buoypar = tpl_grav*tpl_a_t*(t_water-tpl_t_r)
1668
1669!------------------------------------------------------------------------------
1670! End calculations
1671!==============================================================================
1672
1673END FUNCTION flake_buoypar
1674
1675!==============================================================================
1676
1677!==============================================================================
1678! include 'flake_snowdensity.incf'
1679!------------------------------------------------------------------------------
1680
1681REAL (KIND = ireals) FUNCTION flake_snowdensity (h_snow)
1682
1683!------------------------------------------------------------------------------
1684!
1685! Description:
1686!
1687! Computes the snow density,
1688! using an empirical approximation from Heise et al. (2003).
1689!
1690! Declarations:
1691!
1692! Modules used:
1693
1694!_dm Parameters are USEd in module "flake".
1695!_nu USE data_parameters , ONLY : &
1696!_nu ireals, & ! KIND-type parameter for real variables
1697!_nu iintegers ! KIND-type parameter for "normal" integer variables
1698
1699USE flake_parameters , ONLY : &
1700 tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
1701 tpl_rho_s_min , & ! Minimum snow density [kg m^{-3}]
1702 tpl_rho_s_max , & ! Maximum snow density [kg m^{-3}]
1703 tpl_gamma_rho_s , & ! Empirical parameter [kg m^{-4}] in the expression for the snow density
1704 c_small_flk ! A small number
1705
1706use machine, only: kind_phys
1707!==============================================================================
1708
1709IMPLICIT NONE
1710
1711!==============================================================================
1712!
1713! Declarations
1714
1715! Input (function argument)
1716real(kind = kind_phys), INTENT(IN) :: &
1717 h_snow ! Snow thickness [m]
1718
1719!==============================================================================
1720! Start calculations
1721!------------------------------------------------------------------------------
1722
1723! Snow density [kg m^{-3}]
1724
1725! Security. Ensure that the expression in () does not become negative at a very large h_snow.
1726 flake_snowdensity = max( c_small_flk, (1.0 - h_snow*tpl_gamma_rho_s/tpl_rho_w_r) )
1727 flake_snowdensity = min( tpl_rho_s_max, tpl_rho_s_min/flake_snowdensity )
1728
1729!------------------------------------------------------------------------------
1730! End calculations
1731!==============================================================================
1732
1733END FUNCTION flake_snowdensity
1734
1735!==============================================================================
1736
1737!==============================================================================
1738! include 'flake_snowheatconduct.incf'
1739!------------------------------------------------------------------------------
1740
1741REAL (KIND = ireals) FUNCTION flake_snowheatconduct (h_snow)
1742
1743!------------------------------------------------------------------------------
1744!
1745! Description:
1746!
1747! Computes the snow heat conductivity,
1748! using an empirical approximation from Heise et al. (2003).
1749!
1750! Declarations:
1751!
1752! Modules used:
1753
1754!_dm Parameters are USEd in module "flake".
1755!_nu USE data_parameters , ONLY : &
1756!_nu ireals, & ! KIND-type parameter for real variables
1757!_nu iintegers ! KIND-type parameter for "normal" integer variables
1758
1759USE flake_parameters , ONLY : &
1760 tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
1761 tpl_kappa_s_min , & ! Minimum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
1762 tpl_kappa_s_max , & ! Maximum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
1763 tpl_gamma_kappa_s ! Empirical parameter [J m^{-2} s^{-1} K^{-1}] in the expression for kappa_S
1764
1765use machine, only: kind_phys
1766!==============================================================================
1767
1768IMPLICIT NONE
1769
1770!==============================================================================
1771!
1772! Declarations
1773
1774! Input (function argument)
1775real(kind = kind_phys), INTENT(IN) :: &
1776 h_snow ! Snow thickness [m]
1777
1778!==============================================================================
1779! Start calculations
1780!------------------------------------------------------------------------------
1781
1782! Snow heat conductivity [J m^{-1} s^{-1} K^{-1} = kg m s^{-3} K^{-1}]
1783
1784 flake_snowheatconduct = flake_snowdensity( h_snow ) ! Compute snow density
1785 flake_snowheatconduct = min( tpl_kappa_s_max, tpl_kappa_s_min &
1786 + h_snow*tpl_gamma_kappa_s*flake_snowheatconduct/tpl_rho_w_r )
1787
1788!------------------------------------------------------------------------------
1789! End calculations
1790!==============================================================================
1791
1792END FUNCTION flake_snowheatconduct
1793
1794!==============================================================================
1795
1796END MODULE flake
1797
1798!------------------------------------------------------------------------------
1799
1801
1802!------------------------------------------------------------------------------
1803!
1804! Description:
1805!
1806! The main program unit of
1807! the atmospheric surface-layer parameterization scheme "SfcFlx".
1808! "SfcFlx" is used to compute fluxes
1809! of momentum and of sensible and latent heat over lakes.
1810! The surface-layer scheme developed by Mironov (1991) was used as the starting point.
1811! It was modified and further developed to incorporate new results as to
1812! the roughness lenghts for scalar quantities,
1813! heat and mass transfer in free convection,
1814! and the effect of limited fetch on the momentum transfer.
1815! Apart from the momentum flux and sensible and latent heat fluxes,
1816! the long-wave radiation flux from the water surface and
1817! the long-wave radiation flux from the atmosphere can also be computed.
1818! The atmospheric long-wave radiation flux is computed with simple empirical formulae,
1819! where the atmospheric emissivity is taken to be dependent on
1820! the water vapour pressure and cloud fraction.
1821!
1822! A description of SfcFlx is available from the author.
1823! Dmitrii Mironov
1824! German Weather Service, Kaiserleistr. 29/35, D-63067 Offenbach am Main, Germany.
1825! dmitrii.mironov@dwd.de
1826!
1827! Lines embraced with "!_tmp" contain temporary parts of the code.
1828! Lines embraced/marked with "!_dev" may be replaced
1829! as improved parameterizations are developed and tested.
1830! Lines embraced/marked with "!_dm" are DM's comments
1831! that may be helpful to a user.
1832! Lines embraced/marked with "!_dbg" are used
1833! for debugging purposes only.
1834!
1835
1836
1837USE data_parameters , ONLY : &
1838 ireals , & ! KIND-type parameter for real variables
1839 iintegers ! KIND-type parameter for "normal" integer variables
1840
1841USE flake_parameters , ONLY : &
1842 tpl_grav , & ! Acceleration due to gravity [m s^{-2}]
1843 tpl_t_f , & ! Fresh water freezing point [K]
1844 tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
1845 tpl_c_w , & ! Specific heat of water [J kg^{-1} K^{-1}]
1846 tpl_l_f , & ! Latent heat of fusion [J kg^{-1}]
1847 h_ice_min_flk ! Minimum ice thickness [m]
1848
1849use machine, only: kind_phys
1850
1851!==============================================================================
1852
1853IMPLICIT NONE
1854
1855!==============================================================================
1856!
1857! Declarations
1858
1859! Dimensionless constants in the Monin-Obukhov surface-layer
1860! similarity relations and in the expressions for the roughness lengths.
1861real(kind = kind_phys), PARAMETER :: &
1862 c_karman = 0.40 , & ! The von Karman constant
1863! Pr_neutral = 1.0 , & ! Turbulent Prandtl number at neutral static stability
1864 pr_neutral = 0.9 , & ! Turbulent Prandtl number at neutral static stability
1865 sc_neutral = 1.0 , & ! Turbulent Schmidt number at neutral static stability
1866 c_mo_u_stab = 5.0 , & ! Constant of the MO theory (wind, stable stratification)
1867 c_mo_t_stab = 5.0 , & ! Constant of the MO theory (temperature, stable stratification)
1868 c_mo_q_stab = 5.0 , & ! Constant of the MO theory (humidity, stable stratification)
1869 c_mo_u_conv = 15.0 , & ! Constant of the MO theory (wind, convection)
1870 c_mo_t_conv = 15.0 , & ! Constant of the MO theory (temperature, convection)
1871 c_mo_q_conv = 15.0 , & ! Constant of the MO theory (humidity, convection)
1872 c_mo_u_exp = 0.25 , & ! Constant of the MO theory (wind, exponent)
1873 c_mo_t_exp = 0.5 , & ! Constant of the MO theory (temperature, exponent)
1874 c_mo_q_exp = 0.5 , & ! Constant of the MO theory (humidity, exponent)
1875 z0u_ice_rough = 1.0e-03 , & ! Aerodynamic roughness of the ice surface [m] (rough flow)
1876 c_z0u_smooth = 0.1 , & ! Constant in the expression for z0u (smooth flow)
1877 c_z0u_rough = 1.23e-02 , & ! The Charnock constant in the expression for z0u (rough flow)
1878 c_z0u_rough_l = 1.00e-01 , & ! An increased Charnock constant (used as the upper limit)
1879 c_z0u_ftch_f = 0.70 , & ! Factor in the expression for fetch-dependent Charnock parameter
1880 c_z0u_ftch_ex = 0.3333333 , & ! Exponent in the expression for fetch-dependent Charnock parameter
1881 c_z0t_rough_1 = 4.0 , & ! Constant in the expression for z0t (factor)
1882 c_z0t_rough_2 = 3.2 , & ! Constant in the expression for z0t (factor)
1883 c_z0t_rough_3 = 0.5 , & ! Constant in the expression for z0t (exponent)
1884 c_z0q_rough_1 = 4.0 , & ! Constant in the expression for z0q (factor)
1885 c_z0q_rough_2 = 4.2 , & ! Constant in the expression for z0q (factor)
1886 c_z0q_rough_3 = 0.5 , & ! Constant in the expression for z0q (exponent)
1887 c_z0t_ice_b0s = 1.250 , & ! Constant in the expression for z0t over ice
1888 c_z0t_ice_b0t = 0.149 , & ! Constant in the expression for z0t over ice
1889 c_z0t_ice_b1t = -0.550 , & ! Constant in the expression for z0t over ice
1890 c_z0t_ice_b0r = 0.317 , & ! Constant in the expression for z0t over ice
1891 c_z0t_ice_b1r = -0.565 , & ! Constant in the expression for z0t over ice
1892 c_z0t_ice_b2r = -0.183 , & ! Constant in the expression for z0t over ice
1893 c_z0q_ice_b0s = 1.610 , & ! Constant in the expression for z0q over ice
1894 c_z0q_ice_b0t = 0.351 , & ! Constant in the expression for z0q over ice
1895 c_z0q_ice_b1t = -0.628 , & ! Constant in the expression for z0q over ice
1896 c_z0q_ice_b0r = 0.396 , & ! Constant in the expression for z0q over ice
1897 c_z0q_ice_b1r = -0.512 , & ! Constant in the expression for z0q over ice
1898 c_z0q_ice_b2r = -0.180 , & ! Constant in the expression for z0q over ice
1899 re_z0s_ice_t = 2.5 , & ! Threshold value of the surface Reynolds number
1900 ! used to compute z0t and z0q over ice (Andreas 2002)
1901 re_z0u_thresh = 0.1 ! Threshold value of the roughness Reynolds number
1902 ! [value from Zilitinkevich, Grachev, and Fairall (200),
1903 ! currently not used]
1904
1905! Dimensionless constants
1906real(kind = kind_phys), PARAMETER :: &
1907 c_free_conv = 0.14 ! Constant in the expressions for fluxes in free convection
1908
1909! Dimensionless constants
1910real(kind = kind_phys), PARAMETER :: &
1911 c_lwrad_emis = 0.99 ! Surface emissivity with respect to the long-wave radiation
1912
1913! Thermodynamic parameters
1914real(kind = kind_phys), PARAMETER :: &
1915 tpsf_c_stefboltz = 5.67e-08 , & ! The Stefan-Boltzmann constant [W m^{-2} K^{-4}]
1916 tpsf_r_dryair = 2.8705e+02 , & ! Gas constant for dry air [J kg^{-1} K^{-1}]
1917 tpsf_r_watvap = 4.6151e+02 , & ! Gas constant for water vapour [J kg^{-1} K^{-1}]
1918 tpsf_c_a_p = 1.005e+03 , & ! Specific heat of air at constant pressure [J kg^{-1} K^{-1}]
1919 tpsf_l_evap = 2.501e+06 , & ! Specific heat of evaporation [J kg^{-1}]
1920 tpsf_nu_u_a = 1.50e-05 , & ! Kinematic molecular viscosity of air [m^{2} s^{-1}]
1921 tpsf_kappa_t_a = 2.20e-05 , & ! Molecular temperature conductivity of air [m^{2} s^{-1}]
1922 tpsf_kappa_q_a = 2.40e-05 ! Molecular diffusivity of air for water vapour [m^{2} s^{-1}]
1923
1924! Derived thermodynamic parameters
1925real(kind = kind_phys), PARAMETER :: &
1926 tpsf_rd_o_rv = tpsf_r_dryair/tpsf_r_watvap , & ! Ratio of gas constants (Rd/Rv)
1927 tpsf_alpha_q = (1.0-tpsf_rd_o_rv)/tpsf_rd_o_rv ! Diemsnionless ratio
1928
1929! Thermodynamic parameters
1930real(kind = kind_phys), PARAMETER :: &
1931 p_a_ref = 1.0e+05 ! Reference pressure [N m^{-2} = kg m^{-1} s^{-2}]
1932
1933
1934! The variables declared below
1935! are accessible to all program units of the MODULE "SfcFlx"
1936! and to the driving routines that use "SfcFlx".
1937! These are basically the quantities computed by SfcFlx.
1938! Apart from these quantities, there a few local scalars
1939! used by SfcFlx routines mainly for security reasons.
1940! All variables declared below have a suffix "sf".
1941
1942! SfcFlx variables of type REAL
1943
1944! Roughness lengths
1945real(kind = kind_phys) :: &
1946 z0u_sf , & ! Roughness length with respect to wind velocity [m]
1947 z0t_sf , & ! Roughness length with respect to potential temperature [m]
1948 z0q_sf ! Roughness length with respect to specific humidity [m]
1949
1950! Fluxes in the surface air layer
1951real(kind = kind_phys) :: &
1952 u_star_a_sf , & ! Friction velocity [m s^{-1}]
1953 q_mom_a_sf , & ! Momentum flux [N m^{-2}]
1954 q_sens_a_sf , & ! Sensible heat flux [W m^{-2}]
1955 q_lat_a_sf , & ! Laten heat flux [W m^{-2}]
1956 q_watvap_a_sf ! Flux of water vapout [kg m^{-2} s^{-1}]
1957
1958! Security constants
1959real(kind = kind_phys), PARAMETER :: &
1960 u_wind_min_sf = 1.0e-02 , & ! Minimum wind speed [m s^{-1}]
1961 u_star_min_sf = 1.0e-04 , & ! Minimum value of friction velocity [m s^{-1}]
1962 c_accur_sf = 1.0e-07 , & ! A small number (accuracy)
1963 c_small_sf = 1.0e-04 ! A small number (used to compute fluxes)
1964
1965! Useful constants
1966real(kind = kind_phys), PARAMETER :: &
1967 num_1o3_sf = 1.0/3.0 ! 1/3
1968
1969!==============================================================================
1970! Procedures
1971!==============================================================================
1972
1973CONTAINS
1974
1975!==============================================================================
1976! The codes of the SfcFlx procedures are stored in separate "*.incf" files
1977! and are included below.
1978!------------------------------------------------------------------------------
1979
1980!==============================================================================
1981! include 'SfcFlx_lwradatm.incf'
1982!------------------------------------------------------------------------------
1983
1984REAL (KIND = kind_phys) FUNCTION sfcflx_lwradatm (T_a, e_a, cl_tot, cl_low)
1985
1986!------------------------------------------------------------------------------
1987!
1988! Description:
1989!
1990! Computes the long-wave radiation flux from the atmosphere
1991! as function of air temperature, water vapour pressure and cloud fraction.
1992!
1993! Declarations:
1994!
1995! Modules used:
1996
1997!_dm Parameters are USEd in module "SfcFlx".
1998!_nu USE data_parameters , ONLY : &
1999!_nu ireals, & ! KIND-type parameter for real variables
2000!_nu iintegers ! KIND-type parameter for "normal" integer variables
2001
2002!==============================================================================
2003
2004IMPLICIT NONE
2005
2006!==============================================================================
2007!
2008! Declarations
2009
2010! Input (function argument)
2011real(kind = kind_phys), INTENT(IN) :: &
2012 t_a , & ! Air temperature [K]
2013 e_a , & ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
2014 cl_tot , & ! Total cloud cover [0,1]
2015 cl_low ! Lowe-level cloud cover [0,1]
2016
2017
2018! Local parameters
2019
2020! Coefficients in the empirical formulation
2021! developed at the Main Geophysical Observatory (MGO), St. Petersburg, Russia.
2022real(kind = kind_phys), PARAMETER :: &
2023 c_lmmgo_1 = 43.057924 , & ! Empirical coefficient
2024 c_lmmgo_2 = 540.795 ! Empirical coefficient
2025! Temperature-dependent cloud-correction coefficients in the MGO formula
2026INTEGER (KIND = iintegers), PARAMETER :: &
2027 nband_coef = 6_iintegers ! Number of temperature bands
2028real(kind = kind_phys), PARAMETER, DIMENSION (nband_coef) :: &
2029 corr_cl_tot = (/0.70, 0.45, 0.32, &
2030 0.23, 0.18, 0.13/) , & ! Total clouds
2031 corr_cl_low = (/0.76, 0.49, 0.35, &
2032 0.26, 0.20, 0.15/) , & ! Low-level clouds
2033 corr_cl_midhigh = (/0.46, 0.30, 0.21, &
2034 0.15, 0.12, 0.09/) ! Mid- and high-level clouds
2035real(kind = kind_phys), PARAMETER :: &
2036 t_low = 253.15 , & ! Low-limit temperature in the interpolation formula [K]
2037 del_t = 10.0 ! Temperature step in the interpolation formula [K]
2038
2039! Coefficients in the empirical water-vapour correction function
2040! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001).
2041real(kind = kind_phys), PARAMETER :: &
2042 c_watvap_corr_min = 0.6100 , & ! Empirical coefficient (minimum value of the correction function)
2043 c_watvap_corr_max = 0.7320 , & ! Empirical coefficient (maximum value of the correction function)
2044 c_watvap_corr_e = 0.0050 ! Empirical coefficient [(N m^{-2})^{-1/2}]
2045
2046! Local variables of type INTEGER
2047INTEGER (KIND = iintegers) :: &
2048 i ! Loop index
2049
2050! Local variables of type REAL
2051real(kind = kind_phys) :: &
2052 c_cl_tot_corr , & ! The MGO cloud correction coefficient, total clouds
2053 c_cl_low_corr , & ! The MGO cloud correction coefficient, low-level clouds
2054 c_cl_midhigh_corr , & ! The MGO cloud correction coefficient, mid- and high-level clouds
2055 t_corr , & ! Temperature used to compute the MGO cloud correction [K]
2056 f_wvpres_corr , & ! Correction function with respect to water vapour
2057 f_cloud_corr ! Correction function with respect to cloudiness
2058
2059!==============================================================================
2060! Start calculations
2061!------------------------------------------------------------------------------
2062
2063! Water-vapour correction function
2064 f_wvpres_corr = c_watvap_corr_min + c_watvap_corr_e*sqrt(e_a)
2065 f_wvpres_corr = min(f_wvpres_corr, c_watvap_corr_max)
2066
2067! Cloud-correction coefficients using the MGO formulation with linear interpolation
2068IF(t_a.LT.t_low) THEN
2069 c_cl_tot_corr = corr_cl_tot(1)
2070 c_cl_low_corr = corr_cl_low(1)
2071 c_cl_midhigh_corr = corr_cl_midhigh(1)
2072ELSE IF(t_a.GE.t_low+(nband_coef-1_iintegers)*del_t) THEN
2073 c_cl_tot_corr = corr_cl_tot(nband_coef)
2074 c_cl_low_corr = corr_cl_low(nband_coef)
2075 c_cl_midhigh_corr = corr_cl_midhigh(nband_coef)
2076ELSE
2077 t_corr = t_low
2078 DO i=1, nband_coef-1
2079 IF(t_a.GE.t_corr.AND.t_a.LT.t_corr+del_t) THEN
2080 c_cl_tot_corr = (t_a-t_corr)/del_t
2081 c_cl_low_corr = corr_cl_low(i) + (corr_cl_low(i+1)-corr_cl_low(i))*c_cl_tot_corr
2082 c_cl_midhigh_corr = corr_cl_midhigh(i) + (corr_cl_midhigh(i+1)-corr_cl_midhigh(i))*c_cl_tot_corr
2083 c_cl_tot_corr = corr_cl_tot(i) + (corr_cl_tot(i+1)-corr_cl_tot(i))*c_cl_tot_corr
2084 END IF
2085 t_corr = t_corr + del_t
2086 END DO
2087END IF
2088! Cloud correction function
2089IF(cl_low.LT.0.0) THEN ! Total cloud cover only
2090 f_cloud_corr = 1.0 + c_cl_tot_corr*cl_tot*cl_tot
2091ELSE ! Total and low-level cloud cover
2092 f_cloud_corr = (1.0 + c_cl_low_corr*cl_low*cl_low) &
2093 * (1.0 + c_cl_midhigh_corr*(cl_tot*cl_tot-cl_low*cl_low))
2094END IF
2095
2096! Long-wave radiation flux [W m^{-2}]
2097
2098! The MGO formulation
2099!_nu The MGO formulation
2100!_nu SfcFlx_lwradatm = -SfcFlx_lwradatm*c_lwrad_emis &
2101!_nu * (c_lmMGO_1*SQRT(tpsf_C_StefBoltz*T_a**4_iintegers)-c_lmMGO_2)
2102!_nu
2103
2104! "Conventional" formulation
2105! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001)
2106sfcflx_lwradatm = -c_lwrad_emis*tpsf_c_stefboltz*t_a**4_iintegers &
2107 * f_wvpres_corr*f_cloud_corr
2108
2109!------------------------------------------------------------------------------
2110! End calculations
2111!==============================================================================
2112
2113END FUNCTION sfcflx_lwradatm
2114
2115!==============================================================================
2116
2117!==============================================================================
2118! include 'SfcFlx_lwradwsfc.incf'
2119!------------------------------------------------------------------------------
2120
2121REAL (KIND = kind_phys) FUNCTION sfcflx_lwradwsfc (T)
2122
2123!------------------------------------------------------------------------------
2124!
2125! Description:
2126!
2127! Computes the surface long-wave radiation flux
2128! as function of temperature.
2129!
2130! Declarations:
2131!
2132! Modules used:
2133
2134!_dm Parameters are USEd in module "SfcFlx".
2135!_nu USE data_parameters , ONLY : &
2136!_nu ireals, & ! KIND-type parameter for real variables
2137!_nu iintegers ! KIND-type parameter for "normal" integer variables
2138
2139!==============================================================================
2140
2141IMPLICIT NONE
2142
2143!==============================================================================
2144!
2145! Declarations
2146
2147! Input (function argument)
2148real(kind = kind_phys), INTENT(IN) :: &
2149 t ! Temperature [K]
2150
2151!==============================================================================
2152! Start calculations
2153!------------------------------------------------------------------------------
2154
2155! Long-wave radiation flux [W m^{-2}]
2156
2157sfcflx_lwradwsfc = c_lwrad_emis*tpsf_c_stefboltz*t**4_iintegers
2158
2159!------------------------------------------------------------------------------
2160! End calculations
2161!==============================================================================
2162
2163END FUNCTION sfcflx_lwradwsfc
2164
2165!==============================================================================
2166
2167!==============================================================================
2168! include 'SfcFlx_momsenlat.incf'
2169!------------------------------------------------------------------------------
2170
2171SUBROUTINE sfcflx_momsenlat ( height_u, height_tq, fetch, &
2172 U_a, T_a, q_a, T_s, P_a, h_ice, &
2173 Q_momentum, Q_sensible, Q_latent, Q_watvap,&
2174 q_s,rho_a )
2175
2176!------------------------------------------------------------------------------
2177!
2178! Description:
2179!
2180! The SfcFlx routine
2181! where fluxes of momentum and of sensible and latent heat
2182! at the air-water or air-ice (air-snow) interface are computed.
2183!
2184! Lines embraced with "!_tmp" contain temporary parts of the code.
2185! Lines embraced/marked with "!_dev" may be replaced
2186! as improved parameterizations are developed and tested.
2187! Lines embraced/marked with "!_dm" are DM's comments
2188! that may be helpful to a user.
2189! Lines embraced/marked with "!_dbg" are used
2190! for debugging purposes only.
2191!
2192! Declarations:
2193!
2194! Modules used:
2195
2196!_dm Parameters are USEd in module "SfcFlx".
2197!_nu USE data_parameters , ONLY : &
2198!_nu ireals, & ! KIND-type parameter for real variables
2199!_nu iintegers ! KIND-type parameter for "normal" integer variables
2200
2201use machine, only: kind_phys
2202!==============================================================================
2203
2204IMPLICIT NONE
2205
2206!==============================================================================
2207!
2208! Declarations
2209
2210! Input (procedure arguments)
2211
2212real(kind = kind_phys), INTENT(IN) :: &
2213 height_u , & ! Height where wind is measured [m]
2214 height_tq , & ! Height where temperature and humidity are measured [m]
2215 fetch , & ! Typical wind fetch [m]
2216 u_a , & ! Wind speed [m s^{-1}]
2217 t_a , & ! Air temperature [K]
2218 q_a , & ! Air specific humidity [-]
2219 t_s , & ! Surface temperature (water, ice or snow) [K]
2220 p_a , & ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}]
2221 h_ice ! Ice thickness [m]
2222
2223! Output (procedure arguments)
2224
2225real(kind = kind_phys), INTENT(OUT) :: &
2226 q_momentum , & ! Momentum flux [N m^{-2}]
2227 q_sensible , & ! Sensible heat flux [W m^{-2}]
2228 q_latent , & ! Laten heat flux [W m^{-2}]
2229 q_watvap ! Flux of water vapout [kg m^{-2} s^{-1}]
2230
2231
2232! Local parameters of type INTEGER
2233INTEGER (KIND = iintegers), PARAMETER :: &
2234 n_iter_max = 24 ! Maximum number of iterations
2235
2236! Local variables of type LOGICAL
2237LOGICAL :: &
2238 l_conv_visc , & ! Switch, TRUE = viscous free convection, the Nu=C Ra^(1/3) law is used
2239 l_conv_cbl ! Switch, TRUE = CBL scale convective structures define surface fluxes
2240
2241! Local variables of type INTEGER
2242INTEGER (KIND = iintegers) :: &
2243 i , & ! Loop index
2244 n_iter ! Number of iterations performed
2245
2246! Local variables of type REAL
2247real(kind = kind_phys) :: &
2248 rho_a , & ! Air density [kg m^{-3}]
2249 wvpres_s , & ! Saturation water vapour pressure at T=T_s [N m^{-2}]
2250 q_s ! Saturation specific humidity at T=T_s [-]
2251
2252! Local variables of type REAL
2253real(kind = kind_phys) :: &
2254 q_mom_tur , & ! Turbulent momentum flux [N m^{-2}]
2255 q_sen_tur , & ! Turbulent sensible heat flux [W m^{-2}]
2256 q_lat_tur , & ! Turbulent laten heat flux [W m^{-2}]
2257 q_mom_mol , & ! Molecular momentum flux [N m^{-2}]
2258 q_sen_mol , & ! Molecular sensible heat flux [W m^{-2}]
2259 q_lat_mol , & ! Molecular laten heat flux [W m^{-2}]
2260 q_mom_con , & ! Momentum flux in free convection [N m^{-2}]
2261 q_sen_con , & ! Sensible heat flux in free convection [W m^{-2}]
2262 q_lat_con ! Laten heat flux in free convection [W m^{-2}]
2263
2264! Local variables of type REAL
2265real(kind = kind_phys) :: &
2266 par_conv_visc , & ! Viscous convection stability parameter
2267 par_conv_cbl , & ! CBL convection stability parameter
2268 c_z0u_fetch , & ! Fetch-dependent Charnock parameter
2269 u_a_thresh , & ! Threshld value of the wind speed [m s^{-1}]
2270 u_star_thresh , & ! Threshld value of friction velocity [m s^{-1}]
2271 u_star_previter , & ! Friction velocity from previous iteration [m s^{-1}]
2272 u_star_n , & ! Friction velocity at neutral stratification [m s^{-1}]
2273 u_star_st , & ! Friction velocity with due regard for stratification [m s^{-1}]
2274 zol , & ! The z/L ratio, z=height_u
2275 ri , & ! Gradient Richardson number
2276 ri_cr , & ! Critical value of Ri
2277 r_z , & ! Ratio of "height_tq" to "height_u"
2278 fun , & ! A function of generic variable "x"
2279 fun_prime , & ! Derivative of "Fun" with respect to "x"
2280 delta , & ! Relative error
2281 psi_u , & ! The MO stability function for wind profile
2282 psi_t , & ! The MO stability function for temperature profile
2283 psi_q ! The MO stability function for specific humidity profile
2284
2285
2286!==============================================================================
2287! Start calculations
2288!------------------------------------------------------------------------------
2289
2290!write(*,*) 'Inside Flake Flux Module'
2291!write(*,*) height_u,height_tq,fetch,U_a,T_a,q_a,T_s,P_a,h_ice
2292
2293!_dm All fluxes are positive when directed upwards.
2294
2295!------------------------------------------------------------------------------
2296! Compute saturation specific humidity and the air density at T=T_s
2297!------------------------------------------------------------------------------
2298
2299wvpres_s = sfcflx_satwvpres(t_s, h_ice) ! Saturation water vapour pressure at T=T_s
2300q_s = sfcflx_spechum(wvpres_s, p_a) ! Saturation specific humidity at T=T_s
2301rho_a = sfcflx_rhoair(t_s, q_s, p_a) ! Air density at T_s and q_s (surface values)
2302
2303!------------------------------------------------------------------------------
2304! Compute molecular fluxes of momentum and of sensible and latent heat
2305!------------------------------------------------------------------------------
2306
2307!_dm The fluxes are in kinematic units
2308q_mom_mol = -tpsf_nu_u_a*u_a/height_u
2309q_sen_mol = -tpsf_kappa_t_a*(t_a-t_s)/height_tq
2310q_lat_mol = -tpsf_kappa_q_a*(q_a-q_s)/height_tq
2311
2312!------------------------------------------------------------------------------
2313! Compute fluxes in free convection
2314!------------------------------------------------------------------------------
2315
2316par_conv_visc = (t_s-t_a)/t_s*sqrt(tpsf_kappa_t_a) + (q_s-q_a)*tpsf_alpha_q*sqrt(tpsf_kappa_q_a)
2317IF(par_conv_visc.GT.0.0) THEN ! Viscous convection takes place
2318 l_conv_visc = .true.
2319 par_conv_visc = (par_conv_visc*tpl_grav/tpsf_nu_u_a)**num_1o3_sf
2320 q_sen_con = c_free_conv*sqrt(tpsf_kappa_t_a)*par_conv_visc
2321 q_sen_con = q_sen_con*(t_s-t_a)
2322 q_lat_con = c_free_conv*sqrt(tpsf_kappa_q_a)*par_conv_visc
2323 q_lat_con = q_lat_con*(q_s-q_a)
2324ELSE ! No viscous convection, set fluxes to zero
2325 l_conv_visc = .false.
2326 q_sen_con = 0.0
2327 q_lat_con = 0.0
2328END IF
2329q_mom_con = 0.0 ! Momentum flux in free (viscous or CBL-scale) convection is zero
2330
2331!------------------------------------------------------------------------------
2332! Compute turbulent fluxes
2333!------------------------------------------------------------------------------
2334
2335r_z = height_tq/height_u ! Ratio of "height_tq" to "height_u"
2336ri_cr = c_mo_t_stab/c_mo_u_stab**2_iintegers*r_z ! Critical Ri
2337ri = tpl_grav*((t_a-t_s)/t_s+tpsf_alpha_q*(q_a-q_s))/max(u_a,u_wind_min_sf)**2_iintegers
2338ri = ri*height_u/pr_neutral ! Gradient Richardson number
2339
2340turb_fluxes: IF(u_a.LT.u_wind_min_sf.OR.ri.GT.ri_cr-c_small_sf) THEN ! Low wind or Ri>Ri_cr
2341
2342u_star_st = 0.0 ! Set turbulent fluxes to zero
2343q_mom_tur = 0.0
2344q_sen_tur = 0.0
2345q_lat_tur = 0.0
2346
2347ELSE turb_fluxes ! Compute turbulent fluxes using MO similarity
2348
2349! Compute z/L, where z=height_u
2350IF(ri.GE.0.0) THEN ! Stable stratification
2351 zol = sqrt(1.0-4.0*(c_mo_u_stab-r_z*c_mo_t_stab)*ri)
2352 zol = zol - 1.0 + 2.0*c_mo_u_stab*ri
2353 zol = zol/2.0/c_mo_u_stab/c_mo_u_stab/(ri_cr-ri)
2354ELSE ! Convection
2355 n_iter = 0_iintegers
2356 delta = 1.0 ! Set initial error to a large value (as compared to the accuracy)
2357 u_star_previter = ri*max(1.0, sqrt(r_z*c_mo_t_conv/c_mo_u_conv)) ! Initial guess for ZoL
2358 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2359 fun = u_star_previter**2_iintegers*(c_mo_u_conv*u_star_previter-1.0) &
2360 + ri**2_iintegers*(1.0-r_z*c_mo_t_conv*u_star_previter)
2361 fun_prime = 3.0*c_mo_u_conv*u_star_previter**2_iintegers &
2362 - 2.0*u_star_previter - r_z*c_mo_t_conv*ri**2_iintegers
2363 zol = u_star_previter - fun/fun_prime
2364 delta = abs(zol-u_star_previter)/max(c_accur_sf, abs(zol+u_star_previter))
2365 u_star_previter = zol
2366 n_iter = n_iter + 1_iintegers
2367 END DO
2368!_dbg
2369! IF(n_iter.GE.n_iter_max-1_iintegers) &
2370! print*(*,*) 'ZoL: Max No. iters. exceeded (n_iter = ', n_iter, ')!'
2371!_dbg
2372END IF
2373
2374! Compute fetch-dependent Charnock parameter, use "u_star_min_sf"
2375CALL sfcflx_roughness (fetch, u_a, u_star_min_sf, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2376
2377! Threshold value of wind speed
2378u_star_st = u_star_thresh
2379CALL sfcflx_roughness (fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2380IF(zol.GT.0.0) THEN ! MO function in stable stratification
2381 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2382ELSE ! MO function in convection
2383 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2384 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2385 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2386 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2387 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2388END IF
2389u_a_thresh = u_star_thresh/c_karman*(log(height_u/z0u_sf)+psi_u)
2390
2391! Compute friction velocity
2392n_iter = 0_iintegers
2393delta = 1.0 ! Set initial error to a large value (as compared to the accuracy)
2394u_star_previter = u_star_thresh ! Initial guess for friction velocity
2395IF(u_a.LE.u_a_thresh) THEN ! Smooth surface
2396 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2397 CALL sfcflx_roughness (fetch, u_a, min(u_star_thresh, u_star_previter), h_ice, &
2398 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2399 IF(zol.GE.0.0) THEN ! Stable stratification
2400 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2401 fun = log(height_u/z0u_sf) + psi_u
2402 fun_prime = (fun + 1.0 + c_mo_u_stab*zol*min(z0u_sf/height_u, 1.0))/c_karman
2403 fun = fun*u_star_previter/c_karman - u_a
2404 ELSE ! Convection
2405 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2406 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2407 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2408 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2409 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2410 fun = log(height_u/z0u_sf) + psi_u
2411 fun_prime = (fun + 1.0/psi_q)/c_karman
2412 fun = fun*u_star_previter/c_karman - u_a
2413 END IF
2414 u_star_st = u_star_previter - fun/fun_prime
2415 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
2416 u_star_previter = u_star_st
2417 n_iter = n_iter + 1_iintegers
2418 END DO
2419ELSE ! Rough surface
2420 DO WHILE (delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
2421 CALL sfcflx_roughness (fetch, u_a, max(u_star_thresh, u_star_previter), h_ice, &
2422 c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2423 IF(zol.GE.0.0) THEN ! Stable stratification
2424 psi_u = c_mo_u_stab*zol*(1.0-min(z0u_sf/height_u, 1.0))
2425 fun = log(height_u/z0u_sf) + psi_u
2426 fun_prime = (fun - 2.0 - 2.0*c_mo_u_stab*zol*min(z0u_sf/height_u, 1.0))/c_karman
2427 fun = fun*u_star_previter/c_karman - u_a
2428 ELSE ! Convection
2429 psi_t = (1.0-c_mo_u_conv*zol)**c_mo_u_exp
2430 psi_q = (1.0-c_mo_u_conv*zol*min(z0u_sf/height_u, 1.0))**c_mo_u_exp
2431 psi_u = 2.0*(atan(psi_t)-atan(psi_q)) &
2432 + 2.0*log((1.0+psi_q)/(1.0+psi_t)) &
2433 + log((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
2434 fun = log(height_u/z0u_sf) + psi_u
2435 fun_prime = (fun - 2.0/psi_q)/c_karman
2436 fun = fun*u_star_previter/c_karman - u_a
2437 END IF
2438 IF(h_ice.GE.h_ice_min_flk) THEN ! No iteration is required for rough flow over ice
2439 u_star_st = c_karman*u_a/max(c_small_sf, log(height_u/z0u_sf)+psi_u)
2440 u_star_previter = u_star_st
2441 ELSE ! Iterate in case of open water
2442 u_star_st = u_star_previter - fun/fun_prime
2443 END IF
2444 delta = abs((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
2445 u_star_previter = u_star_st
2446 n_iter = n_iter + 1_iintegers
2447 END DO
2448END IF
2449
2450!_dbg
2451! print*(*,*) 'MO stab. func. psi_u = ', psi_u, ' n_iter = ', n_iter
2452! print*(*,*) ' Wind speed = ', U_a, ' u_* = ', u_star_st
2453! print*(*,*) ' Fun = ', Fun
2454!_dbg
2455
2456!_dbg
2457! IF(n_iter.GE.n_iter_max-1_iintegers) &
2458! print*(*,*) 'u_*: Max No. iters. exceeded (n_iter = ', n_iter, ')!'
2459!_dbg
2460
2461! Momentum flux
2462q_mom_tur = -u_star_st*u_star_st
2463
2464! Temperature and specific humidity fluxes
2465CALL sfcflx_roughness (fetch, u_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
2466IF(zol.GE.0.0) THEN ! Stable stratification
2467 psi_t = c_mo_t_stab*r_z*zol*(1.0-min(z0t_sf/height_tq, 1.0))
2468 psi_q = c_mo_q_stab*r_z*zol*(1.0-min(z0q_sf/height_tq, 1.0))
2469!_dbg
2470! print*(*,*) 'STAB: psi_t = ', psi_t, ' psi_q = ', psi_q
2471!_dbg
2472ELSE ! Convection
2473 psi_u = (1.0-c_mo_t_conv*r_z*zol)**c_mo_t_exp
2474 psi_t = (1.0-c_mo_t_conv*r_z*zol*min(z0t_sf/height_tq, 1.0))**c_mo_t_exp
2475! psi_t = 2.0*LOG((1.0+psi_t)/(1.0+psi_u))
2476 psi_t = abs(2.0*log((1.0+psi_t)/(1.0+psi_u)))
2477 psi_u = (1.0-c_mo_q_conv*r_z*zol)**c_mo_q_exp
2478 psi_q = (1.0-c_mo_q_conv*r_z*zol*min(z0q_sf/height_tq, 1.0))**c_mo_q_exp
2479! psi_q = 2.0*LOG((1.0+psi_q)/(1.0+psi_u))
2480 psi_q = abs(2.0*log((1.0+psi_q)/(1.0+psi_u)))
2481! write(0,*) 'psi_q= ',psi_q
2482!_dbg
2483! print*(*,*) 'CONV: psi_t = ', psi_t, ' psi_q = ', psi_q
2484!_dbg
2485END IF
2486q_sen_tur = -(t_a-t_s)*u_star_st*c_karman/pr_neutral &
2487 / max(c_small_sf, log(height_tq/z0t_sf)+psi_t)
2488if(max(c_small_sf, log(height_tq/z0t_sf)+psi_t) .lt. 10e-6) then
2489 write(0,*)'inside flake'
2490 write(0,*) q_sen_tur, t_a, t_s, u_star_st, c_karman, pr_neutral
2491 write(0,*) c_small_sf,height_tq,z0t_sf,psi_t
2492 write(0,*) 'nominator= ', (t_a-t_s)*u_star_st*c_karman/pr_neutral
2493 write(0,*) 'denominator= ',max(c_small_sf, log(height_tq/z0t_sf)+psi_t)
2494endif
2495q_lat_tur = -(q_a-q_s)*u_star_st*c_karman/sc_neutral &
2496 / max(c_small_sf, log(height_tq/z0q_sf)+psi_q)
2497if(q_lat_tur .gt. 6.0e-4) then
2498 q_lat_tur = -(q_a-q_s)*u_star_st*c_karman/3.0 &
2499 / max(c_small_sf, log(height_tq/z0q_sf)+psi_q)
2500 write(0,*) 'Q_lat_tur= ',q_lat_tur
2501 write(0,135) q_a,q_s,u_star_st,c_karman
2502 write(0,136) max(c_small_sf,log(height_tq/z0q_sf)+psi_q),c_small_sf, log(height_tq/z0q_sf),psi_q
2503endif
2504135 format(1x,4(f16.4))
2505136 format(1x,4(f16.4))
2506
2507END IF turb_fluxes
2508
2509!------------------------------------------------------------------------------
2510! Decide between turbulent, molecular, and convective fluxes
2511!------------------------------------------------------------------------------
2512
2513q_momentum = min(q_mom_tur, q_mom_mol, q_mom_con) ! Momentum flux is negative
2514IF(l_conv_visc) THEN ! Convection, take fluxes that are maximal in magnitude
2515 IF(abs(q_sen_tur).GE.abs(q_sen_con)) THEN
2516 q_sensible = q_sen_tur
2517 ELSE
2518 q_sensible = q_sen_con
2519 END IF
2520 IF(abs(q_sensible).LT.abs(q_sen_mol)) THEN
2521 q_sensible = q_sen_mol
2522 END IF
2523 IF(abs(q_lat_tur).GE.abs(q_lat_con)) THEN
2524 q_latent = q_lat_tur
2525 ELSE
2526 q_latent = q_lat_con
2527 END IF
2528 IF(abs(q_latent).LT.abs(q_lat_mol)) THEN
2529 q_latent = q_lat_mol
2530 END IF
2531ELSE ! Stable or neutral stratification, chose fluxes that are maximal in magnitude
2532 IF(abs(q_sen_tur).GE.abs(q_sen_mol)) THEN
2533 q_sensible = q_sen_tur
2534 ELSE
2535 q_sensible = q_sen_mol
2536 END IF
2537 IF(abs(q_lat_tur).GE.abs(q_lat_mol)) THEN
2538 q_latent = q_lat_tur
2539 ELSE
2540 q_latent = q_lat_mol
2541 END IF
2542END IF
2543
2544!------------------------------------------------------------------------------
2545! Set output (notice that fluxes are no longer in kinematic units)
2546!------------------------------------------------------------------------------
2547
2548q_momentum = q_momentum*rho_a
2549!Q_sensible = Q_sensible*rho_a*tpsf_c_a_p
2550!write(0,*) 'Q_sensible= ',Q_sensible
2551
2552q_watvap = q_latent*rho_a
2553
2554!Q_latent = tpsf_L_evap
2555IF(h_ice.GE.h_ice_min_flk) q_latent = q_latent + tpl_l_f ! Add latent heat of fusion over ice
2556!Q_latent = Q_watvap*Q_latent
2557q_latent = q_watvap*tpsf_l_evap
2558if(q_latent .gt. 2000.00) then
2559 write(0,145) 'final Q_watvap= ',q_watvap, 'tpsf_L_evap= ',tpsf_l_evap, 'Q_latent= ', q_latent
2560endif
2561!Q_latent = Q_watvap*Q_latent
2562145 format(a17,e12.5,1x,a13,1x,f10.2,1x,a10,1x,e12.4)
2563! Set "*_sf" variables to make fluxes accessible to driving routines that use "SfcFlx"
2564u_star_a_sf = u_star_st
2565q_mom_a_sf = q_momentum
2566q_sens_a_sf = q_sensible
2567q_lat_a_sf = q_latent
2568q_watvap_a_sf = q_watvap
2569
2570!write(85,127) Q_sensible, Q_watvap, Q_latent
2571 127 format(1x, 3(f16.5,1x))
2572
2573!------------------------------------------------------------------------------
2574! End calculations
2575!==============================================================================
2576
2577END SUBROUTINE sfcflx_momsenlat
2578
2579!==============================================================================
2580
2581!==============================================================================
2582! include 'SfcFlx_rhoair.incf'
2583!------------------------------------------------------------------------------
2584
2585REAL (KIND = kind_phys) FUNCTION sfcflx_rhoair (T, q, P)
2586
2587!------------------------------------------------------------------------------
2588!
2589! Description:
2590!
2591! Computes the air density as function
2592! of temperature, specific humidity and pressure.
2593!
2594! Declarations:
2595!
2596! Modules used:
2597
2598!_dm Parameters are USEd in module "SfcFlx".
2599!_nu USE data_parameters , ONLY : &
2600!_nu ireals, & ! KIND-type parameter for real variables
2601!_nu iintegers ! KIND-type parameter for "normal" integer variables
2602
2603use machine, only: kind_phys
2604!==============================================================================
2605
2606IMPLICIT NONE
2607
2608!==============================================================================
2609!
2610! Declarations
2611
2612! Input (function argument)
2613real(kind = kind_phys), INTENT(IN) :: &
2614 t , & ! Temperature [K]
2615 q , & ! Specific humidity
2616 p ! Pressure [N m^{-2} = kg m^{-1} s^{-2}]
2617
2618!==============================================================================
2619! Start calculations
2620!------------------------------------------------------------------------------
2621
2622! Air density [kg m^{-3}]
2623
2624sfcflx_rhoair = p/tpsf_r_dryair/t/(1.0+(1.0/tpsf_rd_o_rv-1.0)*q)
2625
2626!------------------------------------------------------------------------------
2627! End calculations
2628!==============================================================================
2629
2630END FUNCTION sfcflx_rhoair
2631
2632!==============================================================================
2633
2634!==============================================================================
2635! include 'SfcFlx_roughness.incf'
2636!------------------------------------------------------------------------------
2637
2638SUBROUTINE sfcflx_roughness (fetch, U_a, u_star, h_ice, &
2639 c_z0u_fetch, u_star_thresh, z0u, z0t, z0q)
2640
2641!------------------------------------------------------------------------------
2642!
2643! Description:
2644!
2645! Computes the water-surface or the ice-surface roughness lengths
2646! with respect to wind velocity, potential temperature and specific humidity.
2647!
2648! The water-surface roughness lengths with respect to wind velocity is computed
2649! from the Charnock formula when the surface is aerodynamically rough.
2650! A simple empirical formulation is used to account for the dependence
2651! of the Charnock parameter on the wind fetch.
2652! When the flow is aerodynamically smooth, the roughness length with respect to
2653! wind velocity is proportional to the depth of the viscous sub-layer.
2654! The water-surface roughness lengths for scalars are computed using the power-law
2655! formulations in terms of the roughness Reynolds number (Zilitinkevich et al. 2001).
2656! The ice-surface aerodynamic roughness is taken to be constant.
2657! The ice-surface roughness lengths for scalars
2658! are computed through the power-law formulations
2659! in terms of the roughness Reynolds number (Andreas 2002).
2660!
2661! Declarations:
2662!
2663! Modules used:
2664
2665!_dm Parameters are USEd in module "SfcFlx".
2666!_nu USE data_parameters , ONLY : &
2667!_nu ireals , & ! KIND-type parameter for real variables
2668!_nu iintegers ! KIND-type parameter for "normal" integer variables
2669
2670use machine, only: kind_phys
2671!==============================================================================
2672
2673IMPLICIT NONE
2674
2675!==============================================================================
2676!
2677! Declarations
2678
2679! Input (procedure arguments)
2680real(kind = kind_phys), INTENT(IN) :: &
2681 fetch , & ! Typical wind fetch [m]
2682 u_a , & ! Wind speed [m s^{-1}]
2683 u_star , & ! Friction velocity in the surface air layer [m s^{-1}]
2684 h_ice ! Ice thickness [m]
2685
2686! Output (procedure arguments)
2687real(kind = kind_phys), INTENT(OUT) :: &
2688 c_z0u_fetch , & ! Fetch-dependent Charnock parameter
2689 u_star_thresh , & ! Threshold value of friction velocity [m s^{-1}]
2690 z0u , & ! Roughness length with respect to wind velocity [m]
2691 z0t , & ! Roughness length with respect to potential temperature [m]
2692 z0q ! Roughness length with respect to specific humidity [m]
2693
2694! Local variables of type REAL
2695real(kind = kind_phys) :: &
2696 re_s , & ! Surface Reynolds number
2697 re_s_thresh ! Threshold value of Re_s
2698
2699!==============================================================================
2700! Start calculations
2701!------------------------------------------------------------------------------
2702
2703water_or_ice: IF(h_ice.LT.h_ice_min_flk) THEN ! Water surface
2704
2705! The Charnock parameter as dependent on dimensionless fetch
2706 c_z0u_fetch = max(u_a, u_wind_min_sf)**2_iintegers/tpl_grav/fetch ! Inverse dimensionless fetch
2707 c_z0u_fetch = c_z0u_rough + c_z0u_ftch_f*c_z0u_fetch**c_z0u_ftch_ex
2708 c_z0u_fetch = min(c_z0u_fetch, c_z0u_rough_l) ! Limit Charnock parameter
2709
2710! Threshold value of friction velocity
2711 u_star_thresh = (c_z0u_smooth/c_z0u_fetch*tpl_grav*tpsf_nu_u_a)**num_1o3_sf
2712
2713! Surface Reynolds number and its threshold value
2714 re_s = u_star**3_iintegers/tpsf_nu_u_a/tpl_grav
2715 re_s_thresh = c_z0u_smooth/c_z0u_fetch
2716
2717! Aerodynamic roughness
2718 IF(re_s.LE.re_s_thresh) THEN
2719 z0u = c_z0u_smooth*tpsf_nu_u_a/u_star ! Smooth flow
2720 ELSE
2721 z0u = c_z0u_fetch*u_star*u_star/tpl_grav ! Rough flow
2722 END IF
2723! Roughness for scalars
2724 z0q = c_z0u_fetch*max(re_s, re_s_thresh)
2725 z0t = c_z0t_rough_1*z0q**c_z0t_rough_3 - c_z0t_rough_2
2726 z0q = c_z0q_rough_1*z0q**c_z0q_rough_3 - c_z0q_rough_2
2727 z0t = z0u*exp(-c_karman/pr_neutral*z0t)
2728 z0q = z0u*exp(-c_karman/sc_neutral*z0q)
2729
2730ELSE water_or_ice ! Ice surface
2731
2732! The Charnock parameter is not used over ice, formally set "c_z0u_fetch" to its minimum value
2733 c_z0u_fetch = c_z0u_rough
2734
2735! Threshold value of friction velocity
2736 u_star_thresh = c_z0u_smooth*tpsf_nu_u_a/z0u_ice_rough
2737
2738! Aerodynamic roughness
2739 z0u = max(z0u_ice_rough, c_z0u_smooth*tpsf_nu_u_a/u_star)
2740
2741! Roughness Reynolds number
2742 re_s = max(u_star*z0u/tpsf_nu_u_a, c_accur_sf)
2743
2744! Roughness for scalars
2745 IF(re_s.LE.re_z0s_ice_t) THEN
2746 z0t = c_z0t_ice_b0t + c_z0t_ice_b1t*log(re_s)
2747 z0t = min(z0t, c_z0t_ice_b0s)
2748 z0q = c_z0q_ice_b0t + c_z0q_ice_b1t*log(re_s)
2749 z0q = min(z0q, c_z0q_ice_b0s)
2750 ELSE
2751 z0t = c_z0t_ice_b0r + c_z0t_ice_b1r*log(re_s) + c_z0t_ice_b2r*log(re_s)**2_iintegers
2752 z0q = c_z0q_ice_b0r + c_z0q_ice_b1r*log(re_s) + c_z0q_ice_b2r*log(re_s)**2_iintegers
2753 END IF
2754 z0t = z0u*exp(z0t)
2755 z0q = z0u*exp(z0q)
2756
2757END IF water_or_ice
2758
2759!------------------------------------------------------------------------------
2760! End calculations
2761!==============================================================================
2762
2763END SUBROUTINE sfcflx_roughness
2764
2765!==============================================================================
2766
2767!==============================================================================
2768! include 'SfcFlx_satwvpres.incf'
2769!------------------------------------------------------------------------------
2770
2771REAL (KIND = kind_phys) FUNCTION sfcflx_satwvpres (T, h_ice)
2772
2773!------------------------------------------------------------------------------
2774!
2775! Description:
2776!
2777! Computes saturation water vapour pressure
2778! over the water surface or over the ice surface
2779! as function of temperature.
2780!
2781! Declarations:
2782!
2783! Modules used:
2784
2785!_dm Parameters are USEd in module "SfcFlx".
2786!_nu USE data_parameters , ONLY : &
2787!_nu ireals, & ! KIND-type parameter for real variables
2788!_nu iintegers ! KIND-type parameter for "normal" integer variables
2789
2790!_dm The variable is USEd in module "SfcFlx".
2791!_nu USE flake_parameters , ONLY : &
2792!_nu h_Ice_min_flk ! Minimum ice thickness [m]
2793use machine, only: kind_phys
2794
2795!==============================================================================
2796
2797IMPLICIT NONE
2798
2799!==============================================================================
2800!
2801! Declarations
2802
2803! Input (function argument)
2804real(kind = kind_phys), INTENT(IN) :: &
2805 t , & ! Temperature [K]
2806 h_ice ! Ice thickness [m]
2807
2808! Local parameters
2809real(kind = kind_phys), PARAMETER :: &
2810 b1_vap = 610.780 , & ! Coefficient [N m^{-2} = kg m^{-1} s^{-2}]
2811 b3_vap = 273.160 , & ! Triple point [K]
2812 b2w_vap = 17.26938820 , & ! Coefficient (water)
2813 b2i_vap = 21.87455840 , & ! Coefficient (ice)
2814 b4w_vap = 35.860 , & ! Coefficient (temperature) [K]
2815 b4i_vap = 7.660 ! Coefficient (temperature) [K]
2816
2817!==============================================================================
2818! Start calculations
2819!------------------------------------------------------------------------------
2820
2821! Saturation water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
2822
2823IF(h_ice.LT.h_ice_min_flk) THEN ! Water surface
2824 sfcflx_satwvpres = b1_vap*exp(b2w_vap*(t-b3_vap)/(t-b4w_vap))
2825ELSE ! Ice surface
2826 sfcflx_satwvpres = b1_vap*exp(b2i_vap*(t-b3_vap)/(t-b4i_vap))
2827END IF
2828
2829!------------------------------------------------------------------------------
2830! End calculations
2831!==============================================================================
2832
2833END FUNCTION sfcflx_satwvpres
2834
2835!==============================================================================
2836
2837!==============================================================================
2838! include 'SfcFlx_spechum.incf'
2839!------------------------------------------------------------------------------
2840
2841REAL (KIND = kind_phys) FUNCTION sfcflx_spechum (wvpres, P)
2842
2843!------------------------------------------------------------------------------
2844!
2845! Description:
2846!
2847! Computes specific humidity as function
2848! of water vapour pressure and air pressure.
2849!
2850! Declarations:
2851!
2852! Modules used:
2853
2854!_dm Parameters are USEd in module "SfcFlx".
2855!_nu USE data_parameters , ONLY : &
2856!_nu ireals, & ! KIND-type parameter for real variables
2857!_nu iintegers ! KIND-type parameter for "normal" integer variables
2858
2859use machine, only: kind_phys
2860!==============================================================================
2861
2862IMPLICIT NONE
2863
2864!==============================================================================
2865!
2866! Declarations
2867
2868! Input (function argument)
2869real(kind = kind_phys), INTENT(IN) :: &
2870 wvpres , & ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
2871 p ! Air pressure [N m^{-2} = kg m^{-1} s^{-2}]
2872
2873!==============================================================================
2874! Start calculations
2875!------------------------------------------------------------------------------
2876
2877! Specific humidity
2878
2879sfcflx_spechum = tpsf_rd_o_rv*wvpres/(p-(1.0-tpsf_rd_o_rv)*wvpres)
2880
2881!------------------------------------------------------------------------------
2882! End calculations
2883!==============================================================================
2884
2885END FUNCTION sfcflx_spechum
2886
2887!==============================================================================
2888
2889!==============================================================================
2890! include 'SfcFlx_wvpreswetbulb.incf'
2891!------------------------------------------------------------------------------
2892
2893REAL (KIND = ireals) FUNCTION sfcflx_wvpreswetbulb (T_dry, T_wetbulb, satwvpres_bulb, P)
2894
2895!------------------------------------------------------------------------------
2896!
2897! Description:
2898!
2899! Computes water vapour pressure as function of air temperature,
2900! wet bulb temperature, satururation vapour pressure at wet-bulb temperature,
2901! and air pressure.
2902!
2903! Declarations:
2904!
2905! Modules used:
2906
2907!_dm Parameters are USEd in module "SfcFlx".
2908!_nu USE data_parameters , ONLY : &
2909!_nu ireals, & ! KIND-type parameter for real variables
2910!_nu iintegers ! KIND-type parameter for "normal" integer variables
2911
2912use machine, only: kind_phys
2913!==============================================================================
2914
2915IMPLICIT NONE
2916
2917!==============================================================================
2918!
2919! Declarations
2920
2921! Input (function argument)
2922real(kind = kind_phys), INTENT(IN) :: &
2923 t_dry , & ! Dry air temperature [K]
2924 t_wetbulb , & ! Wet bulb temperature [K]
2925 satwvpres_bulb , & ! Satururation vapour pressure at wet-bulb temperature [N m^{-2}]
2926 p ! Atmospheric pressure [N m^{-2}]
2927
2928!==============================================================================
2929! Start calculations
2930!------------------------------------------------------------------------------
2931
2932! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
2933
2934sfcflx_wvpreswetbulb = satwvpres_bulb &
2935 - tpsf_c_a_p*p/tpsf_l_evap/tpsf_rd_o_rv*(t_dry-t_wetbulb)
2936
2937
2938!------------------------------------------------------------------------------
2939! End calculations
2940!==============================================================================
2941
2942END FUNCTION sfcflx_wvpreswetbulb
2943
2944!==============================================================================
2945
2946END MODULE sfcflx
2947
2948
2950IMPLICIT NONE
2951CONTAINS
2952
2953!------------------------------------------------------------------------------
2954
2955SUBROUTINE flake_interface ( dMsnowdt_in, I_atm_in, Q_atm_lw_in, height_u_in, height_tq_in, &
2956 U_a_in, T_a_in, q_a_in, P_a_in, &
2957
2958 depth_w, fetch, depth_bs, T_bs, par_Coriolis, del_time, &
2959 T_snow_in, T_ice_in, T_mnw_in, T_wML_in, T_bot_in, T_B1_in, &
2960 C_T_in, h_snow_in, h_ice_in, h_ML_in, H_B1_in, T_sfc_p, &
2961 ch, cm, albedo_water, water_extinc, &
2962
2963 T_snow_out, T_ice_out, T_mnw_out, T_wML_out, T_bot_out, &
2964 T_B1_out, C_T_out, h_snow_out, h_ice_out, h_ML_out, &
2965 H_B1_out, T_sfc_n, hflx_out, evap_out, gflx_out, lflx_out, &
2966
2967 T_bot_2_in, T_bot_2_out,ustar, q_sfc, chh, cmm )
2968
2969!------------------------------------------------------------------------------
2970!
2971! Description:
2972!
2973! The FLake interface is
2974! a communication routine between "flake_main"
2975! and a prediction system that uses FLake.
2976! It assigns the FLake variables at the previous time step
2977! to their input values given by the driving model,
2978! calls a number of routines to compute the heat and radiation fluxes,
2979! calls "flake_main",
2980! and returns the updated FLake variables to the driving model.
2981! The "flake_interface" does not contain any Flake physics.
2982! It only serves as a convenient means to organize calls of "flake_main"
2983! and of external routines that compute heat and radiation fluxes.
2984! The interface may (should) be changed so that to provide
2985! the most convenient use of FLake.
2986! Within a 3D atmospheric prediction system,
2987! "flake_main" may be called in a DO loop within "flake_interface"
2988! for each grid-point where a lake is present.
2989! In this way, the driving atmospheric model should call "flake_interface"
2990! only once, passing the FLake variables to "flake_interface" as 2D fields.
2991!
2992! Lines embraced with "!_tmp" contain temporary parts of the code.
2993! These should be removed prior to using FLake in applications.
2994! Lines embraced/marked with "!_dev" may be replaced
2995! as improved parameterizations are developed and tested.
2996! Lines embraced/marked with "!_dm" are DM's comments
2997! that may be helpful to a user.
2998!
2999use machine, only: kind_phys
3000
3001USE data_parameters , ONLY : &
3002 ireals, & ! KIND-type parameter for real variables
3003 iintegers ! KIND-type parameter for "normal" integer variables
3004
3005USE flake_derivedtypes ! Definitions of several derived TYPEs
3006
3007!USE flake_parameters , ONLY : &
3008! tpl_T_f , & ! Fresh water freezing point [K]
3009! tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
3010! h_Snow_min_flk , & ! Minimum snow thickness [m]
3011! h_Ice_min_flk ! Minimum ice thickness [m]
3012
3013USE flake_paramoptic_ref ! Reference values of the optical characteristics
3014 ! of the lake water, lake ice and snow
3015
3016USE flake_albedo_ref ! Reference values the albedo for the lake water, lake ice and snow
3017
3018USE flake , ONLY : &
3019 flake_main , & ! Subroutine, FLake driver
3020 flake_radflux , & ! Subroutine, computes radiation fluxes at various depths
3021 !
3022 t_snow_p_flk, t_snow_n_flk , & ! Temperature at the air-snow interface [K]
3023 t_ice_p_flk, t_ice_n_flk , & ! Temperature at the snow-ice or air-ice interface [K]
3024 t_mnw_p_flk, t_mnw_n_flk , & ! Mean temperature of the water column [K]
3025 t_wml_p_flk, t_wml_n_flk , & ! Mixed-layer temperature [K]
3026 t_bot_p_flk, t_bot_n_flk , & ! Temperature at the water-bottom sediment interface [K]
3027 t_b1_p_flk, t_b1_n_flk , & ! Temperature at the bottom of the upper layer of the sediments [K]
3028 c_t_p_flk, c_t_n_flk , & ! Shape factor (thermocline)
3029 h_snow_p_flk, h_snow_n_flk , & ! Snow thickness [m]
3030 h_ice_p_flk, h_ice_n_flk , & ! Ice thickness [m]
3031 h_ml_p_flk, h_ml_n_flk , & ! Thickness of the mixed-layer [m]
3032 h_b1_p_flk, h_b1_n_flk , & ! Thickness of the upper layer of bottom sediments [m]
3033 !
3034 q_snow_flk , & ! Heat flux through the air-snow interface [W m^{-2}]
3035 q_ice_flk , & ! Heat flux through the snow-ice or air-ice interface [W m^{-2}]
3036 q_w_flk , & ! Heat flux through the ice-water or air-water interface [W m^{-2}]
3037 q_bot_flk , & ! Heat flux through the water-bottom sediment interface [W m^{-2}]
3038 i_atm_flk , & ! Radiation flux at the lower boundary of the atmosphere [W m^{-2}],
3039 ! i.e. the incident radiation flux with no regard for the surface albedo
3040 i_snow_flk , & ! Radiation flux through the air-snow interface [W m^{-2}]
3041 i_ice_flk , & ! Radiation flux through the snow-ice or air-ice interface [W m^{-2}]
3042 i_w_flk , & ! Radiation flux through the ice-water or air-water interface [W m^{-2}]
3043 i_h_flk , & ! Radiation flux through the mixed-layer-thermocline interface [W m^{-2}]
3044 i_bot_flk , & ! Radiation flux through the water-bottom sediment interface [W m^{-2}]
3045 i_intm_0_h_flk , & ! Mean radiation flux over the mixed layer [W m^{-1}]
3046 i_intm_h_d_flk , & ! Mean radiation flux over the thermocline [W m^{-1}]
3047 q_star_flk , & ! A generalized heat flux scale [W m^{-2}]
3048 u_star_w_flk , & ! Friction velocity in the surface layer of lake water [m s^{-1}]
3049 w_star_sfc_flk , & ! Convective velocity scale, using a generalized heat flux scale [m s^{-1}]
3050 dmsnowdt_flk , & ! The rate of snow accumulation [kg m^{-2} s^{-1}]
3051 t_bot_2_in_flk
3052
3053
3054USE sfcflx , ONLY : &
3055 sfcflx_lwradwsfc , & ! Function, returns the surface long-wave radiation flux
3056 sfcflx_momsenlat ! Subroutine, computes fluxes of momentum and of sensible and latent heat
3057
3058!==============================================================================
3059
3060IMPLICIT NONE
3061
3062!==============================================================================
3063!
3064! Declarations
3065
3066! Input (procedure arguments)
3067
3068real(kind = kind_phys), INTENT(IN) :: &
3069 dmsnowdt_in , & ! The rate of snow accumulation [kg m^{-2} s^{-1}]
3070 i_atm_in , & ! Solar radiation flux at the surface [W m^{-2}]
3071 q_atm_lw_in , & ! Long-wave radiation flux from the atmosphere [W m^{-2}]
3072 height_u_in , & ! Height above the lake surface where the wind speed is measured [m]
3073 height_tq_in , & ! Height where temperature and humidity are measured [m]
3074 u_a_in , & ! Wind speed at z=height_u_in [m s^{-1}]
3075 t_a_in , & ! Air temperature at z=height_tq_in [K]
3076 q_a_in , & ! Air specific humidity at z=height_tq_in
3077 p_a_in , & ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}]
3078 ch , &
3079 cm , &
3080 albedo_water, & ! Water surface albedo with respect to the solar radiation
3081 water_extinc
3082
3083real(kind = kind_phys), INTENT(IN) :: &
3084 depth_w , & ! The lake depth [m]
3085 fetch , & ! Typical wind fetch [m]
3086 depth_bs , & ! Depth of the thermally active layer of the bottom sediments [m]
3087 t_bs , & ! Temperature at the outer edge of
3088 ! the thermally active layer of the bottom sediments [K]
3089 par_coriolis , & ! The Coriolis parameter [s^{-1}]
3090 del_time ! The model time step [s]
3091
3092real(kind = kind_phys), INTENT(IN) :: &
3093 t_snow_in , & ! Temperature at the air-snow interface [K]
3094 t_ice_in , & ! Temperature at the snow-ice or air-ice interface [K]
3095 t_mnw_in , & ! Mean temperature of the water column [K]
3096 t_wml_in , & ! Mixed-layer temperature [K]
3097 t_bot_in , & ! Temperature at the water-bottom sediment interface [K]
3098 t_b1_in , & ! Temperature at the bottom of the upper layer of the sediments [K]
3099 c_t_in , & ! Shape factor (thermocline)
3100 h_snow_in , & ! Snow thickness [m]
3101 h_ice_in , & ! Ice thickness [m]
3102 h_ml_in , & ! Thickness of the mixed-layer [m]
3103 h_b1_in , & ! Thickness of the upper layer of bottom sediments [m]
3104 t_sfc_p , & ! Surface temperature at the previous time step [K]
3105 t_bot_2_in
3106
3107! Input/Output (procedure arguments)
3108
3109!REAL (KIND = ireals), INTENT(INOUT) :: &
3110real(kind = kind_phys) :: &
3111 albedo_ice , & ! Ice surface albedo with respect to the solar radiation
3112 albedo_snow ! Snow surface albedo with respect to the solar radiation
3113
3114!TYPE (opticpar_medium), INTENT(INOUT) :: &
3115TYPE (opticpar_medium) :: &
3116 opticpar_water , & ! Optical characteristics of water
3117 opticpar_ice , & ! Optical characteristics of ice
3118 opticpar_snow ! Optical characteristics of snow
3119
3120! Output (procedure arguments)
3121
3122real(kind = kind_phys), INTENT(OUT) :: &
3123 t_snow_out , & ! Temperature at the air-snow interface [K]
3124 t_ice_out , & ! Temperature at the snow-ice or air-ice interface [K]
3125 t_mnw_out , & ! Mean temperature of the water column [K]
3126 t_wml_out , & ! Mixed-layer temperature [K]
3127 t_bot_out , & ! Temperature at the water-bottom sediment interface [K]
3128 t_b1_out , & ! Temperature at the bottom of the upper layer of the sediments [K]
3129 c_t_out , & ! Shape factor (thermocline)
3130 h_snow_out , & ! Snow thickness [m]
3131 h_ice_out , & ! Ice thickness [m]
3132 h_ml_out , & ! Thickness of the mixed-layer [m]
3133 h_b1_out , & ! Thickness of the upper layer of bottom sediments [m]
3134 t_sfc_n , & ! Updated surface temperature [K]
3135 hflx_out , & ! sensibl heat flux
3136 evap_out , & ! Latent heat flux
3137 gflx_out , & ! flux from to water
3138 lflx_out , & ! latent heat flux
3139 t_bot_2_out , & ! Bottom temperature
3140 ustar , &
3141 q_sfc , &
3142 chh , &
3143 cmm
3144
3145! Local variables of type REAL
3146
3147real(kind = kind_phys) :: &
3148 q_momentum , & ! Momentum flux [N m^{-2}]
3149 q_sensible , & ! Sensible heat flux [W m^{-2}]
3150 q_latent , & ! Latent heat flux [W m^{-2}]
3151 q_watvap , & ! Flux of water vapour [kg m^{-2} s^{-1}]
3152 q_w_flux , & ! flux from ice to water
3153 rho_a
3154
3155! ADDED by Shaobo Zhang
3156LOGICAL lflk_botsed_use
3157!REAL (KIND = kind_phys) :: T_bot_2_in, T_bot_2_out
3158real(kind = kind_phys), parameter :: tpl_rho_w_r = 1.0e+03
3159real(kind = kind_phys), parameter :: tpl_t_f = 273.15
3160real(kind = kind_phys), parameter :: h_snow_min_flk = 1.0e-5
3161real(kind = kind_phys), parameter :: h_ice_min_flk = 1.0e-9
3162!==============================================================================
3163! Start calculations
3164!------------------------------------------------------------------------------
3165! lflk_botsed_use = .TRUE.
3166 lflk_botsed_use = .false.
3167!------------------------------------------------------------------------------
3168! Set albedos of the lake water, lake ice and snow
3169!------------------------------------------------------------------------------
3170
3171! Use default value
3172! albedo_water = albedo_water_ref
3173! Use empirical formulation proposed by Mironov and Ritter (2004) for GME
3174!_nu albedo_ice = albedo_whiteice_ref
3175!albedo_ice = EXP(-c_albice_MR*(tpl_T_f-T_sfc_p)/tpl_T_f)
3176!albedo_ice = albedo_whiteice_ref*(1.0-albedo_ice) + albedo_blueice_ref*albedo_ice
3177! Snow is not considered
3178!albedo_snow = albedo_ice
3179albedo_ice = albedo_whiteice_ref
3180!albedo_snow = albedo_ice
3181albedo_snow = albedo_drysnow_ref
3182opticpar_water%extincoef_optic(1) = water_extinc
3183!write(0,*)'albedo= ',albedo_water,albedo_ice,albedo_snow
3184
3185!------------------------------------------------------------------------------
3186! Set optical characteristics of the lake water, lake ice and snow
3187!------------------------------------------------------------------------------
3188
3189! Use default values
3190opticpar_water = opticpar_water_ref
3191opticpar_ice = opticpar_ice_opaque ! Opaque ice
3192opticpar_snow = opticpar_snow_opaque ! Opaque snow
3193
3194!print*,'opticpar = ',opticpar_water, opticpar_ice,opticpar_snow
3195
3196!------------------------------------------------------------------------------
3197! Set initial values
3198!------------------------------------------------------------------------------
3199!print*,'Inter depth_w=',depth_w
3200!print*,'Inter depth_bs=',depth_bs
3201
3202t_snow_p_flk = t_snow_in
3203t_ice_p_flk = t_ice_in
3204t_mnw_p_flk = t_mnw_in
3205t_wml_p_flk = t_wml_in
3206t_bot_p_flk = t_bot_in
3207t_b1_p_flk = t_b1_in
3208c_t_p_flk = c_t_in
3209h_snow_p_flk = h_snow_in
3210h_ice_p_flk = h_ice_in
3211h_ml_p_flk = h_ml_in
3212h_b1_p_flk = h_b1_in
3213t_bot_2_in_flk = t_bot_2_in
3214
3215!write(71,120) T_sfc_p,T_mnw_in,T_wML_in,T_bot_in,T_B1_in,T_bot_2_in
3216 120 format(1x,6(f12.5,1x))
3217!------------------------------------------------------------------------------
3218! Set the rate of snow accumulation
3219!------------------------------------------------------------------------------
3220
3221dmsnowdt_flk = dmsnowdt_in
3222
3223!------------------------------------------------------------------------------
3224! Compute solar radiation fluxes (positive downward)
3225!------------------------------------------------------------------------------
3226
3227i_atm_flk = i_atm_in
3228CALL flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
3229 opticpar_water, opticpar_ice, opticpar_snow, &
3230 depth_bs )
3231
3232!------------------------------------------------------------------------------
3233! Compute long-wave radiation fluxes (positive downward)
3234!------------------------------------------------------------------------------
3235
3236q_w_flk = q_atm_lw_in ! Radiation of the atmosphere
3237q_w_flk = q_w_flk - sfcflx_lwradwsfc(t_sfc_p) ! Radiation of the surface (notice the sign)
3238
3239!------------------------------------------------------------------------------
3240! Compute the surface friction velocity and fluxes of sensible and latent heat
3241!------------------------------------------------------------------------------
3242
3243CALL sfcflx_momsenlat ( height_u_in, height_tq_in, fetch, &
3244 u_a_in, t_a_in, q_a_in, t_sfc_p, p_a_in, h_ice_p_flk, &
3245 q_momentum, q_sensible, q_latent, q_watvap, q_sfc, rho_a )
3246!write(0,*)'tpl_rho_w_r= ', tpl_rho_w_r
3247!write(0,*) 'Q_momentum= ',Q_momentum
3248u_star_w_flk = sqrt(-q_momentum/tpl_rho_w_r)
3249ustar = u_star_w_flk
3250
3251!------------------------------------------------------------------------------
3252! Compute heat fluxes Q_snow_flk, Q_ice_flk, Q_w_flk
3253!------------------------------------------------------------------------------
3254
3255q_w_flk = q_w_flk - q_sensible - q_latent ! Add sensible and latent heat fluxes (notice the signs)
3256IF(h_ice_p_flk.GE.h_ice_min_flk) THEN ! Ice exists
3257 IF(h_snow_p_flk.GE.h_snow_min_flk) THEN ! There is snow above the ice
3258 q_snow_flk = q_w_flk
3259 q_ice_flk = 0.0
3260 q_w_flk = 0.0
3261 ELSE ! No snow above the ice
3262 q_snow_flk = 0.0
3263 q_ice_flk = q_w_flk
3264 q_w_flk = 0.0
3265 END IF
3266ELSE ! No ice-snow cover
3267 q_snow_flk = 0.0
3268 q_ice_flk = 0.0
3269END IF
3270
3271!------------------------------------------------------------------------------
3272! Advance FLake variables
3273!------------------------------------------------------------------------------
3274
3275CALL flake_main ( depth_w, depth_bs, t_bs, par_coriolis, &
3276 opticpar_water%extincoef_optic(1), &
3277 del_time, t_sfc_p, t_sfc_n, t_bot_2_in_flk, &
3278 t_bot_2_out )
3279
3280!------------------------------------------------------------------------------
3281! Set output values
3282!------------------------------------------------------------------------------
3283
3284t_snow_out = t_snow_n_flk
3285t_ice_out = t_ice_n_flk
3286t_mnw_out = t_mnw_n_flk
3287t_wml_out = t_wml_n_flk
3288t_bot_out = t_bot_n_flk
3289t_b1_out = t_b1_n_flk
3290c_t_out = c_t_n_flk
3291h_snow_out = h_snow_n_flk
3292h_ice_out = h_ice_n_flk
3293h_ml_out = h_ml_n_flk
3294h_b1_out = h_b1_n_flk
3295hflx_out = q_sensible
3296evap_out = q_watvap
3297!evap_out = Q_latent
3298gflx_out = q_w_flk
3299lflx_out = q_latent
3300chh = ch * u_a_in * rho_a
3301cmm = cm * u_a_in
3302
3303!write(72,120) T_sfc_n,T_mnw_out,T_wML_out,T_bot_out,T_B1_out,T_bot_2_out
3304!------------------------------------------------------------------------------
3305! End calculations
3306!==============================================================================
3307
3308END SUBROUTINE flake_interface
3309
3310END MODULE module_flake