!=======================================================================
! Current Code Owner: DWD, Ulrich Schaettler
! phone: +49 69 8062 2739
! fax: +49 69 8236 1493
! email: uschaettler@dwd.d400.de
!
! History:
! Version Date Name
! ---------- ---------- ----
! 1.1 1998/03/11 Ulrich Schaettler
! Initial release
! 1.8 1998/08/03 Ulrich Schaettler
! Eliminated intgribf, intgribc, irealgrib, iwlength and put it to data_io.
! 1.10 1998/09/29 Ulrich Schaettler
! Eliminated parameters for grid point and diagnostic calculations.
! !VERSION! !DATE! <Your name>
! <Modification comments>
!
! Code Description:
! Language: Fortran 90.
! Software Standards: "European Standards for Writing and
! Documenting Exchangeable Fortran 90 Code".
!
! reorganize the FLake to module_FLake.F90 by Shaobo Zhang in 2016-7-13
! added a new layer for deep lakes by Shaobo Zhang in 2016-11-15
!
!=======================================================================
!------------------------------------------------------------------------------
MODULE data_parameters
!------------------------------------------------------------------------------
!
! Description:
! Global parameters for the program are defined.
! Actually, scratch that. We'll import them from machine.F instead.
!
use machine, only: ireals=>kind_phys, iintegers=>kind_INTEGER
!=======================================================================
END MODULE data_parameters
!------------------------------------------------------------------------------
MODULE flake_albedo_ref
!------------------------------------------------------------------------------
!
! Description:
!
! This module contains "reference" values of albedo
! for the lake water, lake ice and snow.
! As in "flake_paramoptic_ref", two ice categories, viz. white ice and blue ice,
! and two snow categories, viz. dry snow and melting snow, are used.
!
USE data_parameters, ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Albedo for water, ice and snow.
!REAL (KIND = ireals), PARAMETER :: &
! albedo_water_ref = 0.070 , & ! Water
! albedo_whiteice_ref = 0.600 , & ! White ice
! albedo_blueice_ref = 0.100 , & ! Blue ice
! albedo_drysnow_ref = 0.600 , & ! Dry snow
! albedo_meltingsnow_ref = 0.100 ! Melting snow
! Empirical parameters.
!REAL (KIND = ireals), PARAMETER :: &
! c_albice_MR = 95.60 ! Constant in the interpolation formula for
! the ice albedo (Mironov and Ritter 2004)
! Albedo for water, ice and snow.
REAL (KIND = kind_phys), PARAMETER :: &
albedo_water_ref = 0.07 , & ! Water
albedo_whiteice_ref = 0.60 , & ! White ice
albedo_blueice_ref = 0.10 , & ! Blue ice
! albedo_drysnow_ref = 0.60 , & ! Dry snow
albedo_drysnow_ref = 0.90 , & ! Dry snow
albedo_meltingsnow_ref = 0.10 ! Melting snow
! Empirical parameters.
REAL (KIND = kind_phys), PARAMETER :: &
c_albice_MR = 95.6 ! Constant in the interpolation formula for
! the ice albedo (Mironov and Ritter 2004)
!==============================================================================
END MODULE flake_albedo_ref
!------------------------------------------------------------------------------
MODULE flake_configure
!------------------------------------------------------------------------------
!
! Description:
!
! Switches and reference values of parameters
! that configure the lake model FLake are set.
!
USE data_parameters , ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! changed by Shaobo Zhang
LOGICAL lflk_botsed_use
!LOGICAL, PARAMETER :: &
! lflk_botsed_use = .TRUE. ! .TRUE. indicates that the bottom-sediment scheme is used
! to compute the depth penetrated by the thermal wave,
! the temperature at this depth and the bottom heat flux.
! Otherwise, the heat flux at the water-bottom sediment interface
! is set to zero, the depth penetrated by the thermal wave
! is set to a reference value defined below,
! and the temperature at this depth is set to
! the temperature of maximum density of the fresh water.
!REAL (KIND = ireals), PARAMETER :: &
! rflk_depth_bs_ref = 10.00 ! Reference value of the depth of the thermally active
! layer of bottom sediments [m].
! This value is used to (formally) define
! the depth penetrated by the thermal wave
! in case the bottom-sediment scheme is not used.
REAL (KIND = kind_phys), PARAMETER :: &
rflk_depth_bs_ref = 10.0
!==============================================================================
END MODULE flake_configure
!------------------------------------------------------------------------------
MODULE flake_derivedtypes
!------------------------------------------------------------------------------
!
! Description:
!
! Derived type(s) is(are) defined.
!
USE data_parameters , ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Maximum value of the wave-length bands
! in the exponential decay law for the radiation flux.
! A storage for a ten-band approximation is allocated,
! although a smaller number of bands is actually used.
INTEGER (KIND = iintegers), PARAMETER :: &
nband_optic_max = 10_iintegers
! Define TYPE "opticpar_medium"
TYPE opticpar_medium
INTEGER (KIND = iintegers) :: &
nband_optic ! Number of wave-length bands
REAL (KIND = ireals), DIMENSION (nband_optic_max) :: &
frac_optic , & ! Fractions of total radiation flux
extincoef_optic ! Extinction coefficients
END TYPE opticpar_medium
!==============================================================================
END MODULE flake_derivedtypes
!------------------------------------------------------------------------------
MODULE flake_paramoptic_ref
!------------------------------------------------------------------------------
!
! Description:
!
! This module contains "reference" values of the optical characteristics
! of the lake water, lake ice and snow. These reference values may be used
! if no information about the optical characteristics of the lake in question
! is available. An exponential decay law for the solar radiation flux is assumed.
! In the simplest one-band approximation,
! the extinction coefficient for water is set to a large value,
! leading to the absorption of 95% of the incoming radiation
! within the uppermost 1 m of the lake water.
! The extinction coefficients for ice and snow are taken from
! Launiainen and Cheng (1998). The estimates for the ice correspond
! to the uppermost 0.1 m of the ice layer and to the clear sky conditions
! (see Table 2 in op. cit.).
! Very large values of the extinction coefficients for ice and snow ("opaque")
! can be used to prevent penetration of the solar radiation
! through the snow-ice cover.
!
USE data_parameters, ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
USE flake_derivedtypes, ONLY : &
nband_optic_max , & ! Maximum value of the wave-length bands
opticpar_medium ! Derived TYPE
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
INTEGER (KIND = iintegers), PRIVATE :: & ! Help variable(s)
i ! DO loop index
! Optical characteristics for water, ice and snow.
! The simplest one-band approximation is used as a reference.
!TYPE (opticpar_medium), PARAMETER :: &
! opticpar_water_ref = opticpar_medium(1, & ! Water (reference)
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/3.0, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_water_trans = opticpar_medium(2, & ! Transparent Water (two-band)
! (/0.100, 0.900, (0.0,i=3,nband_optic_max)/), &
! (/2.00, 0.200, (1.E+100,i=3,nband_optic_max)/)) , &
!!_nu opticpar_water_trans = opticpar_medium(1, & ! Transparent Water (one-band)
!!_nu (/1.0, (0.0,i=2,nband_optic_max)/), &
!!_nu (/0.300, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_whiteice_ref = opticpar_medium(1, & ! White ice
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/17.10, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_blueice_ref = opticpar_medium(1, & ! Blue ice
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/8.40, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_drysnow_ref = opticpar_medium(1, & ! Dry snow
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/25.00, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_meltingsnow_ref = opticpar_medium(1, & ! Melting snow
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/15.00, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_ice_opaque = opticpar_medium(1, & ! Opaque ice
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/1.0E+070, (1.E+100,i=2,nband_optic_max)/)) , &
! opticpar_snow_opaque = opticpar_medium(1, & ! Opaque snow
! (/1.0, (0.0,i=2,nband_optic_max)/), &
! (/1.0E+070, (1.E+100,i=2,nband_optic_max)/))
TYPE (opticpar_medium), PARAMETER :: &
opticpar_water_ref = opticpar_medium(1, & ! Water (reference)
(/1., (0.,i=2,nband_optic_max)/), &
(/3., (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_water_trans = opticpar_medium(2, & ! Transparent Water (two-band)
(/0.10, 0.90, (0.,i=3,nband_optic_max)/), &
(/2.0, 0.20, (1.E+10,i=3,nband_optic_max)/)) , &
opticpar_whiteice_ref = opticpar_medium(1, & ! White ice
(/1., (0.,i=2,nband_optic_max)/), &
(/17.1, (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_blueice_ref = opticpar_medium(1, & ! Blue ice
(/1., (0.,i=2,nband_optic_max)/), &
(/8.4, (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_drysnow_ref = opticpar_medium(1, & ! Dry snow
(/1., (0.,i=2,nband_optic_max)/), &
(/25.0, (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_meltingsnow_ref = opticpar_medium(1, & ! Melting snow
(/1., (0.,i=2,nband_optic_max)/), &
(/15.0, (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_ice_opaque = opticpar_medium(1, & ! Opaque ice
(/1., (0.,i=2,nband_optic_max)/), &
(/1.0E+07, (1.E+10,i=2,nband_optic_max)/)) , &
opticpar_snow_opaque = opticpar_medium(1, & ! Opaque snow
(/1., (0.,i=2,nband_optic_max)/), &
(/1.0E+07, (1.E+10,i=2,nband_optic_max)/))
!==============================================================================
END MODULE flake_paramoptic_ref
!------------------------------------------------------------------------------
MODULE flake_parameters
!------------------------------------------------------------------------------
!
! Description:
!
! Values of empirical constants of the lake model FLake
! and of several thermodynamic parameters are set.
!
USE data_parameters , ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Dimensionless constants
! in the equations for the mixed-layer depth
! and for the shape factor with respect to the temperature profile in the thermocline
REAL (KIND = kind_phys), PARAMETER :: &
! c_cbl_1 = 0.170 , & ! Constant in the CBL entrainment equation
! c_cbl_2 = 1.0 , & ! Constant in the CBL entrainment equation
! c_sbl_ZM_n = 0.50 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
! c_sbl_ZM_s = 10.0 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
! c_sbl_ZM_i = 20.0 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
! c_relax_h = 0.0300 , & ! Constant in the relaxation equation for the SBL depth
! c_relax_C = 0.00300 ! Constant in the relaxation equation for the shape factor
! with respect to the temperature profile in the thermocline
c_cbl_1 = 0.17 , & ! Constant in the CBL entrainment equation
c_cbl_2 = 1. , & ! Constant in the CBL entrainment equation
c_sbl_ZM_n = 0.5 , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
c_sbl_ZM_s = 10. , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
c_sbl_ZM_i = 20. , & ! Constant in the ZM1996 equation for the equilibrium SBL depth
c_relax_h = 0.030 , & ! Constant in the relaxation equation for the SBL depth
c_relax_C = 0.0030
! Parameters of the shape functions
! Indices refer to T - thermocline, S - snow, I - ice,
! B1 - upper layer of the bottom sediments, B2 - lower layer of the bottom sediments.
! "pr0" and "pr1" denote zeta derivatives of the corresponding shape function
! at "zeta=0" ad "zeta=1", respectively.
REAL (KIND = kind_phys), PARAMETER :: &
C_T_min = 0.5 , & ! Minimum value of the shape factor C_T (thermocline)
C_T_max = 0.8 , & ! Maximum value of the shape factor C_T (thermocline)
Phi_T_pr0_1 = 40.0/3.0 , & ! Constant in the expression for the T shape-function derivative
Phi_T_pr0_2 = 20.0/3.0 , & ! Constant in the expression for the T shape-function derivative
C_TT_1 = 11.0/18.0 , & ! Constant in the expression for C_TT (thermocline)
C_TT_2 = 7.0/45.0 , & ! Constant in the expression for C_TT (thermocline)
C_B1 = 2.0/3.0 , & ! Shape factor (upper layer of bottom sediments)
C_B2 = 3.0/5.0 , & ! Shape factor (lower layer of bottom sediments)
Phi_B1_pr0 = 2.0 , & ! B1 shape-function derivative
C_S_lin = 0.5 , & ! Shape factor (linear temperature profile in the snow layer)
Phi_S_pr0_lin = 1.0 , & ! S shape-function derivative (linear profile)
C_I_lin = 0.5 , & ! Shape factor (linear temperature profile in the ice layer)
Phi_I_pr0_lin = 1.0 , & ! I shape-function derivative (linear profile)
Phi_I_pr1_lin = 1.0 , & ! I shape-function derivative (linear profile)
Phi_I_ast_MR = 2.0 , & ! Constant in the MR2004 expression for I shape factor
C_I_MR = 1.0/12.0 , & ! Constant in the MR2004 expression for I shape factor
H_Ice_max = 3.0 ! Maximum ice tickness in
! the Mironov and Ritter (2004, MR2004) ice model [m]
! Security constants
REAL (KIND = kind_phys), PARAMETER :: &
h_Snow_min_flk = 1.0E-5 , & ! Minimum snow thickness [m]
h_Ice_min_flk = 1.0E-9 , & ! Minimum ice thickness [m]
h_ML_min_flk = 1.0E-2 , & ! Minimum mixed-layer depth [m]
h_ML_max_flk = 1.0E+3 , & ! Maximum mixed-layer depth [m]
H_B1_min_flk = 1.0E-3 , & ! Minimum thickness of the upper layer of bottom sediments [m]
u_star_min_flk = 1.0E-6 ! Minimum value of the surface friction velocity [m s^{-1}]
! Security constant(s)
REAL (KIND = kind_phys), PARAMETER :: &
c_small_flk = 1.0E-10 ! A small number
! Thermodynamic parameters
REAL (KIND = kind_phys), PARAMETER :: &
tpl_grav = 9.81 , & ! Acceleration due to gravity [m s^{-2}]
tpl_T_r = 277.13 , & ! Temperature of maximum density of fresh water [K]
tpl_T_f = 273.15 , & ! Fresh water freezing point [K]
tpl_a_T = 1.6509E-05 , & ! Constant in the fresh-water equation of state [K^{-2}]
tpl_rho_w_r = 1.0E+03 , & ! Maximum density of fresh water [kg m^{-3}]
tpl_rho_I = 9.1E+02 , & ! Density of ice [kg m^{-3}]
tpl_rho_S_min = 1.0E+02 , & ! Minimum snow density [kg m^{-3}]
tpl_rho_S_max = 4.0E+02 , & ! Maximum snow density [kg m^{-3}]
tpl_Gamma_rho_S = 2.0E+02 , & ! Empirical parameter [kg m^{-4}]
! in the expression for the snow density
tpl_L_f = 3.3E+05 , & ! Latent heat of fusion [J kg^{-1}]
tpl_c_w = 4.2E+03 , & ! Specific heat of water [J kg^{-1} K^{-1}]
tpl_c_I = 2.1E+03 , & ! Specific heat of ice [J kg^{-1} K^{-1}]
tpl_c_S = 2.1E+03 , & ! Specific heat of snow [J kg^{-1} K^{-1}]
tpl_kappa_w = 5.46E-01 , & ! Molecular heat conductivity of water [J m^{-1} s^{-1} K^{-1}]
tpl_kappa_I = 2.29 , & ! Molecular heat conductivity of ice [J m^{-1} s^{-1} K^{-1}]
tpl_kappa_S_min = 0.2 , & ! Minimum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
tpl_kappa_S_max = 1.5 , & ! Maximum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
tpl_Gamma_kappa_S = 1.3 ! Empirical parameter [J m^{-2} s^{-1} K^{-1}]
! in the expression for the snow heat conductivity
!==============================================================================
END MODULE flake_parameters
!------------------------------------------------------------------------------
MODULE flake
!------------------------------------------------------------------------------
!
! Description:
!
! The main program unit of the lake model FLake,
! containing most of the FLake procedures.
! Most FLake variables and local parameters are declared.
!
! FLake (Fresh-water Lake) is a lake model capable of predicting the surface temperature
! in lakes of various depth on the time scales from a few hours to a year.
! The model is based on a two-layer parametric representation of
! the evolving temperature profile, where the structure of the stratified layer between the
! upper mixed layer and the basin bottom, the lake thermocline,
! is described using the concept of self-similarity of the temperature-depth curve.
! The concept was put forward by Kitaigorodskii and Miropolsky (1970)
! to describe the vertical temperature structure of the oceanic seasonal thermocline.
! It has since been successfully used in geophysical applications.
! The concept of self-similarity of the evolving temperature profile
! is also used to describe the vertical structure of the thermally active upper layer
! of bottom sediments and of the ice and snow cover.
!
! The lake model incorporates the heat budget equations
! for the four layers in question, viz., snow, ice, water and bottom sediments,
! developed with due regard for the vertically distributed character
! of solar radiation heating.
! The entrainment equation that incorporates the Zilitinkevich (1975) spin-up term
! is used to compute the depth of a convectively-mixed layer.
! A relaxation-type equation is used
! to compute the wind-mixed layer depth in stable and neutral stratification,
! where a multi-limit formulation for the equilibrium mixed-layer depth
! proposed by Zilitinkevich and Mironov (1996)
! accounts for the effects of the earth's rotation, of the surface buoyancy flux
! and of the static stability in the thermocline.
! The equations for the mixed-layer depth are developed with due regard for
! the volumetric character of the radiation heating.
! Simple thermodynamic arguments are invoked to develop
! the evolution equations for the ice thickness and for the snow thickness.
! The heat flux through the water-bottom sediment interface is computed,
! using a parameterization proposed by Golosov et al. (1998).
! The heat flux trough the air-water interface
! (or through the air-ice or air-snow interface)
! is provided by the driving atmospheric model.
!
! Empirical constants and parameters of the lake model
! are estimated, using independent empirical and numerical data.
! They should not be re-evaluated when the model is applied to a particular lake.
! The only lake-specific parameters are the lake depth,
! the optical characteristics of lake water,
! the temperature at the bottom of the thermally active layer
! of bottom sediments and the depth of that layer.
!
! A detailed description of the lake model is given in
! Mironov, D. V., 2005:
! Parameterization of Lakes in Numerical Weather Prediction.
! Part 1: Description of a Lake Model.
! Manuscript is available from the author.
! Dmitrii Mironov
! German Weather Service, Kaiserleistr. 29/35, D-63067 Offenbach am Main, Germany.
! dmitrii.mironov@dwd.de
!
! Lines embraced with "!_tmp" contain temporary parts of the code.
! Lines embraced/marked with "!_dev" may be replaced
! as improved parameterizations are developed and tested.
! Lines embraced/marked with "!_dm" are DM's comments
! that may be helpful to a user.
! Lines embraced/marked with "!_dbg" are used
! for debugging purposes only.
!
USE data_parameters , ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
!
! The variables declared below
! are accessible to all program units of the MODULE flake.
! Some of them should be USEd by the driving routines that call flake routines.
! These are basically the quantities computed by FLake.
! All variables declared below have a suffix "flk".
! FLake variables of type REAL
! Temperatures at the previous time step ("p") and the updated temperatures ("n")
REAL (KIND = kind_phys) :: &
T_mnw_p_flk, T_mnw_n_flk , & ! Mean temperature of the water column [K]
T_snow_p_flk, T_snow_n_flk , & ! Temperature at the air-snow interface [K]
T_ice_p_flk, T_ice_n_flk , & ! Temperature at the snow-ice or air-ice interface [K]
T_wML_p_flk, T_wML_n_flk , & ! Mixed-layer temperature [K]
T_bot_p_flk, T_bot_n_flk , & ! Temperature at the water-bottom sediment interface [K]
T_B1_p_flk, T_B1_n_flk ! Temperature at the bottom of the upper layer of the sediments [K]
! Thickness of various layers at the previous time step ("p") and the updated values ("n")
REAL (KIND = kind_phys) :: &
h_snow_p_flk, h_snow_n_flk , & ! Snow thickness [m]
h_ice_p_flk, h_ice_n_flk , & ! Ice thickness [m]
h_ML_p_flk, h_ML_n_flk , & ! Thickness of the mixed-layer [m]
H_B1_p_flk, H_B1_n_flk ! Thickness of the upper layer of bottom sediments [m]
! The shape factor(s) at the previous time step ("p") and the updated value(s) ("n")
REAL (KIND = kind_phys) :: &
C_T_p_flk, C_T_n_flk , & ! Shape factor (thermocline)
C_TT_flk , & ! Dimensionless parameter (thermocline)
C_Q_flk , & ! Shape factor with respect to the heat flux (thermocline)
C_I_flk , & ! Shape factor (ice)
C_S_flk ! Shape factor (snow)
! Derivatives of the shape functions
REAL (KIND = kind_phys) :: &
Phi_T_pr0_flk , & ! d\Phi_T(0)/d\zeta (thermocline)
Phi_I_pr0_flk , & ! d\Phi_I(0)/d\zeta_I (ice)
Phi_I_pr1_flk , & ! d\Phi_I(1)/d\zeta_I (ice)
Phi_S_pr0_flk ! d\Phi_S(0)/d\zeta_S (snow)
! Heat and radiation fluxes
REAL (KIND = kind_phys) :: &
Q_snow_flk , & ! Heat flux through the air-snow interface [W m^{-2}]
Q_ice_flk , & ! Heat flux through the snow-ice or air-ice interface [W m^{-2}]
Q_w_flk , & ! Heat flux through the ice-water or air-water interface [W m^{-2}]
Q_bot_flk , & ! Heat flux through the water-bottom sediment interface [W m^{-2}]
I_atm_flk , & ! Radiation flux at the lower boundary of the atmosphere [W m^{-2}],
! i.e. the incident radiation flux with no regard for the surface albedo.
I_snow_flk , & ! Radiation flux through the air-snow interface [W m^{-2}]
I_ice_flk , & ! Radiation flux through the snow-ice or air-ice interface [W m^{-2}]
I_w_flk , & ! Radiation flux through the ice-water or air-water interface [W m^{-2}]
I_h_flk , & ! Radiation flux through the mixed-layer-thermocline interface [W m^{-2}]
I_bot_flk , & ! Radiation flux through the water-bottom sediment interface [W m^{-2}]
I_intm_0_h_flk , & ! Mean radiation flux over the mixed layer [W m^{-1}]
I_intm_h_D_flk , & ! Mean radiation flux over the thermocline [W m^{-1}]
I_intm_D_H_flk , & ! Mean radiation flux over the deeper layer defined by Shaobo Zhang [W m^{-1}]
I_HH_flk , & ! Radiation flux through the bottom of the deeper layer defined by Shaobo Zhang [W m^{-2}]
Q_star_flk ! A generalized heat flux scale [W m^{-2}]
! Velocity scales
REAL (KIND = kind_phys) :: &
u_star_w_flk , & ! Friction velocity in the surface layer of lake water [m s^{-1}]
w_star_sfc_flk ! Convective velocity scale,
! using a generalized heat flux scale [m s^{-1}]
! The rate of snow accumulation
REAL (KIND = kind_phys) :: &
dMsnowdt_flk ! The rate of snow accumulation [kg m^{-2} s^{-1}]
! The secondary layer temp
REAL (KIND = kind_phys) :: &
T_BOT_2_IN_FLK
!==============================================================================
! Procedures
!==============================================================================
CONTAINS
!==============================================================================
! The codes of the FLake procedures are stored in separate "*.incf" files
! and are included below.
!------------------------------------------------------------------------------
!==============================================================================
! include 'flake_radflux.incf'
!------------------------------------------------------------------------------
! changed by Shaobo Zhang
SUBROUTINE flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
opticpar_water, opticpar_ice, opticpar_snow, &
depth_bs )
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the radiation fluxes
! at the snow-ice, ice-water, air-water,
! mixed layer-thermocline and water column-bottom sediment interfaces,
! the mean radiation flux over the mixed layer,
! and the mean radiation flux over the thermocline.
!
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "flake".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
USE flake_derivedtypes ! Definitions of derived TYPEs
USE flake_parameters , ONLY : &
h_Snow_min_flk , & ! Minimum snow thickness [m]
h_Ice_min_flk , & ! Minimum ice thickness [m]
h_ML_min_flk ! Minimum mixed-layer depth [m]
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (procedure arguments)
REAL (KIND = kind_phys), INTENT(IN) :: &
depth_w , & ! The lake depth [m]
depth_bs , & ! The depth_bs added by Shaobo Zhang
albedo_water , & ! Albedo of the water surface
albedo_ice , & ! Albedo of the ice surface
albedo_snow ! Albedo of the snow surface
TYPE (opticpar_medium), INTENT(IN) :: &
opticpar_water , & ! Optical characteristics of water
opticpar_ice , & ! Optical characteristics of ice
opticpar_snow ! Optical characteristics of snow
! Local variables of type INTEGER
INTEGER (KIND = iintegers) :: & ! Help variable(s)
i ! DO loop index
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
IF(h_ice_p_flk.GE.h_Ice_min_flk) THEN ! Ice exists
IF(h_snow_p_flk.GE.h_Snow_min_flk) THEN ! There is snow above the ice
I_snow_flk = I_atm_flk*(1.0-albedo_snow)
I_bot_flk = 0.0
DO i=1, opticpar_snow%nband_optic
I_bot_flk = I_bot_flk + &
opticpar_snow%frac_optic(i)*EXP(-opticpar_snow%extincoef_optic(i)*h_snow_p_flk)
END DO
I_ice_flk = I_snow_flk*I_bot_flk
ELSE ! No snow above the ice
I_snow_flk = I_atm_flk
I_ice_flk = I_atm_flk*(1.0-albedo_ice)
END IF
I_bot_flk = 0.0
DO i=1, opticpar_ice%nband_optic
I_bot_flk = I_bot_flk + &
opticpar_ice%frac_optic(i)*EXP(-opticpar_ice%extincoef_optic(i)*h_ice_p_flk)
END DO
I_w_flk = I_ice_flk*I_bot_flk
ELSE ! No ice-snow cover
I_snow_flk = I_atm_flk
I_ice_flk = I_atm_flk
I_w_flk = I_atm_flk*(1.0-albedo_water)
END IF
IF(h_ML_p_flk.GE.h_ML_min_flk) THEN ! Radiation flux at the bottom of the mixed layer
I_bot_flk = 0.0
DO i=1, opticpar_water%nband_optic
I_bot_flk = I_bot_flk + &
opticpar_water%frac_optic(i)*EXP(-opticpar_water%extincoef_optic(i)*h_ML_p_flk)
! print*,'nband_optic=',opticpar_water%nband_optic
! print*,'Extinction=',opticpar_water%extincoef_optic(i)
END DO
I_h_flk = I_w_flk*I_bot_flk
ELSE ! Mixed-layer depth is less then a minimum value
I_h_flk = I_w_flk
END IF
I_bot_flk = 0.0 ! Radiation flux at the lake bottom
DO i=1, opticpar_water%nband_optic
I_bot_flk = I_bot_flk + &
opticpar_water%frac_optic(i)*EXP(-opticpar_water%extincoef_optic(i)*depth_w)
END DO
I_bot_flk = I_w_flk*I_bot_flk
IF(h_ML_p_flk.GE.h_ML_min_flk) THEN ! Integral-mean radiation flux over the mixed layer
I_intm_0_h_flk = 0.0
DO i=1, opticpar_water%nband_optic
I_intm_0_h_flk = I_intm_0_h_flk + &
opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
(1.0 - EXP(-opticpar_water%extincoef_optic(i)*h_ML_p_flk))
END DO
I_intm_0_h_flk = I_w_flk*I_intm_0_h_flk/h_ML_p_flk
ELSE
I_intm_0_h_flk = I_h_flk
END IF
IF(h_ML_p_flk.LE.depth_w-h_ML_min_flk) THEN ! Integral-mean radiation flux over the thermocline
I_intm_h_D_flk = 0.0
DO i=1, opticpar_water%nband_optic
I_intm_h_D_flk = I_intm_h_D_flk + &
opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
( EXP(-opticpar_water%extincoef_optic(i)*h_ML_p_flk) &
- EXP(-opticpar_water%extincoef_optic(i)*depth_w) )
END DO
I_intm_h_D_flk = I_w_flk*I_intm_h_D_flk/(depth_w-h_ML_p_flk)
ELSE
I_intm_h_D_flk = I_h_flk
END IF
! Added by Shaobo Zhang
IF(depth_bs.GE.h_ML_min_flk) THEN! Integral-mean radiation flux over the deeper layer defined by Shaobo Zhang
I_intm_D_H_flk = 0.0
DO i=1, opticpar_water%nband_optic
I_intm_D_H_flk = I_intm_D_H_flk + &
opticpar_water%frac_optic(i)/opticpar_water%extincoef_optic(i)* &
( EXP(-opticpar_water%extincoef_optic(i)*depth_w) &
- EXP(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs)) )
END DO
I_intm_D_H_flk = I_w_flk*I_intm_D_H_flk/depth_bs
ELSE
I_intm_D_H_flk = I_bot_flk
END IF
! Radiation flux at the bottom of the deeper layer defined by Shaobo Zhang
I_HH_flk = 0.0
DO i=1, opticpar_water%nband_optic
I_HH_flk = I_HH_flk + &
opticpar_water%frac_optic(i)*EXP(-opticpar_water%extincoef_optic(i)*(depth_w+depth_bs))
END DO
I_HH_flk = I_w_flk*I_HH_flk
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END SUBROUTINE flake_radflux
!==============================================================================
!==============================================================================
! include 'flake_main.incf'
!------------------------------------------------------------------------------
SUBROUTINE flake_main ( depthw, depthbs, T_bs, par_Coriolis, &
extincoef_water_typ, &
del_time, T_sfc_p, T_sfc_n, T_bot_2_in, &
T_bot_2_out )
!------------------------------------------------------------------------------
!
! Description:
!
! The main driving routine of the lake model FLake
! where computations are performed.
! Advances the surface temperature
! and other FLake variables one time step.
! At the moment, the Euler explicit scheme is used.
!
! Lines embraced with "!_tmp" contain temporary parts of the code.
! Lines embraced/marked with "!_dev" may be replaced
! as improved parameterizations are developed and tested.
! Lines embraced/marked with "!_dm" are DM's comments
! that may be helpful to a user.
! Lines embraced/marked with "!_dbg" are used
! for debugging purposes only.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "flake".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
USE flake_parameters ! Thermodynamic parameters and dimensionless constants of FLake
USE flake_configure ! Switches and parameters that configure FLake
use machine, only: kind_phys
! ADDED by Shaobo Zhang
! USE mod_dynparam, only : lake_depth_max
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (procedure arguments)
! changed by Shaobo Zhang
REAL (KIND = kind_phys), INTENT(IN) :: &
depthw , & ! The lake depth [m]
depthbs , & ! Depth of the thermally active layer of bottom sediments [m]
T_bs , & ! Temperature at the outer edge of
! the thermally active layer of bottom sediments [K]
par_Coriolis , & ! The Coriolis parameter [s^{-1}]
extincoef_water_typ , & ! "Typical" extinction coefficient of the lake water [m^{-1}],
! used to compute the equilibrium CBL depth
del_time , & ! The model time step [s]
T_sfc_p , & ! Surface temperature at the previous time step [K]
T_bot_2_in
REAL (KIND = kind_phys) :: &
depth_w , & ! The lake depth [m]
depth_bs ! Depth of the thermally active layer of bottom sediments [m]
! Output (procedure arguments)
REAL (KIND = kind_phys), INTENT(OUT) :: &
T_sfc_n , & ! Updated surface temperature [K]
! (equal to the updated value of either T_ice, T_snow or T_wML)
T_bot_2_out
! Local variables of type LOGICAL
LOGICAL :: &
l_ice_create , & ! Switch, .TRUE. = ice does not exist but should be created
l_snow_exists , & ! Switch, .TRUE. = there is snow above the ice
l_ice_meltabove ! Switch, .TRUE. = snow/ice melting from above takes place
! Local variables of type INTEGER
INTEGER (KIND = iintegers) :: &
i ! Loop index
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
d_T_mnw_dt , & ! Time derivative of T_mnw [K s^{-1}]
d_T_ice_dt , & ! Time derivative of T_ice [K s^{-1}]
d_T_bot_dt , & ! Time derivative of T_bot [K s^{-1}]
d_T_B1_dt , & ! Time derivative of T_B1 [K s^{-1}]
d_h_snow_dt , & ! Time derivative of h_snow [m s^{-1}]
d_h_ice_dt , & ! Time derivative of h_ice [m s^{-1}]
d_h_ML_dt , & ! Time derivative of h_ML [m s^{-1}]
d_H_B1_dt , & ! Time derivative of H_B1 [m s^{-1}]
d_h_D_dt , & ! Time derivative of h_D, new defined by Shaobo Zhang
d_T_H_dt , & ! Time derivative of T_H, new defined by Shaobo Zhang
d_C_T_dt ! Time derivative of C_T [s^{-1}]
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
N_T_mean , & ! The mean buoyancy frequency in the thermocline [s^{-1}]
tmp , & ! temperary variable
ZM_h_scale , & ! The ZM96 equilibrium SBL depth scale [m]
conv_equil_h_scale ! The equilibrium CBL depth scale [m]
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
h_ice_threshold , & ! If h_ice<h_ice_threshold, use quasi-equilibrium ice model
flk_str_1 , & ! Help storage variable
flk_str_2 , & ! Help storage variable
R_H_icesnow , & ! Dimensionless ratio, used to store intermediate results
R_rho_c_icesnow , & ! Dimensionless ratio, used to store intermediate results
R_TI_icesnow , & ! Dimensionless ratio, used to store intermediate results
R_Tstar_icesnow ! Dimensionless ratio, used to store intermediate results
! ADDED by Shaobo Zhang
!REAL (KIND = kind_phys) :: T_bot_2_in, T_bot_2_out
REAL (KIND = kind_phys) :: CT, CTT, CQ
REAL (KIND = kind_phys) :: lake_depth_max
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
!_dm
! Security. Set time-rate-of-change of prognostic variables to zero.
! Set prognostic variables to their values at the previous time step.
! (This is to avoid spurious changes of prognostic variables
! when FLake is used within a 3D model, e.g. to avoid spurious generation of ice
! at the neighbouring lake points as noticed by Burkhardt Rockel.)
!_dm
! print*,'T_sfc_p=',T_sfc_p,' T_bot_2_in=',T_bot_2_in
! print*,'h_snow_p_flk=',h_snow_p_flk,' h_ice_p_flk=',h_ice_p_flk
! print*,'T_snow_p_flk=',T_snow_p_flk, ' T_ice_p_flk=',T_ice_p_flk
! print*,'T_wML_p_flk=',T_wML_p_flk, ' T_mnw_p_flk=',T_mnw_p_flk
! print*,'T_bot_p_flk=',T_bot_p_flk, ' T_B1_p_flk=',T_B1_p_flk
! print*,'h_ML_p_flk=',h_ML_p_flk, ' H_B1_p_flk=',H_B1_p_flk
! print*,'C_T_p_flk=',C_T_p_flk
! print*,'depthw= ', depthw
! print*,'depthbs=',depthbs
lake_depth_max = 60.0
depth_w =depthw
depth_bs = depthbs
d_T_mnw_dt = 0.0
d_T_ice_dt = 0.0
d_T_bot_dt = 0.0
d_T_B1_dt = 0.0
d_h_snow_dt = 0.0
d_h_ice_dt = 0.0
d_h_ML_dt = 0.0
d_H_B1_dt = 0.0
d_C_T_dt = 0.0
T_snow_n_flk = T_snow_p_flk
T_ice_n_flk = T_ice_p_flk
T_wML_n_flk = T_wML_p_flk
T_mnw_n_flk = T_mnw_p_flk
T_bot_n_flk = T_bot_p_flk
T_B1_n_flk = T_B1_p_flk
h_snow_n_flk = h_snow_p_flk
h_ice_n_flk = h_ice_p_flk
h_ML_n_flk = h_ML_p_flk
H_B1_n_flk = H_B1_p_flk
C_T_n_flk = C_T_p_flk
!------------------------------------------------------------------------------
! Compute fluxes, using variables from the previous time step.
!------------------------------------------------------------------------------
!_dm
! At this point, the heat and radiation fluxes, namely,
! Q_snow_flk, Q_ice_flk, Q_w_flk,
! I_atm_flk, I_snow_flk, I_ice_flk, I_w_flk, I_h_flk, I_bot_flk,
! the mean radiation flux over the mixed layer, I_intm_0_h_flk,
! and the mean radiation flux over the thermocline, I_intm_h_D_flk,
! should be known.
! They are computed within "flake_interface" (or within the driving model)
! and are available to "flake_main"
! through the above variables declared in the MODULE "flake".
! In case a lake is ice-covered, Q_w_flk is re-computed below.
!_dm
! Heat flux through the ice-water interface
IF(h_ice_p_flk.GE.h_Ice_min_flk) THEN ! Ice exists
IF(h_ML_p_flk.LE.h_ML_min_flk) THEN ! Mixed-layer depth is zero, compute flux
Q_w_flk = -tpl_kappa_w*(T_bot_p_flk-T_wML_p_flk)/depth_w ! Flux with linear T(z)
Phi_T_pr0_flk = Phi_T_pr0_1*C_T_p_flk-Phi_T_pr0_2 ! d\Phi(0)/d\zeta (thermocline)
Q_w_flk = Q_w_flk*MAX(Phi_T_pr0_flk, 1.0) ! Account for an increased d\Phi(0)/d\zeta
ELSE
Q_w_flk = 0.0 ! Mixed-layer depth is greater than zero, set flux to zero
END IF
ELSE ! If Ice doesn't exist, set flux to zero, YWu 2019
Q_w_flk = 0.0
END IF
! A generalized heat flux scale
Q_star_flk = Q_w_flk + I_w_flk + I_h_flk - 2.0*I_intm_0_h_flk
! changed by Shaobo Zhang
! Heat flux through the water-bottom sediment interface
!IF(lflk_botsed_use) THEN
! 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
!ELSE
! Q_bot_flk = 0.0 ! The bottom-sediment scheme is not used
!END IF
IF(lflk_botsed_use) THEN ! The bottom-sediment scheme is used
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
ELSE ! The scheme written by Shaobo Zhang for the deeper layer of a deep lake is used
Q_bot_flk = -tpl_kappa_w*(T_bot_2_in-T_bot_p_flk)/MAX(depth_bs, H_B1_min_flk)*Phi_B1_pr0
END IF
!------------------------------------------------------------------------------
! Check if ice exists or should be created.
! If so, compute the thickness and the temperature of ice and snow.
!------------------------------------------------------------------------------
!_dm
! Notice that a quasi-equilibrium ice-snow model is used
! to avoid numerical instability when the ice is thin.
! This is always the case when new ice is created.
!_dm
!_dev
! The dependence of snow density and of snow heat conductivity
! on the snow thickness is accounted for parametrically.
! That is, the time derivatives of \rho_S and \kappa_S are neglected.
! The exception is the equation for the snow thickness
! in case of snow accumulation and no melting,
! where d\rho_S/dt is incorporated.
! Furthermore, some (presumably small) correction terms incorporating
! the snow density and the snow heat conductivity are dropped out.
! Those terms may be included as better formulations
! for \rho_S and \kappa_S are available.
!_dev
! Default values
l_ice_create = .FALSE.
l_ice_meltabove = .FALSE.
Ice_exist: IF(h_ice_p_flk.LT.h_Ice_min_flk) THEN ! Ice does not exist
l_ice_create = T_wML_p_flk.LE.(tpl_T_f+c_small_flk).AND.Q_w_flk.LT.0.0
IF(l_ice_create) THEN ! Ice does not exist but should be created
d_h_ice_dt = -Q_w_flk/tpl_rho_I/tpl_L_f
h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time ! Advance h_ice
T_ice_n_flk = tpl_T_f + h_ice_n_flk*Q_w_flk/tpl_kappa_I/Phi_I_pr0_lin ! Ice temperature
d_h_snow_dt = dMsnowdt_flk/tpl_rho_S_min
h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
Phi_I_pr1_flk = Phi_I_pr1_lin &
+ Phi_I_ast_MR*MIN(1.0, h_ice_n_flk/H_Ice_max) ! d\Phi_I(1)/d\zeta_I (ice)
R_H_icesnow = Phi_I_pr1_flk/Phi_S_pr0_lin*tpl_kappa_I/flake_snowheatconduct(h_snow_n_flk) &
* h_snow_n_flk/MAX(h_ice_n_flk, h_Ice_min_flk)
T_snow_n_flk = T_ice_n_flk + R_H_icesnow*(T_ice_n_flk-tpl_T_f) ! Snow temperature
END IF
ELSE Ice_exist ! Ice exists
l_snow_exists = h_snow_p_flk.GE.h_Snow_min_flk ! Check if there is snow above the ice
Melting: IF(T_snow_p_flk.GE.(tpl_T_f-c_small_flk)) THEN ! T_sfc = T_f, check for melting from above
! T_snow = T_ice if snow is absent
IF(l_snow_exists) THEN ! There is snow above the ice
flk_str_1 = Q_snow_flk + I_snow_flk - I_ice_flk ! Atmospheric forcing
IF(flk_str_1.GE.0.0) THEN ! Melting of snow and ice from above
l_ice_meltabove = .TRUE.
d_h_snow_dt = (-flk_str_1/tpl_L_f+dMsnowdt_flk)/flake_snowdensity(h_snow_p_flk)
d_h_ice_dt = -(I_ice_flk - I_w_flk - Q_w_flk)/tpl_L_f/tpl_rho_I
END IF
ELSE ! No snow above the ice
flk_str_1 = Q_ice_flk + I_ice_flk - I_w_flk - Q_w_flk ! Atmospheric forcing + heating from the water
IF(flk_str_1.GE.0.0) THEN ! Melting of ice from above, snow accumulation may occur
l_ice_meltabove = .TRUE.
d_h_ice_dt = -flk_str_1/tpl_L_f/tpl_rho_I
d_h_snow_dt = dMsnowdt_flk/tpl_rho_S_min
END IF
END IF
IF(l_ice_meltabove) THEN ! Melting from above takes place
h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time ! Advance h_ice
h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
T_ice_n_flk = tpl_T_f ! Set T_ice to the freezing point
T_snow_n_flk = tpl_T_f ! Set T_snow to the freezing point
END IF
END IF Melting
No_Melting: IF(.NOT.l_ice_meltabove) THEN ! No melting from above
d_h_snow_dt = flake_snowdensity(h_snow_p_flk)
IF(d_h_snow_dt.LT.tpl_rho_S_max) THEN ! Account for d\rho_S/dt
flk_str_1 = h_snow_p_flk*tpl_Gamma_rho_S/tpl_rho_w_r
flk_str_1 = flk_str_1/(1.0-flk_str_1)
ELSE ! Snow density is equal to its maximum value, d\rho_S/dt=0
flk_str_1 = 0.0
END IF
d_h_snow_dt = dMsnowdt_flk/d_h_snow_dt/(1.0+flk_str_1) ! Snow accumulation
h_snow_n_flk = h_snow_p_flk + d_h_snow_dt*del_time ! Advance h_snow
Phi_I_pr0_flk = h_ice_p_flk/H_Ice_max ! h_ice relative to its maximum value
C_I_flk = C_I_lin - C_I_MR*(1.0+Phi_I_ast_MR)*Phi_I_pr0_flk ! Shape factor (ice)
Phi_I_pr1_flk = Phi_I_pr1_lin + Phi_I_ast_MR*Phi_I_pr0_flk ! d\Phi_I(1)/d\zeta_I (ice)
Phi_I_pr0_flk = Phi_I_pr0_lin - Phi_I_pr0_flk ! d\Phi_I(0)/d\zeta_I (ice)
h_ice_threshold = MAX(1.0, 2.0*C_I_flk*tpl_c_I*(tpl_T_f-T_ice_p_flk)/tpl_L_f)
h_ice_threshold = Phi_I_pr0_flk/C_I_flk*tpl_kappa_I/tpl_rho_I/tpl_c_I*h_ice_threshold
h_ice_threshold = SQRT(h_ice_threshold*del_time) ! Threshold value of h_ice
h_ice_threshold = MIN(0.9*H_Ice_max, MAX(h_ice_threshold, h_Ice_min_flk))
! h_ice(threshold) < 0.9*H_Ice_max
IF(h_ice_p_flk.LT.h_ice_threshold) THEN ! Use a quasi-equilibrium ice model
IF(l_snow_exists) THEN ! Use fluxes at the air-snow interface
flk_str_1 = Q_snow_flk + I_snow_flk - I_w_flk
ELSE ! Use fluxes at the air-ice interface
flk_str_1 = Q_ice_flk + I_ice_flk - I_w_flk
END IF
d_h_ice_dt = -(flk_str_1-Q_w_flk)/tpl_L_f/tpl_rho_I
h_ice_n_flk = h_ice_p_flk + d_h_ice_dt *del_time ! Advance h_ice
T_ice_n_flk = tpl_T_f + h_ice_n_flk*flk_str_1/tpl_kappa_I/Phi_I_pr0_flk ! Ice temperature
ELSE ! Use a complete ice model
d_h_ice_dt = tpl_kappa_I*(tpl_T_f-T_ice_p_flk)/h_ice_p_flk*Phi_I_pr0_flk
d_h_ice_dt = (Q_w_flk+d_h_ice_dt)/tpl_L_f/tpl_rho_I
h_ice_n_flk = h_ice_p_flk + d_h_ice_dt*del_time ! Advance h_ice
R_TI_icesnow = tpl_c_I*(tpl_T_f-T_ice_p_flk)/tpl_L_f ! Dimensionless parameter
R_Tstar_icesnow = 1.0 - C_I_flk ! Dimensionless parameter
IF(l_snow_exists) THEN ! There is snow above the ice
R_H_icesnow = Phi_I_pr1_flk/Phi_S_pr0_lin*tpl_kappa_I/flake_snowheatconduct(h_snow_p_flk) &
* h_snow_p_flk/h_ice_p_flk
R_rho_c_icesnow = flake_snowdensity(h_snow_p_flk)*tpl_c_S/tpl_rho_I/tpl_c_I
!_dev
!_dm
! These terms should be included as an improved understanding of the snow scheme is gained,
! of the effect of snow density in particular.
!_dm
!_nu R_Tstar_icesnow = R_Tstar_icesnow &
!_nu + (1.0+C_S_lin*h_snow_p_flk/h_ice_p_flk)*R_H_icesnow*R_rho_c_icesnow
!_dev
R_Tstar_icesnow = R_Tstar_icesnow*R_TI_icesnow ! Dimensionless parameter
!_dev
!_nu R_Tstar_icesnow = R_Tstar_icesnow &
!_nu + (1.0-R_rho_c_icesnow)*tpl_c_I*T_ice_p_flk/tpl_L_f
!_dev
flk_str_2 = Q_snow_flk+I_snow_flk-I_w_flk ! Atmospheric fluxes
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
d_T_ice_dt = -(1.0-2.0*C_S_lin)*R_H_icesnow*(tpl_T_f-T_ice_p_flk) &
* tpl_c_S*dMsnowdt_flk ! Effect of snow accumulation
ELSE ! No snow above the ice
R_Tstar_icesnow = R_Tstar_icesnow*R_TI_icesnow ! Dimensionless parameter
flk_str_2 = Q_ice_flk+I_ice_flk-I_w_flk ! Atmospheric fluxes
flk_str_1 = C_I_flk*h_ice_p_flk
d_T_ice_dt = 0.0
END IF
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 &
* (1.0-R_Tstar_icesnow) ! Add flux due to heat conduction
d_T_ice_dt = d_T_ice_dt - R_Tstar_icesnow*Q_w_flk ! Add flux from water to ice
d_T_ice_dt = d_T_ice_dt + flk_str_2 ! Add atmospheric fluxes
d_T_ice_dt = d_T_ice_dt/tpl_rho_I/tpl_c_I ! Total forcing
d_T_ice_dt = d_T_ice_dt/flk_str_1 ! dT_ice/dt
T_ice_n_flk = T_ice_p_flk + d_T_ice_dt*del_time ! Advance T_ice
END IF
Phi_I_pr1_flk = MIN(1.0, h_ice_n_flk/H_Ice_max) ! h_ice relative to its maximum value
Phi_I_pr1_flk = Phi_I_pr1_lin + Phi_I_ast_MR*Phi_I_pr1_flk ! d\Phi_I(1)/d\zeta_I (ice)
R_H_icesnow = Phi_I_pr1_flk/Phi_S_pr0_lin*tpl_kappa_I/flake_snowheatconduct(h_snow_n_flk) &
*h_snow_n_flk/MAX(h_ice_n_flk, h_Ice_min_flk)
T_snow_n_flk = T_ice_n_flk + R_H_icesnow*(T_ice_n_flk-tpl_T_f) ! Snow temperature
END IF No_Melting
END IF Ice_exist
! Security, limit h_ice by its maximum value
h_ice_n_flk = MIN(h_ice_n_flk, H_Ice_max)
! Security, limit the ice and snow temperatures by the freezing point
T_snow_n_flk = MIN(T_snow_n_flk, tpl_T_f)
T_ice_n_flk = MIN(T_ice_n_flk, tpl_T_f)
!_tmp
! Security, avoid too low values (these constraints are used for debugging purposes)
! T_snow_n_flk = MAX(T_snow_n_flk, 73.15)
! T_ice_n_flk = MAX(T_ice_n_flk, 73.15)
! The lowest natural temperature ever
! directly recorded at ground level on Earth is -89.2 C (-128.6 F; 184.0 K)
! at the Soviet Vostok Station in Antarctica on 21 July, 1983 by ground
! measurements.
! https://en.wikipedia.org/wiki/Lowest_temperature_recorded_on_Earth
T_snow_n_flk = MAX(T_snow_n_flk, 184.0)
T_ice_n_flk = MAX(T_ice_n_flk, 184.0)
!_tmp
! Remove too thin ice and/or snow
IF(h_ice_n_flk.LT.h_Ice_min_flk) THEN ! Check ice
h_ice_n_flk = 0.0 ! Ice is too thin, remove it, and
T_ice_n_flk = tpl_T_f ! set T_ice to the freezing point.
h_snow_n_flk = 0.0 ! Remove snow when there is no ice, and
T_snow_n_flk = tpl_T_f ! set T_snow to the freezing point.
l_ice_create = .FALSE. ! "Exotic" case, ice has been created but proved to be too thin
ELSE IF(h_snow_n_flk.LT.h_Snow_min_flk) THEN ! Ice exists, check snow
h_snow_n_flk = 0.0 ! Snow is too thin, remove it,
T_snow_n_flk = T_ice_n_flk ! and set the snow temperature equal to the ice temperature.
END IF
!------------------------------------------------------------------------------
! Compute the mean temperature of the water column.
!------------------------------------------------------------------------------
IF(l_ice_create) Q_w_flk = 0.0 ! Ice has just been created, set Q_w to zero
d_T_mnw_dt = (Q_w_flk - Q_bot_flk + I_w_flk - I_bot_flk)/tpl_rho_w_r/tpl_c_w/depth_w
!print*,'d_T_mnw_dt= ',d_T_mnw_dt
!print*,'Q_w_flk=',Q_w_flk
!print*,'Q_bot_flk=',Q_bot_flk
!print*,'I_w_flk=',I_w_flk
!print*,'I_bot_flk=',I_bot_flk
!print*,'tpl_rho_w_r=',tpl_rho_w_r
!print*,'tpl_c_w=',tpl_c_w
!print*,'depth_w=',depth_w
T_mnw_n_flk = T_mnw_p_flk + d_T_mnw_dt*del_time ! Advance T_mnw
T_mnw_n_flk = MAX(T_mnw_n_flk, tpl_T_f) ! Limit T_mnw by the freezing point
!------------------------------------------------------------------------------
! Compute the mixed-layer depth, the mixed-layer temperature,
! the bottom temperature and the shape factor
! with respect to the temperature profile in the thermocline.
! Different formulations are used, depending on the regime of mixing.
!------------------------------------------------------------------------------
HTC_Water: IF(h_ice_n_flk.GE.h_Ice_min_flk) THEN ! Ice exists
T_mnw_n_flk = MIN(T_mnw_n_flk, tpl_T_r) ! Limit the mean temperature under the ice by T_r
T_wML_n_flk = tpl_T_f ! The mixed-layer temperature is equal to the freezing point
IF(l_ice_create) THEN ! Ice has just been created
IF(h_ML_p_flk.GE.depth_w-h_ML_min_flk) THEN ! h_ML=D when ice is created
h_ML_n_flk = 0.0 ! Set h_ML to zero
C_T_n_flk = C_T_min ! Set C_T to its minimum value
ELSE ! h_ML<D when ice is created
h_ML_n_flk = h_ML_p_flk ! h_ML remains unchanged
C_T_n_flk = C_T_p_flk ! C_T (thermocline) remains unchanged
END IF
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)
! Update the bottom temperature
ELSE IF(T_bot_p_flk.LT.tpl_T_r) THEN ! Ice exists and T_bot < T_r, molecular heat transfer
h_ML_n_flk = h_ML_p_flk ! h_ML remains unchanged
C_T_n_flk = C_T_p_flk ! C_T (thermocline) remains unchanged
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)
! Update the bottom temperature
ELSE ! Ice exists and T_bot = T_r, convection due to bottom heating
T_bot_n_flk = tpl_T_r ! T_bot is equal to the temperature of maximum density
IF(h_ML_p_flk.GE.c_small_flk) THEN ! h_ML > 0
C_T_n_flk = C_T_p_flk ! C_T (thermocline) remains unchanged
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)
h_ML_n_flk = MAX(h_ML_n_flk, 0.0) ! Update the mixed-layer depth
ELSE ! h_ML = 0
h_ML_n_flk = h_ML_p_flk ! h_ML remains unchanged
C_T_n_flk = (T_wML_n_flk-T_mnw_n_flk)/(T_wML_n_flk-T_bot_n_flk)
C_T_n_flk = MIN(C_T_max, MAX(C_T_n_flk, C_T_min)) ! Update the shape factor (thermocline)
END IF
END IF
T_bot_n_flk = MIN(T_bot_n_flk, tpl_T_r) ! Security, limit the bottom temperature by T_r
ELSE HTC_Water ! Open water
! Generalised buoyancy flux scale and convective velocity scale
flk_str_1 = flake_buoypar(T_wML_p_flk)*Q_star_flk/tpl_rho_w_r/tpl_c_w
IF(flk_str_1.LT.0.0) THEN
w_star_sfc_flk = (-flk_str_1*h_ML_p_flk)**(1.0/3.0) ! Convection
ELSE
w_star_sfc_flk = 0.0 ! Neutral or stable stratification
END IF
!_dm
! The equilibrium depth of the CBL due to surface cooling with the volumetric heating
! is not computed as a solution to the transcendental equation.
! Instead, an algebraic formula is used
! that interpolates between the two asymptotic limits.
!_dm
conv_equil_h_scale = -Q_w_flk/MAX(I_w_flk, c_small_flk)
IF(conv_equil_h_scale.GT.0.0 .AND. conv_equil_h_scale.LT.1.0 &
.AND. T_wML_p_flk.GT.tpl_T_r) THEN ! The equilibrium CBL depth scale is only used above T_r
conv_equil_h_scale = SQRT(6.0*conv_equil_h_scale) &
+ 2.0*conv_equil_h_scale/(1.0-conv_equil_h_scale)
conv_equil_h_scale = MIN(depth_w, conv_equil_h_scale/extincoef_water_typ)
! print*,'extincoef_water_typ=',extincoef_water_typ
ELSE
conv_equil_h_scale = 0.0 ! Set the equilibrium CBL depth to zero
END IF
! Mean buoyancy frequency in the thermocline
N_T_mean = flake_buoypar(0.5*(T_wML_p_flk+T_bot_p_flk))*(T_wML_p_flk-T_bot_p_flk)
IF(h_ML_p_flk.LE.depth_w-h_ML_min_flk) THEN
tmp = N_T_mean/(depth_w-h_ML_p_flk)
N_T_mean = SQRT(abs(tmp)) ! Compute N
ELSE
N_T_mean = 0.0 ! h_ML=D, set N to zero
END IF
! The rate of change of C_T
d_C_T_dt = MAX(w_star_sfc_flk, u_star_w_flk, u_star_min_flk)**2_iintegers
d_C_T_dt = N_T_mean*(depth_w-h_ML_p_flk)**2_iintegers &
/ c_relax_C/d_C_T_dt ! Relaxation time scale for C_T
d_C_T_dt = (C_T_max-C_T_min)/MAX(d_C_T_dt, c_small_flk) ! Rate-of-change of C_T
! Compute the shape factor and the mixed-layer depth,
! using different formulations for convection and wind mixing
C_TT_flk = C_TT_1*C_T_p_flk-C_TT_2 ! C_TT, using C_T at the previous time step
C_Q_flk = 2.0*C_TT_flk/C_T_p_flk ! C_Q using C_T at the previous time step
Mixing_regime: IF(flk_str_1.LT.0.0) THEN ! Convective mixing
C_T_n_flk = C_T_p_flk + d_C_T_dt*del_time ! Update C_T, assuming dh_ML/dt>0
C_T_n_flk = MIN(C_T_max, MAX(C_T_n_flk, C_T_min)) ! Limit C_T
d_C_T_dt = (C_T_n_flk-C_T_p_flk)/del_time ! Re-compute dC_T/dt
IF(h_ML_p_flk.LE.depth_w-h_ML_min_flk) THEN ! Compute dh_ML/dt
IF(h_ML_p_flk.LE.h_ML_min_flk) THEN ! Use a reduced entrainment equation (spin-up)
d_h_ML_dt = c_cbl_1/c_cbl_2*MAX(w_star_sfc_flk, c_small_flk)
!_dbg
! print*, ' FLake: reduced entrainment eq. D_time*d_h_ML_dt = ', d_h_ML_dt*del_time
! print*, ' w_* = ', w_star_sfc_flk
! print*, ' \beta*Q_* = ', flk_str_1
!_dbg
ELSE ! Use a complete entrainment equation
R_H_icesnow = depth_w/h_ML_p_flk
R_rho_c_icesnow = R_H_icesnow-1.0
R_TI_icesnow = C_T_p_flk/C_TT_flk
R_Tstar_icesnow = (R_TI_icesnow/2.0-1.0)*R_rho_c_icesnow + 1.0
d_h_ML_dt = -Q_star_flk*(R_Tstar_icesnow*(1.0+c_cbl_1)-1.0) - Q_bot_flk
d_h_ML_dt = d_h_ML_dt/tpl_rho_w_r/tpl_c_w ! Q_* and Q_b flux terms
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
d_h_ML_dt = d_h_ML_dt + flk_str_2 ! Add dC_T/dt term
flk_str_2 = I_bot_flk + (R_TI_icesnow-1.0)*I_h_flk - R_TI_icesnow*I_intm_h_D_flk
flk_str_2 = flk_str_2 + (R_TI_icesnow-2.0)*R_rho_c_icesnow*(I_h_flk-I_intm_0_h_flk)
flk_str_2 = flk_str_2/tpl_rho_w_r/tpl_c_w
d_h_ML_dt = d_h_ML_dt + flk_str_2 ! Add radiation terms
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)
flk_str_2 = flk_str_2 + C_T_p_flk*(T_wML_p_flk-T_bot_p_flk)
d_h_ML_dt = d_h_ML_dt/flk_str_2 ! dh_ML/dt = r.h.s.
END IF
!_dm
! Notice that dh_ML/dt may appear to be negative
! (e.g. due to buoyancy loss to bottom sediments and/or
! the effect of volumetric radiation heating),
! although a negative generalized buoyancy flux scale indicates
! that the equilibrium CBL depth has not yet been reached
! and convective deepening of the mixed layer should take place.
! Physically, this situation reflects an approximate character of the lake model.
! Using the self-similar temperature profile in the thermocline,
! there is always communication between the mixed layer, the thermocline
! and the lake bottom. As a result, the rate of change of the CBL depth
! is always dependent on the bottom heat flux and the radiation heating of the thermocline.
! In reality, convective mixed-layer deepening may be completely decoupled
! from the processes underneath. In order to account for this fact,
! the rate of CBL deepening is set to a small value
! if dh_ML/dt proves to be negative.
! This is "double insurance" however,
! as a negative dh_ML/dt is encountered very rarely.
!_dm
!_dbg
! IF(d_h_ML_dt.LT.0.0) THEN
! print*, 'FLake: negative d_h_ML_dt during convection, = ', d_h_ML_dt
! print*, ' d_h_ML_dt*del_time = ', MAX(d_h_ML_dt, c_small_flk)*del_time
! print*, ' u_* = ', u_star_w_flk
! print*, ' w_* = ', w_star_sfc_flk
! print*, ' h_CBL_eqi = ', conv_equil_h_scale
! print*, ' ZM scale = ', ZM_h_scale
! print*, ' h_ML_p_flk = ', h_ML_p_flk
! END IF
! print*, 'FLake: Convection, = ', d_h_ML_dt
! print*, ' Q_* = ', Q_star_flk
! print*, ' \beta*Q_* = ', flk_str_1
!_dbg
d_h_ML_dt = MAX(d_h_ML_dt, c_small_flk)
h_ML_n_flk = h_ML_p_flk + d_h_ML_dt*del_time ! Update h_ML
h_ML_n_flk = MAX(h_ML_min_flk, MIN(h_ML_n_flk, depth_w)) ! Security, limit h_ML
ELSE ! Mixing down to the lake bottom
h_ML_n_flk = depth_w
END IF
ELSE Mixing_regime ! Wind mixing
d_h_ML_dt = MAX(u_star_w_flk, u_star_min_flk) ! The surface friction velocity
ZM_h_scale = (ABS(par_Coriolis)/c_sbl_ZM_n + N_T_mean/c_sbl_ZM_i)*d_h_ML_dt**2_iintegers
ZM_h_scale = ZM_h_scale + flk_str_1/c_sbl_ZM_s
ZM_h_scale = MAX(ZM_h_scale, c_small_flk)
ZM_h_scale = d_h_ML_dt**3_iintegers/ZM_h_scale
ZM_h_scale = MAX(h_ML_min_flk, MIN(ZM_h_scale, h_ML_max_flk)) ! The ZM96 SBL depth scale
ZM_h_scale = MAX(ZM_h_scale, conv_equil_h_scale) ! Equilibrium mixed-layer depth
!_dm
! In order to avoid numerical discretization problems,
! an analytical solution to the evolution equation
! for the wind-mixed layer depth is used.
! That is, an exponential relaxation formula is applied
! over the time interval equal to the model time step.
!_dm
d_h_ML_dt = c_relax_h*d_h_ML_dt/ZM_h_scale*del_time
h_ML_n_flk = ZM_h_scale - (ZM_h_scale-h_ML_p_flk)*EXP(-d_h_ML_dt) ! Update h_ML
h_ML_n_flk = MAX(h_ML_min_flk, MIN(h_ML_n_flk, depth_w)) ! Limit h_ML
d_h_ML_dt = (h_ML_n_flk-h_ML_p_flk)/del_time ! Re-compute dh_ML/dt
IF(h_ML_n_flk.LE.h_ML_p_flk) &
d_C_T_dt = -d_C_T_dt ! Mixed-layer retreat or stationary state, dC_T/dt<0
C_T_n_flk = C_T_p_flk + d_C_T_dt*del_time ! Update C_T
C_T_n_flk = MIN(C_T_max, MAX(C_T_n_flk, C_T_min)) ! Limit C_T
d_C_T_dt = (C_T_n_flk-C_T_p_flk)/del_time ! Re-compute dC_T/dt
!_dbg
! print*, 'FLake: wind mixing: d_h_ML_dt*del_time = ', d_h_ML_dt*del_time
! print*, ' h_CBL_eqi = ', conv_equil_h_scale
! print*, ' ZM scale = ', ZM_h_scale
! print*, ' w_* = ', w_star_sfc_flk
! print*, ' u_* = ', u_star_w_flk
! print*, ' h_ML_p_flk = ', h_ML_p_flk
!_dbg
END IF Mixing_regime
! Compute the time-rate-of-change of the the bottom temperature,
! depending on the sign of dh_ML/dt
! Update the bottom temperature and the mixed-layer temperature
IF(h_ML_n_flk.LE.depth_w-h_ML_min_flk) THEN ! Mixing did not reach the bottom
IF(h_ML_n_flk.GT.h_ML_p_flk) THEN ! Mixed-layer deepening
R_H_icesnow = h_ML_p_flk/depth_w
R_rho_c_icesnow = 1.0-R_H_icesnow
R_TI_icesnow = 0.5*C_T_p_flk*R_rho_c_icesnow+C_TT_flk*(2.0*R_H_icesnow-1.0)
R_Tstar_icesnow = (0.5+C_TT_flk-C_Q_flk)/R_TI_icesnow
R_TI_icesnow = (1.0-C_T_p_flk*R_rho_c_icesnow)/R_TI_icesnow
d_T_bot_dt = (Q_w_flk-Q_bot_flk+I_w_flk-I_bot_flk)/tpl_rho_w_r/tpl_c_w
d_T_bot_dt = d_T_bot_dt - C_T_p_flk*(T_wML_p_flk-T_bot_p_flk)*d_h_ML_dt
d_T_bot_dt = d_T_bot_dt*R_Tstar_icesnow/depth_w ! Q+I fluxes and dh_ML/dt term
flk_str_2 = I_intm_h_D_flk - (1.0-C_Q_flk)*I_h_flk - C_Q_flk*I_bot_flk
flk_str_2 = flk_str_2*R_TI_icesnow/(depth_w-h_ML_p_flk)/tpl_rho_w_r/tpl_c_w
d_T_bot_dt = d_T_bot_dt + flk_str_2 ! Add radiation-flux term
flk_str_2 = (1.0-C_TT_2*R_TI_icesnow)/C_T_p_flk
flk_str_2 = flk_str_2*(T_wML_p_flk-T_bot_p_flk)*d_C_T_dt
d_T_bot_dt = d_T_bot_dt + flk_str_2 ! Add dC_T/dt term
ELSE ! Mixed-layer retreat or stationary state
d_T_bot_dt = 0.0 ! dT_bot/dt=0
END IF
T_bot_n_flk = T_bot_p_flk + d_T_bot_dt*del_time ! Update T_bot
T_bot_n_flk = MAX(T_bot_n_flk, tpl_T_f) ! Security, limit T_bot by the freezing point
flk_str_2 = (T_bot_n_flk-tpl_T_r)*flake_buoypar(T_mnw_n_flk)
IF(flk_str_2.LT.0.0) T_bot_n_flk = tpl_T_r ! Security, avoid T_r crossover
T_wML_n_flk = C_T_n_flk*(1.0-h_ML_n_flk/depth_w)
T_wML_n_flk = (T_mnw_n_flk-T_bot_n_flk*T_wML_n_flk)/(1.0-T_wML_n_flk)
T_wML_n_flk = MAX(T_wML_n_flk, tpl_T_f) ! Security, limit T_wML by the freezing point
! print*,'mark 2 T_wML_n_flk=',T_wML_n_flk
ELSE ! Mixing down to the lake bottom
h_ML_n_flk = depth_w
T_wML_n_flk = T_mnw_n_flk
T_bot_n_flk = T_mnw_n_flk
C_T_n_flk = C_T_min
! print*,'mark 3 T_wML_n_flk=',T_wML_n_flk
END IF
END IF HTC_Water
!------------------------------------------------------------------------------
! Compute the depth of the upper layer of bottom sediments
! and the temperature at that depth.
!------------------------------------------------------------------------------
Use_sediment: IF(lflk_botsed_use) THEN ! The bottom-sediment scheme is used
IF(H_B1_p_flk.GE.depth_bs-H_B1_min_flk) THEN ! No T(z) maximum (no thermal wave)
H_B1_p_flk = 0.0 ! Set H_B1_p to zero
T_B1_p_flk = T_bot_p_flk ! Set T_B1_p to the bottom temperature
END IF
flk_str_1 = 2.0*Phi_B1_pr0/(1.0-C_B1)*tpl_kappa_w/tpl_rho_w_r/tpl_c_w*del_time
h_ice_threshold = SQRT(flk_str_1) ! Threshold value of H_B1
h_ice_threshold = MIN(0.9*depth_bs, h_ice_threshold) ! Limit H_B1
flk_str_2 = C_B2/(1.0-C_B2)*(T_bs-T_B1_p_flk)/(depth_bs-H_B1_p_flk)
IF(H_B1_p_flk.LT.h_ice_threshold) THEN ! Use a truncated equation for H_B1(t)
H_B1_n_flk = SQRT(H_B1_p_flk**2_iintegers+flk_str_1) ! Advance H_B1
d_H_B1_dt = (H_B1_n_flk-H_B1_p_flk)/del_time ! Re-compute dH_B1/dt
ELSE ! Use a full equation for H_B1(t)
flk_str_1 = (Q_bot_flk+I_bot_flk)/H_B1_p_flk/tpl_rho_w_r/tpl_c_w
flk_str_1 = flk_str_1 - (1.0-C_B1)*(T_bot_n_flk-T_bot_p_flk)/del_time
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
d_H_B1_dt = flk_str_1/d_H_B1_dt
H_B1_n_flk = H_B1_p_flk + d_H_B1_dt*del_time ! Advance H_B1
END IF
d_T_B1_dt = flk_str_2*d_H_B1_dt
T_B1_n_flk = T_B1_p_flk + d_T_B1_dt*del_time ! Advance T_B1
!_dbg
! print*, 'BS module: '
! print*, ' Q_bot = ', Q_bot_flk
! print*, ' d_H_B1_dt = ', d_H_B1_dt
! print*, ' d_T_B1_dt = ', d_T_B1_dt
! print*, ' H_B1 = ', H_B1_n_flk
! print*, ' T_bot = ', T_bot_n_flk
! print*, ' T_B1 = ', T_B1_n_flk
! print*, ' T_bs = ', T_bs
!_dbg
!_nu
! Use a very simplistic procedure, where only the upper layer profile is used,
! H_B1 is always set to depth_bs, and T_B1 is always set to T_bs.
! Then, the time derivatives are zero, and the sign of the bottom heat flux depends on
! whether T_bot is smaller or greater than T_bs.
! This is, of course, an oversimplified scheme.
!_nu d_H_B1_dt = 0.0
!_nu d_T_B1_dt = 0.0
!_nu H_B1_n_flk = H_B1_p_flk + d_H_B1_dt*del_time ! Advance H_B1
!_nu T_B1_n_flk = T_B1_p_flk + d_T_B1_dt*del_time ! Advance T_B1
!_nu
l_snow_exists = H_B1_n_flk.GE.depth_bs-H_B1_min_flk & ! H_B1 reached depth_bs, or
.OR. H_B1_n_flk.LT.H_B1_min_flk & ! H_B1 decreased to zero, or
.OR.(T_bot_n_flk-T_B1_n_flk)*(T_bs-T_B1_n_flk).LE.0.0 ! there is no T(z) maximum
IF(l_snow_exists) THEN
H_B1_n_flk = depth_bs ! Set H_B1 to the depth of the thermally active layer
T_B1_n_flk = T_bs ! Set T_B1 to the climatological temperature
END IF
! changed by Shaobo Zhang
!ELSE Use_sediment ! The bottom-sediment scheme is not used
!
! H_B1_n_flk = rflk_depth_bs_ref ! H_B1 is set to a reference value
! T_B1_n_flk = tpl_T_r ! T_B1 is set to the temperature of maximum density
ELSE Use_sediment ! The scheme written by Shaobo Zhang for the deeper layer of a deep lake is used
if ( abs(T_bot_p_flk-T_bot_2_in) .lt. 0.01 ) then
depth_bs = depthbs
depth_w = depthw
T_bot_2_out = T_bot_2_in
else
CT = 0.65
CTT = 11.0/18.0 * CT - 7.0/45.0
CQ = 2.0 * CTT / CT
! compute d_h_D_dt
flk_str_1 = ( Q_bot_flk * CQ + I_bot_flk - I_intm_D_H_flk/depth_bs ) /tpl_rho_w_r/tpl_c_w
flk_str_1 = flk_str_1 - depth_bs * (0.5-CTT) * (T_bot_n_flk-T_bot_p_flk)/del_time
flk_str_1 = flk_str_1 - CTT/CT*( (Q_bot_flk+I_bot_flk-I_HH_flk)/tpl_rho_w_r/tpl_c_w - &
depth_bs * ( 1.0 - CT ) * (T_bot_n_flk-T_bot_p_flk)/del_time )
flk_str_2 = CTT * (T_bot_p_flk-T_bot_2_in)
if(abs(flk_str_2)<0.01) then
d_h_D_dt = 0.0
else
d_h_D_dt = flk_str_1/flk_str_2
endif
! compute d_T_H_dt
flk_str_1 = (Q_bot_flk+I_bot_flk-I_HH_flk)/tpl_rho_w_r/tpl_c_w
flk_str_1 = flk_str_1 - depth_bs*(1.0-CT)*(T_bot_n_flk-T_bot_p_flk)/del_time
flk_str_1 = flk_str_1 - CT * (T_bot_p_flk-T_bot_2_in) * d_h_D_dt
flk_str_2 = CT * depth_bs
d_T_H_dt = flk_str_1/flk_str_2
! update T_bot_2_out
T_bot_2_out = T_bot_2_in + d_T_H_dt * del_time
if ( (T_bot_2_out-tpl_T_r)*(T_bot_n_flk-tpl_T_r) .le. 0.0 ) then
T_bot_2_out = tpl_T_r
else if ( (T_bot_n_flk-tpl_T_r)/(T_bot_2_out-tpl_T_r) .le. 1.0 ) then
T_bot_2_out = T_bot_n_flk
end if
! update depth_w
flk_str_2 = depth_w + depth_bs ! The total depth of the deep lake
depth_w = depth_w - d_h_D_dt * del_time ! Advance depth_bs
! depth_w = max (lake_depth_max, min(depth_w, 0.4*flk_str_2)) ! Limit depth_bs
depth_w = max (max(lake_depth_max,0.3*flk_str_2), &
min(depth_w, min(0.4*flk_str_2, 80.0))) ! Limit depth_bs
depth_bs = flk_str_2 - depth_w ! Update depth_w
end if
END IF Use_sediment
!------------------------------------------------------------------------------
! Impose additional constraints.
!------------------------------------------------------------------------------
! In case of unstable stratification, force mixing down to the bottom
flk_str_2 = (T_wML_n_flk-T_bot_n_flk)*flake_buoypar(T_mnw_n_flk)
IF(flk_str_2.LT.0.0) THEN
!_dbg
! print*, 'FLake: inverse (unstable) stratification !!! '
! print*, ' Mixing down to the bottom is forced.'
! print*, ' T_wML_p, T_wML_n ', T_wML_p_flk-tpl_T_f, T_wML_n_flk-tpl_T_f
! print*, ' T_mnw_p, T_mnw_n ', T_mnw_p_flk-tpl_T_f, T_mnw_n_flk-tpl_T_f
! print*, ' T_bot_p, T_bot_n ', T_bot_p_flk-tpl_T_f, T_bot_n_flk-tpl_T_f
! print*, ' h_ML_p, h_ML_n ', h_ML_p_flk, h_ML_n_flk
! print*, ' C_T_p, C_T_n ', C_T_p_flk, C_T_n_flk
!_dbg
h_ML_n_flk = depth_w
T_wML_n_flk = T_mnw_n_flk
T_bot_n_flk = T_mnw_n_flk
C_T_n_flk = C_T_min
! print*,'mark 4 T_wML_n_flk=',T_wML_n_flk
END IF
!------------------------------------------------------------------------------
! Update the surface temperature.
!------------------------------------------------------------------------------
IF(h_snow_n_flk.GE.h_Snow_min_flk) THEN
T_sfc_n = T_snow_n_flk ! Snow exists, use the snow temperature
ELSE IF(h_ice_n_flk.GE.h_Ice_min_flk) THEN
T_sfc_n = T_ice_n_flk ! Ice exists but there is no snow, use the ice temperature
ELSE
T_sfc_n = T_wML_n_flk ! No ice-snow cover, use the mixed-layer temperature
END IF
if(T_sfc_n.lt.200.0 .or. T_sfc_n.gt.350.0) then
print*,'h_snow_n_flk=',h_snow_n_flk,' h_Snow_min_flk=',h_Snow_min_flk
print*,'h_ice_n_flk=',h_ice_n_flk, ' h_Ice_min_flk= ',h_Ice_min_flk
print*,'T_wML_n_flk=',T_wML_n_flk
print*,'T_sfc_n=',T_sfc_n
endif
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END SUBROUTINE flake_main
!==============================================================================
!==============================================================================
! include 'flake_buoypar.incf'
!------------------------------------------------------------------------------
REAL (KIND = ireals) FUNCTION flake_buoypar (T_water)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the buoyancy parameter,
! using a quadratic equation of state for the fresh-water.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "flake".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
USE flake_parameters , ONLY : &
tpl_grav , & ! Acceleration due to gravity [m s^{-2}]
tpl_T_r , & ! Temperature of maximum density of fresh water [K]
tpl_a_T ! Constant in the fresh-water equation of state [K^{-2}]
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T_water ! Water temperature [K]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Buoyancy parameter [m s^{-2} K^{-1}]
flake_buoypar = tpl_grav*tpl_a_T*(T_water-tpl_T_r)
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION flake_buoypar
!==============================================================================
!==============================================================================
! include 'flake_snowdensity.incf'
!------------------------------------------------------------------------------
REAL (KIND = ireals) FUNCTION flake_snowdensity (h_snow)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the snow density,
! using an empirical approximation from Heise et al. (2003).
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "flake".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
USE flake_parameters , ONLY : &
tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
tpl_rho_S_min , & ! Minimum snow density [kg m^{-3}]
tpl_rho_S_max , & ! Maximum snow density [kg m^{-3}]
tpl_Gamma_rho_S , & ! Empirical parameter [kg m^{-4}] in the expression for the snow density
c_small_flk ! A small number
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
h_snow ! Snow thickness [m]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Snow density [kg m^{-3}]
! Security. Ensure that the expression in () does not become negative at a very large h_snow.
flake_snowdensity = MAX( c_small_flk, (1.0 - h_snow*tpl_Gamma_rho_S/tpl_rho_w_r) )
flake_snowdensity = MIN( tpl_rho_S_max, tpl_rho_S_min/flake_snowdensity )
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION flake_snowdensity
!==============================================================================
!==============================================================================
! include 'flake_snowheatconduct.incf'
!------------------------------------------------------------------------------
REAL (KIND = ireals) FUNCTION flake_snowheatconduct (h_snow)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the snow heat conductivity,
! using an empirical approximation from Heise et al. (2003).
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "flake".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
USE flake_parameters , ONLY : &
tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
tpl_kappa_S_min , & ! Minimum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
tpl_kappa_S_max , & ! Maximum molecular heat conductivity of snow [J m^{-1} s^{-1} K^{-1}]
tpl_Gamma_kappa_S ! Empirical parameter [J m^{-2} s^{-1} K^{-1}] in the expression for kappa_S
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
h_snow ! Snow thickness [m]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Snow heat conductivity [J m^{-1} s^{-1} K^{-1} = kg m s^{-3} K^{-1}]
flake_snowheatconduct = flake_snowdensity( h_snow ) ! Compute snow density
flake_snowheatconduct = MIN( tpl_kappa_S_max, tpl_kappa_S_min &
+ h_snow*tpl_Gamma_kappa_S*flake_snowheatconduct/tpl_rho_w_r )
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION flake_snowheatconduct
!==============================================================================
END MODULE flake
!------------------------------------------------------------------------------
MODULE SfcFlx
!------------------------------------------------------------------------------
!
! Description:
!
! The main program unit of
! the atmospheric surface-layer parameterization scheme "SfcFlx".
! "SfcFlx" is used to compute fluxes
! of momentum and of sensible and latent heat over lakes.
! The surface-layer scheme developed by Mironov (1991) was used as the starting point.
! It was modified and further developed to incorporate new results as to
! the roughness lenghts for scalar quantities,
! heat and mass transfer in free convection,
! and the effect of limited fetch on the momentum transfer.
! Apart from the momentum flux and sensible and latent heat fluxes,
! the long-wave radiation flux from the water surface and
! the long-wave radiation flux from the atmosphere can also be computed.
! The atmospheric long-wave radiation flux is computed with simple empirical formulae,
! where the atmospheric emissivity is taken to be dependent on
! the water vapour pressure and cloud fraction.
!
! A description of SfcFlx is available from the author.
! Dmitrii Mironov
! German Weather Service, Kaiserleistr. 29/35, D-63067 Offenbach am Main, Germany.
! dmitrii.mironov@dwd.de
!
! Lines embraced with "!_tmp" contain temporary parts of the code.
! Lines embraced/marked with "!_dev" may be replaced
! as improved parameterizations are developed and tested.
! Lines embraced/marked with "!_dm" are DM's comments
! that may be helpful to a user.
! Lines embraced/marked with "!_dbg" are used
! for debugging purposes only.
!
USE data_parameters , ONLY : &
ireals , & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
USE flake_parameters , ONLY : &
tpl_grav , & ! Acceleration due to gravity [m s^{-2}]
tpl_T_f , & ! Fresh water freezing point [K]
tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
tpl_c_w , & ! Specific heat of water [J kg^{-1} K^{-1}]
tpl_L_f , & ! Latent heat of fusion [J kg^{-1}]
h_Ice_min_flk ! Minimum ice thickness [m]
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Dimensionless constants in the Monin-Obukhov surface-layer
! similarity relations and in the expressions for the roughness lengths.
REAL (KIND = kind_phys), PARAMETER :: &
c_Karman = 0.40 , & ! The von Karman constant
! Pr_neutral = 1.0 , & ! Turbulent Prandtl number at neutral static stability
Pr_neutral = 0.9 , & ! Turbulent Prandtl number at neutral static stability
Sc_neutral = 1.0 , & ! Turbulent Schmidt number at neutral static stability
c_MO_u_stab = 5.0 , & ! Constant of the MO theory (wind, stable stratification)
c_MO_t_stab = 5.0 , & ! Constant of the MO theory (temperature, stable stratification)
c_MO_q_stab = 5.0 , & ! Constant of the MO theory (humidity, stable stratification)
c_MO_u_conv = 15.0 , & ! Constant of the MO theory (wind, convection)
c_MO_t_conv = 15.0 , & ! Constant of the MO theory (temperature, convection)
c_MO_q_conv = 15.0 , & ! Constant of the MO theory (humidity, convection)
c_MO_u_exp = 0.25 , & ! Constant of the MO theory (wind, exponent)
c_MO_t_exp = 0.5 , & ! Constant of the MO theory (temperature, exponent)
c_MO_q_exp = 0.5 , & ! Constant of the MO theory (humidity, exponent)
z0u_ice_rough = 1.0E-03 , & ! Aerodynamic roughness of the ice surface [m] (rough flow)
c_z0u_smooth = 0.1 , & ! Constant in the expression for z0u (smooth flow)
c_z0u_rough = 1.23E-02 , & ! The Charnock constant in the expression for z0u (rough flow)
c_z0u_rough_L = 1.00E-01 , & ! An increased Charnock constant (used as the upper limit)
c_z0u_ftch_f = 0.70 , & ! Factor in the expression for fetch-dependent Charnock parameter
c_z0u_ftch_ex = 0.3333333 , & ! Exponent in the expression for fetch-dependent Charnock parameter
c_z0t_rough_1 = 4.0 , & ! Constant in the expression for z0t (factor)
c_z0t_rough_2 = 3.2 , & ! Constant in the expression for z0t (factor)
c_z0t_rough_3 = 0.5 , & ! Constant in the expression for z0t (exponent)
c_z0q_rough_1 = 4.0 , & ! Constant in the expression for z0q (factor)
c_z0q_rough_2 = 4.2 , & ! Constant in the expression for z0q (factor)
c_z0q_rough_3 = 0.5 , & ! Constant in the expression for z0q (exponent)
c_z0t_ice_b0s = 1.250 , & ! Constant in the expression for z0t over ice
c_z0t_ice_b0t = 0.149 , & ! Constant in the expression for z0t over ice
c_z0t_ice_b1t = -0.550 , & ! Constant in the expression for z0t over ice
c_z0t_ice_b0r = 0.317 , & ! Constant in the expression for z0t over ice
c_z0t_ice_b1r = -0.565 , & ! Constant in the expression for z0t over ice
c_z0t_ice_b2r = -0.183 , & ! Constant in the expression for z0t over ice
c_z0q_ice_b0s = 1.610 , & ! Constant in the expression for z0q over ice
c_z0q_ice_b0t = 0.351 , & ! Constant in the expression for z0q over ice
c_z0q_ice_b1t = -0.628 , & ! Constant in the expression for z0q over ice
c_z0q_ice_b0r = 0.396 , & ! Constant in the expression for z0q over ice
c_z0q_ice_b1r = -0.512 , & ! Constant in the expression for z0q over ice
c_z0q_ice_b2r = -0.180 , & ! Constant in the expression for z0q over ice
Re_z0s_ice_t = 2.5 , & ! Threshold value of the surface Reynolds number
! used to compute z0t and z0q over ice (Andreas 2002)
Re_z0u_thresh = 0.1 ! Threshold value of the roughness Reynolds number
! [value from Zilitinkevich, Grachev, and Fairall (200),
! currently not used]
! Dimensionless constants
REAL (KIND = kind_phys), PARAMETER :: &
c_free_conv = 0.14 ! Constant in the expressions for fluxes in free convection
! Dimensionless constants
REAL (KIND = kind_phys), PARAMETER :: &
c_lwrad_emis = 0.99 ! Surface emissivity with respect to the long-wave radiation
! Thermodynamic parameters
REAL (KIND = kind_phys), PARAMETER :: &
tpsf_C_StefBoltz = 5.67E-08 , & ! The Stefan-Boltzmann constant [W m^{-2} K^{-4}]
tpsf_R_dryair = 2.8705E+02 , & ! Gas constant for dry air [J kg^{-1} K^{-1}]
tpsf_R_watvap = 4.6151E+02 , & ! Gas constant for water vapour [J kg^{-1} K^{-1}]
tpsf_c_a_p = 1.005E+03 , & ! Specific heat of air at constant pressure [J kg^{-1} K^{-1}]
tpsf_L_evap = 2.501E+06 , & ! Specific heat of evaporation [J kg^{-1}]
tpsf_nu_u_a = 1.50E-05 , & ! Kinematic molecular viscosity of air [m^{2} s^{-1}]
tpsf_kappa_t_a = 2.20E-05 , & ! Molecular temperature conductivity of air [m^{2} s^{-1}]
tpsf_kappa_q_a = 2.40E-05 ! Molecular diffusivity of air for water vapour [m^{2} s^{-1}]
! Derived thermodynamic parameters
REAL (KIND = kind_phys), PARAMETER :: &
tpsf_Rd_o_Rv = tpsf_R_dryair/tpsf_R_watvap , & ! Ratio of gas constants (Rd/Rv)
tpsf_alpha_q = (1.0-tpsf_Rd_o_Rv)/tpsf_Rd_o_Rv ! Diemsnionless ratio
! Thermodynamic parameters
REAL (KIND = kind_phys), PARAMETER :: &
P_a_ref = 1.0E+05 ! Reference pressure [N m^{-2} = kg m^{-1} s^{-2}]
! The variables declared below
! are accessible to all program units of the MODULE "SfcFlx"
! and to the driving routines that use "SfcFlx".
! These are basically the quantities computed by SfcFlx.
! Apart from these quantities, there a few local scalars
! used by SfcFlx routines mainly for security reasons.
! All variables declared below have a suffix "sf".
! SfcFlx variables of type REAL
! Roughness lengths
REAL (KIND = kind_phys) :: &
z0u_sf , & ! Roughness length with respect to wind velocity [m]
z0t_sf , & ! Roughness length with respect to potential temperature [m]
z0q_sf ! Roughness length with respect to specific humidity [m]
! Fluxes in the surface air layer
REAL (KIND = kind_phys) :: &
u_star_a_sf , & ! Friction velocity [m s^{-1}]
Q_mom_a_sf , & ! Momentum flux [N m^{-2}]
Q_sens_a_sf , & ! Sensible heat flux [W m^{-2}]
Q_lat_a_sf , & ! Laten heat flux [W m^{-2}]
Q_watvap_a_sf ! Flux of water vapout [kg m^{-2} s^{-1}]
! Security constants
REAL (KIND = kind_phys), PARAMETER :: &
u_wind_min_sf = 1.0E-02 , & ! Minimum wind speed [m s^{-1}]
u_star_min_sf = 1.0E-04 , & ! Minimum value of friction velocity [m s^{-1}]
c_accur_sf = 1.0E-07 , & ! A small number (accuracy)
c_small_sf = 1.0E-04 ! A small number (used to compute fluxes)
! Useful constants
REAL (KIND = kind_phys), PARAMETER :: &
num_1o3_sf = 1.0/3.0 ! 1/3
!==============================================================================
! Procedures
!==============================================================================
CONTAINS
!==============================================================================
! The codes of the SfcFlx procedures are stored in separate "*.incf" files
! and are included below.
!------------------------------------------------------------------------------
!==============================================================================
! include 'SfcFlx_lwradatm.incf'
!------------------------------------------------------------------------------
REAL (KIND = kind_phys) FUNCTION SfcFlx_lwradatm (T_a, e_a, cl_tot, cl_low)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the long-wave radiation flux from the atmosphere
! as function of air temperature, water vapour pressure and cloud fraction.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T_a , & ! Air temperature [K]
e_a , & ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
cl_tot , & ! Total cloud cover [0,1]
cl_low ! Lowe-level cloud cover [0,1]
! Local parameters
! Coefficients in the empirical formulation
! developed at the Main Geophysical Observatory (MGO), St. Petersburg, Russia.
REAL (KIND = kind_phys), PARAMETER :: &
c_lmMGO_1 = 43.057924 , & ! Empirical coefficient
c_lmMGO_2 = 540.795 ! Empirical coefficient
! Temperature-dependent cloud-correction coefficients in the MGO formula
INTEGER (KIND = iintegers), PARAMETER :: &
nband_coef = 6_iintegers ! Number of temperature bands
REAL (KIND = kind_phys), PARAMETER, DIMENSION (nband_coef) :: &
corr_cl_tot = (/0.70, 0.45, 0.32, &
0.23, 0.18, 0.13/) , & ! Total clouds
corr_cl_low = (/0.76, 0.49, 0.35, &
0.26, 0.20, 0.15/) , & ! Low-level clouds
corr_cl_midhigh = (/0.46, 0.30, 0.21, &
0.15, 0.12, 0.09/) ! Mid- and high-level clouds
REAL (KIND = kind_phys), PARAMETER :: &
T_low = 253.15 , & ! Low-limit temperature in the interpolation formula [K]
del_T = 10.0 ! Temperature step in the interpolation formula [K]
! Coefficients in the empirical water-vapour correction function
! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001).
REAL (KIND = kind_phys), PARAMETER :: &
c_watvap_corr_min = 0.6100 , & ! Empirical coefficient (minimum value of the correction function)
c_watvap_corr_max = 0.7320 , & ! Empirical coefficient (maximum value of the correction function)
c_watvap_corr_e = 0.0050 ! Empirical coefficient [(N m^{-2})^{-1/2}]
! Local variables of type INTEGER
INTEGER (KIND = iintegers) :: &
i ! Loop index
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
c_cl_tot_corr , & ! The MGO cloud correction coefficient, total clouds
c_cl_low_corr , & ! The MGO cloud correction coefficient, low-level clouds
c_cl_midhigh_corr , & ! The MGO cloud correction coefficient, mid- and high-level clouds
T_corr , & ! Temperature used to compute the MGO cloud correction [K]
f_wvpres_corr , & ! Correction function with respect to water vapour
f_cloud_corr ! Correction function with respect to cloudiness
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Water-vapour correction function
f_wvpres_corr = c_watvap_corr_min + c_watvap_corr_e*SQRT(e_a)
f_wvpres_corr = MIN(f_wvpres_corr, c_watvap_corr_max)
! Cloud-correction coefficients using the MGO formulation with linear interpolation
IF(T_a.LT.T_low) THEN
c_cl_tot_corr = corr_cl_tot(1)
c_cl_low_corr = corr_cl_low(1)
c_cl_midhigh_corr = corr_cl_midhigh(1)
ELSE IF(T_a.GE.T_low+(nband_coef-1_iintegers)*del_T) THEN
c_cl_tot_corr = corr_cl_tot(nband_coef)
c_cl_low_corr = corr_cl_low(nband_coef)
c_cl_midhigh_corr = corr_cl_midhigh(nband_coef)
ELSE
T_corr = T_low
DO i=1, nband_coef-1
IF(T_a.GE.T_corr.AND.T_a.LT.T_corr+del_T) THEN
c_cl_tot_corr = (T_a-T_corr)/del_T
c_cl_low_corr = corr_cl_low(i) + (corr_cl_low(i+1)-corr_cl_low(i))*c_cl_tot_corr
c_cl_midhigh_corr = corr_cl_midhigh(i) + (corr_cl_midhigh(i+1)-corr_cl_midhigh(i))*c_cl_tot_corr
c_cl_tot_corr = corr_cl_tot(i) + (corr_cl_tot(i+1)-corr_cl_tot(i))*c_cl_tot_corr
END IF
T_corr = T_corr + del_T
END DO
END IF
! Cloud correction function
IF(cl_low.LT.0.0) THEN ! Total cloud cover only
f_cloud_corr = 1.0 + c_cl_tot_corr*cl_tot*cl_tot
ELSE ! Total and low-level cloud cover
f_cloud_corr = (1.0 + c_cl_low_corr*cl_low*cl_low) &
* (1.0 + c_cl_midhigh_corr*(cl_tot*cl_tot-cl_low*cl_low))
END IF
! Long-wave radiation flux [W m^{-2}]
! The MGO formulation
!_nu The MGO formulation
!_nu SfcFlx_lwradatm = -SfcFlx_lwradatm*c_lwrad_emis &
!_nu * (c_lmMGO_1*SQRT(tpsf_C_StefBoltz*T_a**4_iintegers)-c_lmMGO_2)
!_nu
! "Conventional" formulation
! (see Fung et al. 1984, Zapadka and Wozniak 2000, Zapadka et al. 2001)
SfcFlx_lwradatm = -c_lwrad_emis*tpsf_C_StefBoltz*T_a**4_iintegers &
* f_wvpres_corr*f_cloud_corr
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_lwradatm
!==============================================================================
!==============================================================================
! include 'SfcFlx_lwradwsfc.incf'
!------------------------------------------------------------------------------
REAL (KIND = kind_phys) FUNCTION SfcFlx_lwradwsfc (T)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the surface long-wave radiation flux
! as function of temperature.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T ! Temperature [K]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Long-wave radiation flux [W m^{-2}]
SfcFlx_lwradwsfc = c_lwrad_emis*tpsf_C_StefBoltz*T**4_iintegers
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_lwradwsfc
!==============================================================================
!==============================================================================
! include 'SfcFlx_momsenlat.incf'
!------------------------------------------------------------------------------
SUBROUTINE SfcFlx_momsenlat ( height_u, height_tq, fetch, &
U_a, T_a, q_a, T_s, P_a, h_ice, &
Q_momentum, Q_sensible, Q_latent, Q_watvap,&
q_s,rho_a )
!------------------------------------------------------------------------------
!
! Description:
!
! The SfcFlx routine
! where fluxes of momentum and of sensible and latent heat
! at the air-water or air-ice (air-snow) interface are computed.
!
! Lines embraced with "!_tmp" contain temporary parts of the code.
! Lines embraced/marked with "!_dev" may be replaced
! as improved parameterizations are developed and tested.
! Lines embraced/marked with "!_dm" are DM's comments
! that may be helpful to a user.
! Lines embraced/marked with "!_dbg" are used
! for debugging purposes only.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (procedure arguments)
REAL (KIND = kind_phys), INTENT(IN) :: &
height_u , & ! Height where wind is measured [m]
height_tq , & ! Height where temperature and humidity are measured [m]
fetch , & ! Typical wind fetch [m]
U_a , & ! Wind speed [m s^{-1}]
T_a , & ! Air temperature [K]
q_a , & ! Air specific humidity [-]
T_s , & ! Surface temperature (water, ice or snow) [K]
P_a , & ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}]
h_ice ! Ice thickness [m]
! Output (procedure arguments)
REAL (KIND = kind_phys), INTENT(OUT) :: &
Q_momentum , & ! Momentum flux [N m^{-2}]
Q_sensible , & ! Sensible heat flux [W m^{-2}]
Q_latent , & ! Laten heat flux [W m^{-2}]
Q_watvap ! Flux of water vapout [kg m^{-2} s^{-1}]
! Local parameters of type INTEGER
INTEGER (KIND = iintegers), PARAMETER :: &
n_iter_max = 24 ! Maximum number of iterations
! Local variables of type LOGICAL
LOGICAL :: &
l_conv_visc , & ! Switch, TRUE = viscous free convection, the Nu=C Ra^(1/3) law is used
l_conv_cbl ! Switch, TRUE = CBL scale convective structures define surface fluxes
! Local variables of type INTEGER
INTEGER (KIND = iintegers) :: &
i , & ! Loop index
n_iter ! Number of iterations performed
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
rho_a , & ! Air density [kg m^{-3}]
wvpres_s , & ! Saturation water vapour pressure at T=T_s [N m^{-2}]
q_s ! Saturation specific humidity at T=T_s [-]
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
Q_mom_tur , & ! Turbulent momentum flux [N m^{-2}]
Q_sen_tur , & ! Turbulent sensible heat flux [W m^{-2}]
Q_lat_tur , & ! Turbulent laten heat flux [W m^{-2}]
Q_mom_mol , & ! Molecular momentum flux [N m^{-2}]
Q_sen_mol , & ! Molecular sensible heat flux [W m^{-2}]
Q_lat_mol , & ! Molecular laten heat flux [W m^{-2}]
Q_mom_con , & ! Momentum flux in free convection [N m^{-2}]
Q_sen_con , & ! Sensible heat flux in free convection [W m^{-2}]
Q_lat_con ! Laten heat flux in free convection [W m^{-2}]
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
par_conv_visc , & ! Viscous convection stability parameter
par_conv_cbl , & ! CBL convection stability parameter
c_z0u_fetch , & ! Fetch-dependent Charnock parameter
U_a_thresh , & ! Threshld value of the wind speed [m s^{-1}]
u_star_thresh , & ! Threshld value of friction velocity [m s^{-1}]
u_star_previter , & ! Friction velocity from previous iteration [m s^{-1}]
u_star_n , & ! Friction velocity at neutral stratification [m s^{-1}]
u_star_st , & ! Friction velocity with due regard for stratification [m s^{-1}]
ZoL , & ! The z/L ratio, z=height_u
Ri , & ! Gradient Richardson number
Ri_cr , & ! Critical value of Ri
R_z , & ! Ratio of "height_tq" to "height_u"
Fun , & ! A function of generic variable "x"
Fun_prime , & ! Derivative of "Fun" with respect to "x"
Delta , & ! Relative error
psi_u , & ! The MO stability function for wind profile
psi_t , & ! The MO stability function for temperature profile
psi_q ! The MO stability function for specific humidity profile
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
!write(*,*) 'Inside Flake Flux Module'
!write(*,*) height_u,height_tq,fetch,U_a,T_a,q_a,T_s,P_a,h_ice
!_dm All fluxes are positive when directed upwards.
!------------------------------------------------------------------------------
! Compute saturation specific humidity and the air density at T=T_s
!------------------------------------------------------------------------------
wvpres_s = SfcFlx_satwvpres(T_s, h_ice) ! Saturation water vapour pressure at T=T_s
q_s = SfcFlx_spechum (wvpres_s, P_a) ! Saturation specific humidity at T=T_s
rho_a = SfcFlx_rhoair(T_s, q_s, P_a) ! Air density at T_s and q_s (surface values)
!------------------------------------------------------------------------------
! Compute molecular fluxes of momentum and of sensible and latent heat
!------------------------------------------------------------------------------
!_dm The fluxes are in kinematic units
Q_mom_mol = -tpsf_nu_u_a*U_a/height_u
Q_sen_mol = -tpsf_kappa_t_a*(T_a-T_s)/height_tq
Q_lat_mol = -tpsf_kappa_q_a*(q_a-q_s)/height_tq
!------------------------------------------------------------------------------
! Compute fluxes in free convection
!------------------------------------------------------------------------------
par_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)
IF(par_conv_visc.GT.0.0) THEN ! Viscous convection takes place
l_conv_visc = .TRUE.
par_conv_visc = (par_conv_visc*tpl_grav/tpsf_nu_u_a)**num_1o3_sf
Q_sen_con = c_free_conv*SQRT(tpsf_kappa_t_a)*par_conv_visc
Q_sen_con = Q_sen_con*(T_s-T_a)
Q_lat_con = c_free_conv*SQRT(tpsf_kappa_q_a)*par_conv_visc
Q_lat_con = Q_lat_con*(q_s-q_a)
ELSE ! No viscous convection, set fluxes to zero
l_conv_visc = .FALSE.
Q_sen_con = 0.0
Q_lat_con = 0.0
END IF
Q_mom_con = 0.0 ! Momentum flux in free (viscous or CBL-scale) convection is zero
!------------------------------------------------------------------------------
! Compute turbulent fluxes
!------------------------------------------------------------------------------
R_z = height_tq/height_u ! Ratio of "height_tq" to "height_u"
Ri_cr = c_MO_t_stab/c_MO_u_stab**2_iintegers*R_z ! Critical Ri
Ri = tpl_grav*((T_a-T_s)/T_s+tpsf_alpha_q*(q_a-q_s))/MAX(U_a,u_wind_min_sf)**2_iintegers
Ri = Ri*height_u/Pr_neutral ! Gradient Richardson number
Turb_Fluxes: IF(U_a.LT.u_wind_min_sf.OR.Ri.GT.Ri_cr-c_small_sf) THEN ! Low wind or Ri>Ri_cr
u_star_st = 0.0 ! Set turbulent fluxes to zero
Q_mom_tur = 0.0
Q_sen_tur = 0.0
Q_lat_tur = 0.0
ELSE Turb_Fluxes ! Compute turbulent fluxes using MO similarity
! Compute z/L, where z=height_u
IF(Ri.GE.0.0) THEN ! Stable stratification
ZoL = SQRT(1.0-4.0*(c_MO_u_stab-R_z*c_MO_t_stab)*Ri)
ZoL = ZoL - 1.0 + 2.0*c_MO_u_stab*Ri
ZoL = ZoL/2.0/c_MO_u_stab/c_MO_u_stab/(Ri_cr-Ri)
ELSE ! Convection
n_iter = 0_iintegers
Delta = 1.0 ! Set initial error to a large value (as compared to the accuracy)
u_star_previter = Ri*MAX(1.0, SQRT(R_z*c_MO_t_conv/c_MO_u_conv)) ! Initial guess for ZoL
DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
Fun = u_star_previter**2_iintegers*(c_MO_u_conv*u_star_previter-1.0) &
+ Ri**2_iintegers*(1.0-R_z*c_MO_t_conv*u_star_previter)
Fun_prime = 3.0*c_MO_u_conv*u_star_previter**2_iintegers &
- 2.0*u_star_previter - R_z*c_MO_t_conv*Ri**2_iintegers
ZoL = u_star_previter - Fun/Fun_prime
Delta = ABS(ZoL-u_star_previter)/MAX(c_accur_sf, ABS(ZoL+u_star_previter))
u_star_previter = ZoL
n_iter = n_iter + 1_iintegers
END DO
!_dbg
! IF(n_iter.GE.n_iter_max-1_iintegers) &
! print*(*,*) 'ZoL: Max No. iters. exceeded (n_iter = ', n_iter, ')!'
!_dbg
END IF
! Compute fetch-dependent Charnock parameter, use "u_star_min_sf"
CALL SfcFlx_roughness (fetch, U_a, u_star_min_sf, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
! Threshold value of wind speed
u_star_st = u_star_thresh
CALL SfcFlx_roughness (fetch, U_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
IF(ZoL.GT.0.0) THEN ! MO function in stable stratification
psi_u = c_MO_u_stab*ZoL*(1.0-MIN(z0u_sf/height_u, 1.0))
ELSE ! MO function in convection
psi_t = (1.0-c_MO_u_conv*ZoL)**c_MO_u_exp
psi_q = (1.0-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.0))**c_MO_u_exp
psi_u = 2.0*(ATAN(psi_t)-ATAN(psi_q)) &
+ 2.0*LOG((1.0+psi_q)/(1.0+psi_t)) &
+ LOG((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
END IF
U_a_thresh = u_star_thresh/c_Karman*(LOG(height_u/z0u_sf)+psi_u)
! Compute friction velocity
n_iter = 0_iintegers
Delta = 1.0 ! Set initial error to a large value (as compared to the accuracy)
u_star_previter = u_star_thresh ! Initial guess for friction velocity
IF(U_a.LE.U_a_thresh) THEN ! Smooth surface
DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
CALL SfcFlx_roughness (fetch, U_a, MIN(u_star_thresh, u_star_previter), h_ice, &
c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
IF(ZoL.GE.0.0) THEN ! Stable stratification
psi_u = c_MO_u_stab*ZoL*(1.0-MIN(z0u_sf/height_u, 1.0))
Fun = LOG(height_u/z0u_sf) + psi_u
Fun_prime = (Fun + 1.0 + c_MO_u_stab*ZoL*MIN(z0u_sf/height_u, 1.0))/c_Karman
Fun = Fun*u_star_previter/c_Karman - U_a
ELSE ! Convection
psi_t = (1.0-c_MO_u_conv*ZoL)**c_MO_u_exp
psi_q = (1.0-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.0))**c_MO_u_exp
psi_u = 2.0*(ATAN(psi_t)-ATAN(psi_q)) &
+ 2.0*LOG((1.0+psi_q)/(1.0+psi_t)) &
+ LOG((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
Fun = LOG(height_u/z0u_sf) + psi_u
Fun_prime = (Fun + 1.0/psi_q)/c_Karman
Fun = Fun*u_star_previter/c_Karman - U_a
END IF
u_star_st = u_star_previter - Fun/Fun_prime
Delta = ABS((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
u_star_previter = u_star_st
n_iter = n_iter + 1_iintegers
END DO
ELSE ! Rough surface
DO WHILE (Delta.GT.c_accur_sf.AND.n_iter.LT.n_iter_max)
CALL SfcFlx_roughness (fetch, U_a, MAX(u_star_thresh, u_star_previter), h_ice, &
c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
IF(ZoL.GE.0.0) THEN ! Stable stratification
psi_u = c_MO_u_stab*ZoL*(1.0-MIN(z0u_sf/height_u, 1.0))
Fun = LOG(height_u/z0u_sf) + psi_u
Fun_prime = (Fun - 2.0 - 2.0*c_MO_u_stab*ZoL*MIN(z0u_sf/height_u, 1.0))/c_Karman
Fun = Fun*u_star_previter/c_Karman - U_a
ELSE ! Convection
psi_t = (1.0-c_MO_u_conv*ZoL)**c_MO_u_exp
psi_q = (1.0-c_MO_u_conv*ZoL*MIN(z0u_sf/height_u, 1.0))**c_MO_u_exp
psi_u = 2.0*(ATAN(psi_t)-ATAN(psi_q)) &
+ 2.0*LOG((1.0+psi_q)/(1.0+psi_t)) &
+ LOG((1.0+psi_q*psi_q)/(1.0+psi_t*psi_t))
Fun = LOG(height_u/z0u_sf) + psi_u
Fun_prime = (Fun - 2.0/psi_q)/c_Karman
Fun = Fun*u_star_previter/c_Karman - U_a
END IF
IF(h_ice.GE.h_Ice_min_flk) THEN ! No iteration is required for rough flow over ice
u_star_st = c_Karman*U_a/MAX(c_small_sf, LOG(height_u/z0u_sf)+psi_u)
u_star_previter = u_star_st
ELSE ! Iterate in case of open water
u_star_st = u_star_previter - Fun/Fun_prime
END IF
Delta = ABS((u_star_st-u_star_previter)/(u_star_st+u_star_previter))
u_star_previter = u_star_st
n_iter = n_iter + 1_iintegers
END DO
END IF
!_dbg
! print*(*,*) 'MO stab. func. psi_u = ', psi_u, ' n_iter = ', n_iter
! print*(*,*) ' Wind speed = ', U_a, ' u_* = ', u_star_st
! print*(*,*) ' Fun = ', Fun
!_dbg
!_dbg
! IF(n_iter.GE.n_iter_max-1_iintegers) &
! print*(*,*) 'u_*: Max No. iters. exceeded (n_iter = ', n_iter, ')!'
!_dbg
! Momentum flux
Q_mom_tur = -u_star_st*u_star_st
! Temperature and specific humidity fluxes
CALL SfcFlx_roughness (fetch, U_a, u_star_st, h_ice, c_z0u_fetch, u_star_thresh, z0u_sf, z0t_sf, z0q_sf)
IF(ZoL.GE.0.0) THEN ! Stable stratification
psi_t = c_MO_t_stab*R_z*ZoL*(1.0-MIN(z0t_sf/height_tq, 1.0))
psi_q = c_MO_q_stab*R_z*ZoL*(1.0-MIN(z0q_sf/height_tq, 1.0))
!_dbg
! print*(*,*) 'STAB: psi_t = ', psi_t, ' psi_q = ', psi_q
!_dbg
ELSE ! Convection
psi_u = (1.0-c_MO_t_conv*R_z*ZoL)**c_MO_t_exp
psi_t = (1.0-c_MO_t_conv*R_z*ZoL*MIN(z0t_sf/height_tq, 1.0))**c_MO_t_exp
! psi_t = 2.0*LOG((1.0+psi_t)/(1.0+psi_u))
psi_t = abs(2.0*LOG((1.0+psi_t)/(1.0+psi_u)))
psi_u = (1.0-c_MO_q_conv*R_z*ZoL)**c_MO_q_exp
psi_q = (1.0-c_MO_q_conv*R_z*ZoL*MIN(z0q_sf/height_tq, 1.0))**c_MO_q_exp
! psi_q = 2.0*LOG((1.0+psi_q)/(1.0+psi_u))
psi_q = abs(2.0*LOG((1.0+psi_q)/(1.0+psi_u)))
! write(0,*) 'psi_q= ',psi_q
!_dbg
! print*(*,*) 'CONV: psi_t = ', psi_t, ' psi_q = ', psi_q
!_dbg
END IF
Q_sen_tur = -(T_a-T_s)*u_star_st*c_Karman/Pr_neutral &
/ MAX(c_small_sf, LOG(height_tq/z0t_sf)+psi_t)
if(MAX(c_small_sf, LOG(height_tq/z0t_sf)+psi_t) .lt. 10E-6) then
write(0,*)'inside flake'
write(0,*) Q_sen_tur, T_a, T_s, u_star_st, c_Karman, Pr_neutral
write(0,*) c_small_sf,height_tq,z0t_sf,psi_t
write(0,*) 'nominator= ', (T_a-T_s)*u_star_st*c_Karman/Pr_neutral
write(0,*) 'denominator= ',MAX(c_small_sf, LOG(height_tq/z0t_sf)+psi_t)
endif
Q_lat_tur = -(q_a-q_s)*u_star_st*c_Karman/Sc_neutral &
/ MAX(c_small_sf, LOG(height_tq/z0q_sf)+psi_q)
if(Q_lat_tur .gt. 6.0E-4) then
Q_lat_tur = -(q_a-q_s)*u_star_st*c_Karman/3.0 &
/ MAX(c_small_sf, LOG(height_tq/z0q_sf)+psi_q)
write(0,*) 'Q_lat_tur= ',Q_lat_tur
write(0,135) q_a,q_s,u_star_st,c_Karman
write(0,136) MAX(c_small_sf,LOG(height_tq/z0q_sf)+psi_q),c_small_sf, LOG(height_tq/z0q_sf),psi_q
endif
135 format(1x,4(f16.4))
136 format(1x,4(f16.4))
END IF Turb_Fluxes
!------------------------------------------------------------------------------
! Decide between turbulent, molecular, and convective fluxes
!------------------------------------------------------------------------------
Q_momentum = MIN(Q_mom_tur, Q_mom_mol, Q_mom_con) ! Momentum flux is negative
IF(l_conv_visc) THEN ! Convection, take fluxes that are maximal in magnitude
IF(ABS(Q_sen_tur).GE.ABS(Q_sen_con)) THEN
Q_sensible = Q_sen_tur
ELSE
Q_sensible = Q_sen_con
END IF
IF(ABS(Q_sensible).LT.ABS(Q_sen_mol)) THEN
Q_sensible = Q_sen_mol
END IF
IF(ABS(Q_lat_tur).GE.ABS(Q_lat_con)) THEN
Q_latent = Q_lat_tur
ELSE
Q_latent = Q_lat_con
END IF
IF(ABS(Q_latent).LT.ABS(Q_lat_mol)) THEN
Q_latent = Q_lat_mol
END IF
ELSE ! Stable or neutral stratification, chose fluxes that are maximal in magnitude
IF(ABS(Q_sen_tur).GE.ABS(Q_sen_mol)) THEN
Q_sensible = Q_sen_tur
ELSE
Q_sensible = Q_sen_mol
END IF
IF(ABS(Q_lat_tur).GE.ABS(Q_lat_mol)) THEN
Q_latent = Q_lat_tur
ELSE
Q_latent = Q_lat_mol
END IF
END IF
!------------------------------------------------------------------------------
! Set output (notice that fluxes are no longer in kinematic units)
!------------------------------------------------------------------------------
Q_momentum = Q_momentum*rho_a
!Q_sensible = Q_sensible*rho_a*tpsf_c_a_p
!write(0,*) 'Q_sensible= ',Q_sensible
Q_watvap = Q_latent*rho_a
!Q_latent = tpsf_L_evap
IF(h_ice.GE.h_Ice_min_flk) Q_latent = Q_latent + tpl_L_f ! Add latent heat of fusion over ice
!Q_latent = Q_watvap*Q_latent
Q_latent = Q_watvap*tpsf_L_evap
if(Q_latent .gt. 2000.00) then
write(0,145) 'final Q_watvap= ',Q_watvap, 'tpsf_L_evap= ',tpsf_L_evap, 'Q_latent= ', Q_latent
endif
!Q_latent = Q_watvap*Q_latent
145 format(A17,E12.5,1x,A13,1x,f10.2,1x,A10,1x,E12.4)
! Set "*_sf" variables to make fluxes accessible to driving routines that use "SfcFlx"
u_star_a_sf = u_star_st
Q_mom_a_sf = Q_momentum
Q_sens_a_sf = Q_sensible
Q_lat_a_sf = Q_latent
Q_watvap_a_sf = Q_watvap
!write(85,127) Q_sensible, Q_watvap, Q_latent
127 format(1x, 3(f16.5,1x))
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END SUBROUTINE SfcFlx_momsenlat
!==============================================================================
!==============================================================================
! include 'SfcFlx_rhoair.incf'
!------------------------------------------------------------------------------
REAL (KIND = kind_phys) FUNCTION SfcFlx_rhoair (T, q, P)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the air density as function
! of temperature, specific humidity and pressure.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T , & ! Temperature [K]
q , & ! Specific humidity
P ! Pressure [N m^{-2} = kg m^{-1} s^{-2}]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Air density [kg m^{-3}]
SfcFlx_rhoair = P/tpsf_R_dryair/T/(1.0+(1.0/tpsf_Rd_o_Rv-1.0)*q)
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_rhoair
!==============================================================================
!==============================================================================
! include 'SfcFlx_roughness.incf'
!------------------------------------------------------------------------------
SUBROUTINE SfcFlx_roughness (fetch, U_a, u_star, h_ice, &
c_z0u_fetch, u_star_thresh, z0u, z0t, z0q)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes the water-surface or the ice-surface roughness lengths
! with respect to wind velocity, potential temperature and specific humidity.
!
! The water-surface roughness lengths with respect to wind velocity is computed
! from the Charnock formula when the surface is aerodynamically rough.
! A simple empirical formulation is used to account for the dependence
! of the Charnock parameter on the wind fetch.
! When the flow is aerodynamically smooth, the roughness length with respect to
! wind velocity is proportional to the depth of the viscous sub-layer.
! The water-surface roughness lengths for scalars are computed using the power-law
! formulations in terms of the roughness Reynolds number (Zilitinkevich et al. 2001).
! The ice-surface aerodynamic roughness is taken to be constant.
! The ice-surface roughness lengths for scalars
! are computed through the power-law formulations
! in terms of the roughness Reynolds number (Andreas 2002).
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals , & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (procedure arguments)
REAL (KIND = kind_phys), INTENT(IN) :: &
fetch , & ! Typical wind fetch [m]
U_a , & ! Wind speed [m s^{-1}]
u_star , & ! Friction velocity in the surface air layer [m s^{-1}]
h_ice ! Ice thickness [m]
! Output (procedure arguments)
REAL (KIND = kind_phys), INTENT(OUT) :: &
c_z0u_fetch , & ! Fetch-dependent Charnock parameter
u_star_thresh , & ! Threshold value of friction velocity [m s^{-1}]
z0u , & ! Roughness length with respect to wind velocity [m]
z0t , & ! Roughness length with respect to potential temperature [m]
z0q ! Roughness length with respect to specific humidity [m]
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
Re_s , & ! Surface Reynolds number
Re_s_thresh ! Threshold value of Re_s
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
Water_or_Ice: IF(h_ice.LT.h_Ice_min_flk) THEN ! Water surface
! The Charnock parameter as dependent on dimensionless fetch
c_z0u_fetch = MAX(U_a, u_wind_min_sf)**2_iintegers/tpl_grav/fetch ! Inverse dimensionless fetch
c_z0u_fetch = c_z0u_rough + c_z0u_ftch_f*c_z0u_fetch**c_z0u_ftch_ex
c_z0u_fetch = MIN(c_z0u_fetch, c_z0u_rough_L) ! Limit Charnock parameter
! Threshold value of friction velocity
u_star_thresh = (c_z0u_smooth/c_z0u_fetch*tpl_grav*tpsf_nu_u_a)**num_1o3_sf
! Surface Reynolds number and its threshold value
Re_s = u_star**3_iintegers/tpsf_nu_u_a/tpl_grav
Re_s_thresh = c_z0u_smooth/c_z0u_fetch
! Aerodynamic roughness
IF(Re_s.LE.Re_s_thresh) THEN
z0u = c_z0u_smooth*tpsf_nu_u_a/u_star ! Smooth flow
ELSE
z0u = c_z0u_fetch*u_star*u_star/tpl_grav ! Rough flow
END IF
! Roughness for scalars
z0q = c_z0u_fetch*MAX(Re_s, Re_s_thresh)
z0t = c_z0t_rough_1*z0q**c_z0t_rough_3 - c_z0t_rough_2
z0q = c_z0q_rough_1*z0q**c_z0q_rough_3 - c_z0q_rough_2
z0t = z0u*EXP(-c_Karman/Pr_neutral*z0t)
z0q = z0u*EXP(-c_Karman/Sc_neutral*z0q)
ELSE Water_or_Ice ! Ice surface
! The Charnock parameter is not used over ice, formally set "c_z0u_fetch" to its minimum value
c_z0u_fetch = c_z0u_rough
! Threshold value of friction velocity
u_star_thresh = c_z0u_smooth*tpsf_nu_u_a/z0u_ice_rough
! Aerodynamic roughness
z0u = MAX(z0u_ice_rough, c_z0u_smooth*tpsf_nu_u_a/u_star)
! Roughness Reynolds number
Re_s = MAX(u_star*z0u/tpsf_nu_u_a, c_accur_sf)
! Roughness for scalars
IF(Re_s.LE.Re_z0s_ice_t) THEN
z0t = c_z0t_ice_b0t + c_z0t_ice_b1t*LOG(Re_s)
z0t = MIN(z0t, c_z0t_ice_b0s)
z0q = c_z0q_ice_b0t + c_z0q_ice_b1t*LOG(Re_s)
z0q = MIN(z0q, c_z0q_ice_b0s)
ELSE
z0t = c_z0t_ice_b0r + c_z0t_ice_b1r*LOG(Re_s) + c_z0t_ice_b2r*LOG(Re_s)**2_iintegers
z0q = c_z0q_ice_b0r + c_z0q_ice_b1r*LOG(Re_s) + c_z0q_ice_b2r*LOG(Re_s)**2_iintegers
END IF
z0t = z0u*EXP(z0t)
z0q = z0u*EXP(z0q)
END IF Water_or_Ice
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END SUBROUTINE SfcFlx_roughness
!==============================================================================
!==============================================================================
! include 'SfcFlx_satwvpres.incf'
!------------------------------------------------------------------------------
REAL (KIND = kind_phys) FUNCTION SfcFlx_satwvpres (T, h_ice)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes saturation water vapour pressure
! over the water surface or over the ice surface
! as function of temperature.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
!_dm The variable is USEd in module "SfcFlx".
!_nu USE flake_parameters , ONLY : &
!_nu h_Ice_min_flk ! Minimum ice thickness [m]
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T , & ! Temperature [K]
h_ice ! Ice thickness [m]
! Local parameters
REAL (KIND = kind_phys), PARAMETER :: &
b1_vap = 610.780 , & ! Coefficient [N m^{-2} = kg m^{-1} s^{-2}]
b3_vap = 273.160 , & ! Triple point [K]
b2w_vap = 17.26938820 , & ! Coefficient (water)
b2i_vap = 21.87455840 , & ! Coefficient (ice)
b4w_vap = 35.860 , & ! Coefficient (temperature) [K]
b4i_vap = 7.660 ! Coefficient (temperature) [K]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Saturation water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
IF(h_ice.LT.h_Ice_min_flk) THEN ! Water surface
SfcFlx_satwvpres = b1_vap*EXP(b2w_vap*(T-b3_vap)/(T-b4w_vap))
ELSE ! Ice surface
SfcFlx_satwvpres = b1_vap*EXP(b2i_vap*(T-b3_vap)/(T-b4i_vap))
END IF
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_satwvpres
!==============================================================================
!==============================================================================
! include 'SfcFlx_spechum.incf'
!------------------------------------------------------------------------------
REAL (KIND = kind_phys) FUNCTION SfcFlx_spechum (wvpres, P)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes specific humidity as function
! of water vapour pressure and air pressure.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
wvpres , & ! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
P ! Air pressure [N m^{-2} = kg m^{-1} s^{-2}]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Specific humidity
SfcFlx_spechum = tpsf_Rd_o_Rv*wvpres/(P-(1.0-tpsf_Rd_o_Rv)*wvpres)
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_spechum
!==============================================================================
!==============================================================================
! include 'SfcFlx_wvpreswetbulb.incf'
!------------------------------------------------------------------------------
REAL (KIND = ireals) FUNCTION SfcFlx_wvpreswetbulb (T_dry, T_wetbulb, satwvpres_bulb, P)
!------------------------------------------------------------------------------
!
! Description:
!
! Computes water vapour pressure as function of air temperature,
! wet bulb temperature, satururation vapour pressure at wet-bulb temperature,
! and air pressure.
!
! Declarations:
!
! Modules used:
!_dm Parameters are USEd in module "SfcFlx".
!_nu USE data_parameters , ONLY : &
!_nu ireals, & ! KIND-type parameter for real variables
!_nu iintegers ! KIND-type parameter for "normal" integer variables
use machine, only: kind_phys
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (function argument)
REAL (KIND = kind_phys), INTENT(IN) :: &
T_dry , & ! Dry air temperature [K]
T_wetbulb , & ! Wet bulb temperature [K]
satwvpres_bulb , & ! Satururation vapour pressure at wet-bulb temperature [N m^{-2}]
P ! Atmospheric pressure [N m^{-2}]
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! Water vapour pressure [N m^{-2} = kg m^{-1} s^{-2}]
SfcFlx_wvpreswetbulb = satwvpres_bulb &
- tpsf_c_a_p*P/tpsf_L_evap/tpsf_Rd_o_Rv*(T_dry-T_wetbulb)
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END FUNCTION SfcFlx_wvpreswetbulb
!==============================================================================
END MODULE SfcFlx
MODULE module_FLake
IMPLICIT NONE
CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE flake_interface ( dMsnowdt_in, I_atm_in, Q_atm_lw_in, height_u_in, height_tq_in, &
U_a_in, T_a_in, q_a_in, P_a_in, &
depth_w, fetch, depth_bs, T_bs, par_Coriolis, del_time, &
T_snow_in, T_ice_in, T_mnw_in, T_wML_in, T_bot_in, T_B1_in, &
C_T_in, h_snow_in, h_ice_in, h_ML_in, H_B1_in, T_sfc_p, &
ch, cm, albedo_water, water_extinc, &
T_snow_out, T_ice_out, T_mnw_out, T_wML_out, T_bot_out, &
T_B1_out, C_T_out, h_snow_out, h_ice_out, h_ML_out, &
H_B1_out, T_sfc_n, hflx_out, evap_out, gflx_out, lflx_out, &
T_bot_2_in, T_bot_2_out,ustar, q_sfc, chh, cmm )
!------------------------------------------------------------------------------
!
! Description:
!
! The FLake interface is
! a communication routine between "flake_main"
! and a prediction system that uses FLake.
! It assigns the FLake variables at the previous time step
! to their input values given by the driving model,
! calls a number of routines to compute the heat and radiation fluxes,
! calls "flake_main",
! and returns the updated FLake variables to the driving model.
! The "flake_interface" does not contain any Flake physics.
! It only serves as a convenient means to organize calls of "flake_main"
! and of external routines that compute heat and radiation fluxes.
! The interface may (should) be changed so that to provide
! the most convenient use of FLake.
! Within a 3D atmospheric prediction system,
! "flake_main" may be called in a DO loop within "flake_interface"
! for each grid-point where a lake is present.
! In this way, the driving atmospheric model should call "flake_interface"
! only once, passing the FLake variables to "flake_interface" as 2D fields.
!
! Lines embraced with "!_tmp" contain temporary parts of the code.
! These should be removed prior to using FLake in applications.
! Lines embraced/marked with "!_dev" may be replaced
! as improved parameterizations are developed and tested.
! Lines embraced/marked with "!_dm" are DM's comments
! that may be helpful to a user.
!
use machine, only: kind_phys
USE data_parameters , ONLY : &
ireals, & ! KIND-type parameter for real variables
iintegers ! KIND-type parameter for "normal" integer variables
USE flake_derivedtypes ! Definitions of several derived TYPEs
!USE flake_parameters , ONLY : &
! tpl_T_f , & ! Fresh water freezing point [K]
! tpl_rho_w_r , & ! Maximum density of fresh water [kg m^{-3}]
! h_Snow_min_flk , & ! Minimum snow thickness [m]
! h_Ice_min_flk ! Minimum ice thickness [m]
USE flake_paramoptic_ref ! Reference values of the optical characteristics
! of the lake water, lake ice and snow
USE flake_albedo_ref ! Reference values the albedo for the lake water, lake ice and snow
USE flake , ONLY : &
flake_main , & ! Subroutine, FLake driver
flake_radflux , & ! Subroutine, computes radiation fluxes at various depths
!
T_snow_p_flk, T_snow_n_flk , & ! Temperature at the air-snow interface [K]
T_ice_p_flk, T_ice_n_flk , & ! Temperature at the snow-ice or air-ice interface [K]
T_mnw_p_flk, T_mnw_n_flk , & ! Mean temperature of the water column [K]
T_wML_p_flk, T_wML_n_flk , & ! Mixed-layer temperature [K]
T_bot_p_flk, T_bot_n_flk , & ! Temperature at the water-bottom sediment interface [K]
T_B1_p_flk, T_B1_n_flk , & ! Temperature at the bottom of the upper layer of the sediments [K]
C_T_p_flk, C_T_n_flk , & ! Shape factor (thermocline)
h_snow_p_flk, h_snow_n_flk , & ! Snow thickness [m]
h_ice_p_flk, h_ice_n_flk , & ! Ice thickness [m]
h_ML_p_flk, h_ML_n_flk , & ! Thickness of the mixed-layer [m]
H_B1_p_flk, H_B1_n_flk , & ! Thickness of the upper layer of bottom sediments [m]
!
Q_snow_flk , & ! Heat flux through the air-snow interface [W m^{-2}]
Q_ice_flk , & ! Heat flux through the snow-ice or air-ice interface [W m^{-2}]
Q_w_flk , & ! Heat flux through the ice-water or air-water interface [W m^{-2}]
Q_bot_flk , & ! Heat flux through the water-bottom sediment interface [W m^{-2}]
I_atm_flk , & ! Radiation flux at the lower boundary of the atmosphere [W m^{-2}],
! i.e. the incident radiation flux with no regard for the surface albedo
I_snow_flk , & ! Radiation flux through the air-snow interface [W m^{-2}]
I_ice_flk , & ! Radiation flux through the snow-ice or air-ice interface [W m^{-2}]
I_w_flk , & ! Radiation flux through the ice-water or air-water interface [W m^{-2}]
I_h_flk , & ! Radiation flux through the mixed-layer-thermocline interface [W m^{-2}]
I_bot_flk , & ! Radiation flux through the water-bottom sediment interface [W m^{-2}]
I_intm_0_h_flk , & ! Mean radiation flux over the mixed layer [W m^{-1}]
I_intm_h_D_flk , & ! Mean radiation flux over the thermocline [W m^{-1}]
Q_star_flk , & ! A generalized heat flux scale [W m^{-2}]
u_star_w_flk , & ! Friction velocity in the surface layer of lake water [m s^{-1}]
w_star_sfc_flk , & ! Convective velocity scale, using a generalized heat flux scale [m s^{-1}]
dMsnowdt_flk , & ! The rate of snow accumulation [kg m^{-2} s^{-1}]
T_bot_2_in_flk
USE SfcFlx , ONLY : &
SfcFlx_lwradwsfc , & ! Function, returns the surface long-wave radiation flux
SfcFlx_momsenlat ! Subroutine, computes fluxes of momentum and of sensible and latent heat
!==============================================================================
IMPLICIT NONE
!==============================================================================
!
! Declarations
! Input (procedure arguments)
REAL (KIND = kind_phys), INTENT(IN) :: &
dMsnowdt_in , & ! The rate of snow accumulation [kg m^{-2} s^{-1}]
I_atm_in , & ! Solar radiation flux at the surface [W m^{-2}]
Q_atm_lw_in , & ! Long-wave radiation flux from the atmosphere [W m^{-2}]
height_u_in , & ! Height above the lake surface where the wind speed is measured [m]
height_tq_in , & ! Height where temperature and humidity are measured [m]
U_a_in , & ! Wind speed at z=height_u_in [m s^{-1}]
T_a_in , & ! Air temperature at z=height_tq_in [K]
q_a_in , & ! Air specific humidity at z=height_tq_in
P_a_in , & ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}]
ch , &
cm , &
albedo_water, & ! Water surface albedo with respect to the solar radiation
water_extinc
REAL (KIND = kind_phys), INTENT(IN) :: &
depth_w , & ! The lake depth [m]
fetch , & ! Typical wind fetch [m]
depth_bs , & ! Depth of the thermally active layer of the bottom sediments [m]
T_bs , & ! Temperature at the outer edge of
! the thermally active layer of the bottom sediments [K]
par_Coriolis , & ! The Coriolis parameter [s^{-1}]
del_time ! The model time step [s]
REAL (KIND = kind_phys), INTENT(IN) :: &
T_snow_in , & ! Temperature at the air-snow interface [K]
T_ice_in , & ! Temperature at the snow-ice or air-ice interface [K]
T_mnw_in , & ! Mean temperature of the water column [K]
T_wML_in , & ! Mixed-layer temperature [K]
T_bot_in , & ! Temperature at the water-bottom sediment interface [K]
T_B1_in , & ! Temperature at the bottom of the upper layer of the sediments [K]
C_T_in , & ! Shape factor (thermocline)
h_snow_in , & ! Snow thickness [m]
h_ice_in , & ! Ice thickness [m]
h_ML_in , & ! Thickness of the mixed-layer [m]
H_B1_in , & ! Thickness of the upper layer of bottom sediments [m]
T_sfc_p , & ! Surface temperature at the previous time step [K]
T_bot_2_in
! Input/Output (procedure arguments)
!REAL (KIND = ireals), INTENT(INOUT) :: &
REAL (KIND = kind_phys) :: &
albedo_ice , & ! Ice surface albedo with respect to the solar radiation
albedo_snow ! Snow surface albedo with respect to the solar radiation
!TYPE (opticpar_medium), INTENT(INOUT) :: &
TYPE (opticpar_medium) :: &
opticpar_water , & ! Optical characteristics of water
opticpar_ice , & ! Optical characteristics of ice
opticpar_snow ! Optical characteristics of snow
! Output (procedure arguments)
REAL (KIND = kind_phys), INTENT(OUT) :: &
T_snow_out , & ! Temperature at the air-snow interface [K]
T_ice_out , & ! Temperature at the snow-ice or air-ice interface [K]
T_mnw_out , & ! Mean temperature of the water column [K]
T_wML_out , & ! Mixed-layer temperature [K]
T_bot_out , & ! Temperature at the water-bottom sediment interface [K]
T_B1_out , & ! Temperature at the bottom of the upper layer of the sediments [K]
C_T_out , & ! Shape factor (thermocline)
h_snow_out , & ! Snow thickness [m]
h_ice_out , & ! Ice thickness [m]
h_ML_out , & ! Thickness of the mixed-layer [m]
H_B1_out , & ! Thickness of the upper layer of bottom sediments [m]
T_sfc_n , & ! Updated surface temperature [K]
hflx_out , & ! sensibl heat flux
evap_out , & ! Latent heat flux
gflx_out , & ! flux from to water
lflx_out , & ! latent heat flux
T_bot_2_out , & ! Bottom temperature
ustar , &
q_sfc , &
chh , &
cmm
! Local variables of type REAL
REAL (KIND = kind_phys) :: &
Q_momentum , & ! Momentum flux [N m^{-2}]
Q_sensible , & ! Sensible heat flux [W m^{-2}]
Q_latent , & ! Latent heat flux [W m^{-2}]
Q_watvap , & ! Flux of water vapour [kg m^{-2} s^{-1}]
Q_w_flux , & ! flux from ice to water
rho_a
! ADDED by Shaobo Zhang
LOGICAL lflk_botsed_use
!REAL (KIND = kind_phys) :: T_bot_2_in, T_bot_2_out
REAL (KIND = kind_phys), parameter :: tpl_rho_w_r = 1.0E+03
REAL (KIND = kind_phys), parameter :: tpl_T_f = 273.15
REAL (KIND = kind_phys), parameter :: h_Snow_min_flk = 1.0E-5
REAL (KIND = kind_phys), parameter :: h_Ice_min_flk = 1.0E-9
!==============================================================================
! Start calculations
!------------------------------------------------------------------------------
! lflk_botsed_use = .TRUE.
lflk_botsed_use = .FALSE.
!------------------------------------------------------------------------------
! Set albedos of the lake water, lake ice and snow
!------------------------------------------------------------------------------
! Use default value
! albedo_water = albedo_water_ref
! Use empirical formulation proposed by Mironov and Ritter (2004) for GME
!_nu albedo_ice = albedo_whiteice_ref
!albedo_ice = EXP(-c_albice_MR*(tpl_T_f-T_sfc_p)/tpl_T_f)
!albedo_ice = albedo_whiteice_ref*(1.0-albedo_ice) + albedo_blueice_ref*albedo_ice
! Snow is not considered
!albedo_snow = albedo_ice
albedo_ice = albedo_whiteice_ref
!albedo_snow = albedo_ice
albedo_snow = albedo_drysnow_ref
opticpar_water%extincoef_optic(1) = water_extinc
!write(0,*)'albedo= ',albedo_water,albedo_ice,albedo_snow
!------------------------------------------------------------------------------
! Set optical characteristics of the lake water, lake ice and snow
!------------------------------------------------------------------------------
! Use default values
opticpar_water = opticpar_water_ref
opticpar_ice = opticpar_ice_opaque ! Opaque ice
opticpar_snow = opticpar_snow_opaque ! Opaque snow
!print*,'opticpar = ',opticpar_water, opticpar_ice,opticpar_snow
!------------------------------------------------------------------------------
! Set initial values
!------------------------------------------------------------------------------
!print*,'Inter depth_w=',depth_w
!print*,'Inter depth_bs=',depth_bs
T_snow_p_flk = T_snow_in
T_ice_p_flk = T_ice_in
T_mnw_p_flk = T_mnw_in
T_wML_p_flk = T_wML_in
T_bot_p_flk = T_bot_in
T_B1_p_flk = T_B1_in
C_T_p_flk = C_T_in
h_snow_p_flk = h_snow_in
h_ice_p_flk = h_ice_in
h_ML_p_flk = h_ML_in
H_B1_p_flk = H_B1_in
T_bot_2_in_flk = T_bot_2_in
!write(71,120) T_sfc_p,T_mnw_in,T_wML_in,T_bot_in,T_B1_in,T_bot_2_in
120 format(1x,6(f12.5,1x))
!------------------------------------------------------------------------------
! Set the rate of snow accumulation
!------------------------------------------------------------------------------
dMsnowdt_flk = dMsnowdt_in
!------------------------------------------------------------------------------
! Compute solar radiation fluxes (positive downward)
!------------------------------------------------------------------------------
I_atm_flk = I_atm_in
CALL flake_radflux ( depth_w, albedo_water, albedo_ice, albedo_snow, &
opticpar_water, opticpar_ice, opticpar_snow, &
depth_bs )
!------------------------------------------------------------------------------
! Compute long-wave radiation fluxes (positive downward)
!------------------------------------------------------------------------------
Q_w_flk = Q_atm_lw_in ! Radiation of the atmosphere
Q_w_flk = Q_w_flk - SfcFlx_lwradwsfc(T_sfc_p) ! Radiation of the surface (notice the sign)
!------------------------------------------------------------------------------
! Compute the surface friction velocity and fluxes of sensible and latent heat
!------------------------------------------------------------------------------
CALL SfcFlx_momsenlat ( height_u_in, height_tq_in, fetch, &
U_a_in, T_a_in, q_a_in, T_sfc_p, P_a_in, h_ice_p_flk, &
Q_momentum, Q_sensible, Q_latent, Q_watvap, q_sfc, rho_a )
!write(0,*)'tpl_rho_w_r= ', tpl_rho_w_r
!write(0,*) 'Q_momentum= ',Q_momentum
u_star_w_flk = SQRT(-Q_momentum/tpl_rho_w_r)
ustar = u_star_w_flk
!------------------------------------------------------------------------------
! Compute heat fluxes Q_snow_flk, Q_ice_flk, Q_w_flk
!------------------------------------------------------------------------------
Q_w_flk = Q_w_flk - Q_sensible - Q_latent ! Add sensible and latent heat fluxes (notice the signs)
IF(h_ice_p_flk.GE.h_Ice_min_flk) THEN ! Ice exists
IF(h_snow_p_flk.GE.h_Snow_min_flk) THEN ! There is snow above the ice
Q_snow_flk = Q_w_flk
Q_ice_flk = 0.0
Q_w_flk = 0.0
ELSE ! No snow above the ice
Q_snow_flk = 0.0
Q_ice_flk = Q_w_flk
Q_w_flk = 0.0
END IF
ELSE ! No ice-snow cover
Q_snow_flk = 0.0
Q_ice_flk = 0.0
END IF
!------------------------------------------------------------------------------
! Advance FLake variables
!------------------------------------------------------------------------------
CALL flake_main ( depth_w, depth_bs, T_bs, par_Coriolis, &
opticpar_water%extincoef_optic(1), &
del_time, T_sfc_p, T_sfc_n, T_bot_2_in_flk, &
T_bot_2_out )
!------------------------------------------------------------------------------
! Set output values
!------------------------------------------------------------------------------
T_snow_out = T_snow_n_flk
T_ice_out = T_ice_n_flk
T_mnw_out = T_mnw_n_flk
T_wML_out = T_wML_n_flk
T_bot_out = T_bot_n_flk
T_B1_out = T_B1_n_flk
C_T_out = C_T_n_flk
h_snow_out = h_snow_n_flk
h_ice_out = h_ice_n_flk
h_ML_out = h_ML_n_flk
H_B1_out = H_B1_n_flk
hflx_out = Q_sensible
evap_out = Q_watvap
!evap_out = Q_latent
gflx_out = Q_w_flk
lflx_out = Q_latent
chh = ch * U_a_in * rho_a
cmm = cm * U_a_in
!write(72,120) T_sfc_n,T_mnw_out,T_wML_out,T_bot_out,T_B1_out,T_bot_2_out
!------------------------------------------------------------------------------
! End calculations
!==============================================================================
END SUBROUTINE flake_interface
END MODULE module_FLake