diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index d5bc98322..1e04a9d44 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -24,8 +24,8 @@ end subroutine GFS_surface_composites_pre_finalize !> \section arg_table_GFS_surface_composites_pre_run Argument Table !! \htmlinclude GFS_surface_composites_pre_run.html !! - subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cplwav2atm, & - landfrac, lakefrac, oceanfrac, & + subroutine GFS_surface_composites_pre_run (im, lkm, frac_grid, flag_cice, cplflx, cplwav2atm, & + landfrac, lakefrac, lakedepth, oceanfrac, & frland, dry, icy, lake, ocean, wet, cice, cimin, zorl, zorlo, zorll, zorl_wat, & zorl_lnd, zorl_ice, snowd, snowd_wat, snowd_lnd, snowd_ice, tprcp, tprcp_wat, & tprcp_lnd, tprcp_ice, uustar, uustar_wat, uustar_lnd, uustar_ice, & @@ -38,12 +38,12 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cpl implicit none ! Interface variables - integer, intent(in ) :: im + integer, intent(in ) :: im, lkm logical, intent(in ) :: frac_grid, cplflx, cplwav2atm logical, dimension(im), intent(in ) :: flag_cice logical, dimension(im), intent(inout) :: dry, icy, lake, ocean, wet real(kind=kind_phys), intent(in ) :: cimin - real(kind=kind_phys), dimension(im), intent(in ) :: landfrac, lakefrac, oceanfrac + real(kind=kind_phys), dimension(im), intent(in ) :: landfrac, lakefrac, lakedepth, oceanfrac real(kind=kind_phys), dimension(im), intent(inout) :: cice real(kind=kind_phys), dimension(im), intent( out) :: frland real(kind=kind_phys), dimension(im), intent(in ) :: zorl, snowd, tprcp, uustar, weasd, qss, hflx @@ -182,6 +182,19 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cpl endif enddo +! to prepare to separate lake from ocean under water category + do i = 1, im + if(lkm == 1) then + if(lakefrac(i) .ge. 0.15 .and. lakedepth(i) .gt. 1.0) then + lake(i) = .true. + else + lake(i) = .false. + endif + else + lake(i) = .false. + endif + enddo + ! Assign sea ice temperature to interstitial variable do i = 1, im tice(i) = tisfc(i) diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 6c9fb5ba0..9b297ca38 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -9,6 +9,14 @@ type = integer intent = in optional = F +[lkm] + standard_name = flag_for_lake_surface_scheme + long_name = flag for lake surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F [frac_grid] standard_name = flag_for_fractional_grid long_name = flag for fractional grid @@ -59,6 +67,15 @@ kind = kind_phys intent = in optional = F +[lakedepth] + standard_name = lake_depth + long_name = lake depth + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F [oceanfrac] standard_name = sea_area_fraction long_name = fraction of horizontal grid area occupied by ocean diff --git a/physics/GFS_time_vary_pre.fv3.F90 b/physics/GFS_time_vary_pre.fv3.F90 index 98a0f6697..5f72a6b27 100644 --- a/physics/GFS_time_vary_pre.fv3.F90 +++ b/physics/GFS_time_vary_pre.fv3.F90 @@ -65,8 +65,8 @@ end subroutine GFS_time_vary_pre_finalize !> \section arg_table_GFS_time_vary_pre_run Argument Table !! \htmlinclude GFS_time_vary_pre_run.html !! - subroutine GFS_time_vary_pre_run (jdat, idat, dtp, lsm, lsm_noahmp, nsswr, & - nslwr, nhfrad, idate, debug, me, master, nscyc, sec, phour, zhour, fhour, & + subroutine GFS_time_vary_pre_run (jdat, idat, dtp, lkm, lsm, lsm_noahmp, nsswr, & + nslwr, nhfrad, idate, debug, me, master, nscyc, sec, phour, zhour, fhour, & kdt, julian, yearlen, ipt, lprnt, lssav, lsswr, lslwr, solhr, errmsg, errflg) use machine, only: kind_phys @@ -75,7 +75,7 @@ subroutine GFS_time_vary_pre_run (jdat, idat, dtp, lsm, lsm_noahmp, nsswr, & integer, intent(in) :: idate(4) integer, intent(in) :: jdat(1:8), idat(1:8) - integer, intent(in) :: lsm, lsm_noahmp, & + integer, intent(in) :: lkm, lsm, lsm_noahmp, & nsswr, nslwr, me, & master, nscyc, nhfrad logical, intent(in) :: debug @@ -121,7 +121,8 @@ subroutine GFS_time_vary_pre_run (jdat, idat, dtp, lsm, lsm_noahmp, nsswr, & fhour = (sec + dtp)/con_hr kdt = nint((sec + dtp)/dtp) - if(lsm == lsm_noahmp) then + if(lsm == lsm_noahmp .or. lkm == 1) then +! flake need this too !GJF* These calculations were originally in GFS_physics_driver.F90 for ! NoahMP. They were moved to this routine since they only depend ! on time (not space). Note that this code is included as-is from diff --git a/physics/GFS_time_vary_pre.fv3.meta b/physics/GFS_time_vary_pre.fv3.meta index 14081f8e4..04f7f1529 100644 --- a/physics/GFS_time_vary_pre.fv3.meta +++ b/physics/GFS_time_vary_pre.fv3.meta @@ -70,6 +70,14 @@ kind = kind_phys intent = in optional = F +[lkm] + standard_name = flag_for_lake_surface_scheme + long_name = flag for lake surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F [lsm] standard_name = flag_for_land_surface_scheme long_name = flag for land surface model diff --git a/physics/flake.F90 b/physics/flake.F90 new file mode 100644 index 000000000..2c2e7218c --- /dev/null +++ b/physics/flake.F90 @@ -0,0 +1,3281 @@ + +!======================================================================= +! 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! +! +! +! 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. +! + +IMPLICIT NONE + +!======================================================================= +! Global (i.e. public) Declarations: +! Parameters for the Program: + + INTEGER, PARAMETER :: & + ireals = SELECTED_REAL_KIND (12,200), & + ! number of desired significant digits for + ! real variables + ! corresponds to 8 byte real variables + + iintegers = KIND (1) + ! kind-type parameter of the integer values + ! corresponds to the default integers + +!======================================================================= + +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_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_iceRi_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_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)) +!_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) +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) + +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 + +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 + +! 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.9,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, & + + 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 + 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}] + rho_a + +! ADDED by Shaobo Zhang +LOGICAL lflk_botsed_use +!REAL (KIND = kind_phys) :: T_bot_2_in, T_bot_2_out + +!============================================================================== +! Start calculations +!------------------------------------------------------------------------------ + lflk_botsed_use = .TRUE. +!------------------------------------------------------------------------------ +! 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 +opticpar_water%extincoef_optic(1) = water_extinc +!print*,'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 ) + +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 +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 diff --git a/physics/flake_driver.F90 b/physics/flake_driver.F90 new file mode 100644 index 000000000..b882c7404 --- /dev/null +++ b/physics/flake_driver.F90 @@ -0,0 +1,393 @@ +!> \file flake_driver.F90 +!! This file contains the flake scheme driver. + +!> This module contains the CCPP-compliant flake scheme driver. + module flake_driver + + implicit none + + private + + public :: flake_driver_init, flake_driver_run, flake_driver_finalize + + contains + +!> \section arg_table_flake_driver_init Argument Table +!! \htmlinclude flake_driver_init.html +!! + subroutine flake_driver_init (errmsg, errflg) + + implicit none + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + end subroutine flake_driver_init + +!> \section arg_table_flake_driver_finalize Argument Table +!! \htmlinclude flake_driver_finalize.html +!! + subroutine flake_driver_finalize (errmsg, errflg) + + implicit none + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + end subroutine flake_driver_finalize + +!> \section arg_table_flake_driver_run Argument Table +!! \htmlinclude flake_driver_run.html +!! + SUBROUTINE flake_driver_run ( & +! ---- Inputs + im, ps, t1, q1, wind, & + dlwflx, dswsfc, weasd, lakedepth, & + lake, xlat, delt, zlvl, elev, & + wet, flag_iter, yearlen, julian, imon, & +! ---- in/outs + snwdph, hice, tsurf, fice, T_sfc, hflx, evap, & + ustar, qsfc, ch, cm, chh, cmm, & + errmsg, errflg ) + +!============================================================================== +! +! Declarations +! use module_flake_ini, only:flake_init + use module_FLake +! use flake_albedo_ref +! use data_parameters +! use flake_derivedtypes +! use flake_paramoptic_ref +! use flake_parameters + use machine , only : kind_phys +! use funcphys, only : fpvs +! use physcons, only : grav => con_g, cp => con_cp, & +! & hvap => con_hvap, rd => con_rd, & +! & eps => con_eps, epsm1 => con_epsm1, & +! & rvrdm1 => con_fvirt + +!============================================================================== +IMPLICIT NONE + + integer, intent(in) :: im, imon,yearlen +! integer, dimension(im), intent(in) :: islmsk + + real (kind=kind_phys), dimension(im), intent(in) :: ps, wind, & + & t1, q1, dlwflx, dswsfc, zlvl, elev + + real (kind=kind_phys), intent(in) :: delt + + real (kind=kind_phys), dimension(im), intent(in) :: & + & xlat, weasd, lakedepth + + real (kind=kind_phys),dimension(im),intent(inout) :: & + & snwdph, hice, tsurf, t_sfc, hflx, evap, fice, ustar, qsfc, & + & ch, cm, chh, cmm + + real (kind=kind_phys), intent(in) :: julian + + logical, dimension(im), intent(in) :: flag_iter, wet, lake + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + +! --- locals + + real (kind=kind_phys) , parameter :: lake_pct_min = 0.1 + + real (kind=kind_phys), dimension(im) :: & + T_snow , & ! Temperature at the air-snow interface [K] + T_ice , & ! Temperature at the snow-ice or air-ice interface [K] + T_mnw , & ! Mean temperature of the water column [K] + T_wML , & ! Mixed-layer temperature [K] + T_bot , & ! Temperature at the water-bottom sediment interface [K] + T_B1 , & ! Temperature at the upper layer of the sediments [K] + C_T , & ! Shape factor (thermocline) + fetch , & ! Typical wind fetch [m] + h_ML , & ! Thickness of the mixed-layer [m] + H_B1 , & ! Thickness of the upper layer of bottom sediments [m] + w_albedo , & ! + w_extinc + +! Input (procedure arguments) + +REAL (KIND = kind_phys) :: & + + 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}] + +REAL (KIND = kind_phys) :: & + depth_w , & ! The lake depth [m] + fetch_in , & ! Typical wind fetch [m] + depth_bs_in , & ! Depth of the thermally active layer of the bottom sediments [m] + T_bs_in , & ! 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) :: & + 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_in , & ! Surface temperature at the previous time step [K] + ch_in , & + cm_in , & + albedo_water , & + water_extinc + +REAL (KIND = kind_phys) :: & + 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_out , & ! surface temperature [K] + T_sfc_n , & ! Updated surface temperature [K] + u_star , & + q_sfc , & + chh_out , & + cmm_out + +REAL (KIND = kind_phys) :: & + Q_momentum , & ! Momentum flux [N m^{-2}] + Q_SHT_flx , & ! Sensible heat flux [W m^{-2}] + Q_LHT_flx , & ! Latent heat flux [W m^{-2}] + Q_watvap ! Flux of water vapour [kg m^{-2} s^{-1}] + +REAL (KIND = kind_phys) :: & + lake_depth_max, T_bot_2_in, T_bot_2_out, dxlat,tb,tr,tt,temp,Kbar, DelK + +INTEGER :: i,ipr,iter + +LOGICAL :: lflk_botsed_use +logical :: flag(im) +CHARACTER(LEN=*), PARAMETER :: FMT2 = "(1x,8(F12.4,1x))" + +!============================================================================== +! Start calculations +!------------------------------------------------------------------------------ +! FLake_write need to assign original value to make the model somooth + + lake_depth_max = 60.0 + ipr = min(im,10) + +! --- ... set flag for lake points + + do i = 1, im + flag(i) = (wet(i) .and. flag_iter(i)) + enddo + + Kbar=3.5 + DelK=3.0 + + do i = 1, im + if (flag(i)) then + if( lake(i) ) then + T_ice(i) = 273.15 + T_snow(i) = 273.15 + fetch(i) = 2.0E+03 + C_T(i) = 0.50 + + dxlat = 57.29578*abs(xlat(i)) + tt = 29.275+0.0813*dxlat-0.0052*dxlat*dxlat-0.0038*elev(i)+273.15 + tb = 29.075-0.7566*dxlat+0.0051*dxlat*dxlat-0.0038*elev(i)+273.15 +! if(fice(i).le.0.0) then +! h_ice(i) = 0.0 +! h_snow(i)= 0.0 +! endif + if(snwdph(i).gt.0.0 .or. hice(i).gt.0.0) then + if(tsurf(i).lt.T_ice(i)) then + T_sfc(i) = T_ice(i) + else + T_sfc(i) = tsurf(i) + endif + else +! if(tsurf(i).lt.tt) then +! T_sfc(i) = tt +! else +! T_sfc(i) = tsurf(i) +! endif + T_sfc(i) = 0.2*tt + 0.8* tsurf(i) + endif + + T_bot(i) = tb + T_B1(i) = tb + +! if(lakedepth(i).lt.10.0) then +! T_bot(i) = T_sfc(i) +! T_B1(i) = T_bot(i) +! endif + + T_mnw(i) = C_T(i)*T_sfc(i)+(1-C_T(i))*T_bot(i) + T_wML(i) = C_T(i)*T_sfc(i)+(1-C_T(i))*T_bot(i) + h_ML(i) = C_T(i)* min ( lakedepth(i), lake_depth_max ) + H_B1(i) = min ( lakedepth(i),4.0) + hflx(i) = 0.0 + evap(i) = 0.0 + +! compute albedo as a function of julian day and latitute + temp = 2*3.14159265*(julian-1)/float(yearlen) + temp = 0.006918-0.399912*cos(temp)+0.070257*sin(temp)- & + 0.006758*cos(2.0*temp)+0.000907*sin(2.0*temp) - & + 0.002697*cos(3.0*temp)+0.00148*sin(3.0*temp) + w_albedo(I) = 0.06/cos((xlat(i)-temp)/1.2) +! w_albedo(I) = 0.06 +! compute water extinction coefficient as a function of julian day + if(julian.lt.90 .or. julian .gt. 333) then + w_extinc(i) = Kbar-Kbar/DelK + else + w_extinc(i) = Kbar+Kbar/DelK*sin(2*3.14159265*(julian-151)/244) + endif +! w_extinc(i) = 3.0 + +! write(65,1002) julian,xlat(i),w_albedo(I),w_extinc(i),lakedepth(i),elev(i),tb,tt,tsurf(i),T_sfc(i) +! print 1002 julian,xlat(i),w_albedo(I),w_extinc(i),lakedepth(i),elev(i),tb,tt,tsurf(i),T_sfc(i) +! print*,'inside flake driver' +! print*, julian,xlat(i),w_albedo(I),w_extinc(i),lakedepth(i),elev(i),tb,tt,tsurf(i),T_sfc(i) + + endif !lake fraction and depth + endif !flag + enddo + 1001 format ( 'At icount=', i5, ' x = ', f5.2,5x, 'y = ', & + 1p, e12.3) +! 1002 format ( ' julian= ',F6.2,1x,5(F8.4,1x),3(f11.4,1x)) + 1002 format (I4,1x,3(f8.4,1x),6(f11.4,1x)) + + +! +! call lake interface + do i=1,im + if (flag(i)) then + if( lake(i) ) then + dMsnowdt_in = weasd(i)/delt + I_atm_in = dswsfc(i) + Q_atm_lw_in = dlwflx(i) + height_u_in = zlvl(i) + height_tq_in = zlvl(i) + U_a_in = wind(i) + T_a_in = t1(i) + q_a_in = q1(i) + P_a_in = ps(i) + ch_in = ch(i) + cm_in = cm(i) + albedo_water= w_albedo(i) + water_extinc= w_extinc(i) + + depth_w = min ( lakedepth(i), lake_depth_max ) + depth_bs_in = max ( 4.0, min ( depth_w * 0.2, 10.0 ) ) + fetch_in = fetch(i) + T_bs_in = T_bot(i) + par_Coriolis = 2 * 7.2921 / 100000. * sin ( xlat(i) ) + del_time = delt + + do iter=1,10 !interation loop + T_snow_in = T_snow(i) + T_ice_in = T_ice(i) + T_mnw_in = T_mnw(i) + T_wML_in = T_wML(i) + T_bot_in = T_bot(i) + T_B1_in = T_B1(i) + C_T_in = C_T(i) + h_snow_in = snwdph(i) + h_ice_in = hice(i) + h_ML_in = h_ML(i) + H_B1_in = H_B1(i) + T_sfc_in = T_sfc(i) + + T_bot_2_in = T_bot(i) + Q_SHT_flx = hflx(i) + Q_watvap = evap(i) + +!------------------------------------------------------------------------------ +! Set the rate of snow accumulation +!------------------------------------------------------------------------------ + + CALL 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_in, depth_bs_in, T_bs_in, 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_in, & + ch_in, cm_in, 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_out, Q_SHT_flx, Q_watvap, & +! + T_bot_2_in, T_bot_2_out,u_star, q_sfc,chh_out,cmm_out ) + +!------------------------------------------------------------------------------ +! Update output and values for previous time step +! + T_snow(i) = T_snow_out + T_ice(i) = T_ice_out + T_mnw(i) = T_mnw_out + T_wML(i) = T_wML_out + T_sfc(i) = T_sfc_out + Tsurf(i) = T_sfc_out + T_bot(i) = T_bot_out + T_B1(i) = T_B1_out + C_T(i) = C_T_out + h_ML(i) = h_ML_out + H_B1(i) = H_B1_out + ustar(i) = u_star + qsfc(i) = q_sfc + chh(i) = chh_out + cmm(i) = cmm_out + snwdph(i) = h_snow_out + hice(i) = h_ice_out + evap(i) = Q_watvap + hflx(i) = Q_SHT_flx + + if(hice(i) .gt. 0.0 .or. snwdph(i) .gt. 0.0) then + fice(i) = 1.0 + else + fice(i) = 0.0 + endif + enddo !iter loop + endif !endif of lake + endif !endif of flag + + ENDDO + + 125 format(1x,i2,1x,i2,1x,i2,1x,6(1x,f14.8)) + 126 format(1x,i2,1x,i2,1x,6(1x,f14.8)) + 127 format(1x,i2,2(1x,f16.9)) +!------------------------------------------------------------------------------ +! End calculations +!============================================================================== + +END SUBROUTINE flake_driver_run + +!--------------------------------- + end module flake_driver diff --git a/physics/flake_driver.meta b/physics/flake_driver.meta new file mode 100644 index 000000000..a40016010 --- /dev/null +++ b/physics/flake_driver.meta @@ -0,0 +1,346 @@ +[ccpp-arg-table] + name = flake_driver_init + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F + +######################################################################## +[ccpp-arg-table] + name = flake_driver_finalize + type = scheme +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F + +######################################################################## +[ccpp-arg-table] + name = flake_driver_run + type = scheme +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[ps] + standard_name = surface_air_pressure + long_name = surface pressure + units = Pa + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[t1] + standard_name = air_temperature_at_lowest_model_layer + long_name = mean temperature at lowest model layer + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[q1] + standard_name = water_vapor_specific_humidity_at_lowest_model_layer + long_name = water vapor specific humidity at lowest model layer + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[wind] + standard_name = wind_speed_at_lowest_model_layer + long_name = wind speed at lowest model level + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[dlwflx] + standard_name = surface_downwelling_longwave_flux_absorbed_by_ground_over_ocean + long_name = total sky surface downward longwave flux absorbed by the ground over ocean + units = W m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[dswsfc] + standard_name = surface_downwelling_shortwave_flux + long_name = surface downwelling shortwave flux at current time + units = W m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[weasd] + standard_name = water_equivalent_accumulated_snow_depth_over_ocean + long_name = water equiv of acc snow depth over ocean + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lakedepth] + standard_name = lake_depth + long_name = lake depth + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[lake] + standard_name = flag_nonzero_lake_surface_fraction + long_name = flag indicating presence of some lake surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[xlat] + standard_name = latitude + long_name = latitude + units = radian + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[delt] + standard_name = time_step_for_dynamics + long_name = dynamics time step + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[zlvl] + standard_name = height_above_ground_at_lowest_model_layer + long_name = layer 1 height above ground (not MSL) + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[elev] + standard_name = orography + long_name = orography + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[wet] + standard_name = flag_nonzero_wet_surface_fraction + long_name = flag indicating presence of some ocean or lake surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[flag_iter] + standard_name = flag_for_iteration + long_name = flag for iteration + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F +[yearlen] + standard_name = number_of_days_in_year + long_name = number of days in a year + units = days + dimensions = () + type = integer + intent = in + optional = F +[julian] + standard_name = julian_day + long_name = julian day + units = days + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[imon] + standard_name = forecast_month + long_name = current forecast month + units = none + dimensions = () + type = integer + intent = in + optional = F +[snwdph] + standard_name = surface_snow_thickness_water_equivalent_over_ocean + long_name = water equivalent snow depth over ocean + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[hice] + standard_name = sea_ice_thickness + long_name = sea ice thickness + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tsurf] + standard_name = surface_skin_temperature_after_iteration_over_ocean + long_name = surface skin temperature after iteration over ocean + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[fice] + standard_name = sea_ice_concentration + long_name = ice fraction over open water + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[t_sfc] + standard_name = surface_skin_temperature_over_ocean_interstitial + long_name = surface skin temperature over ocean (temporary use as interstitial) + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_ocean + long_name = kinematic surface upward sensible heat flux over ocean + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[evap] + standard_name = kinematic_surface_upward_latent_heat_flux_over_ocean + long_name = kinematic surface upward latent heat flux over ocean + units = kg kg-1 m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[ustar] + standard_name = surface_friction_velocity_over_ocean + long_name = surface friction velocity over ocean + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qsfc] + standard_name = surface_specific_humidity_over_ocean + long_name = surface air saturation specific humidity over ocean + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[ch] + standard_name = surface_drag_coefficient_for_heat_and_moisture_in_air_over_ocean + long_name = surface exchange coeff heat & moisture over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cm] + standard_name = surface_drag_coefficient_for_momentum_in_air_over_ocean + long_name = surface exchange coeff for momentum over ocean + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[chh] + standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_ocean + long_name = thermal exchange coefficient over ocean + units = kg m-2 s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cmm] + standard_name = surface_drag_wind_speed_for_momentum_in_air_over_ocean + long_name = momentum exchange coefficient over ocean + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F diff --git a/physics/sfc_ocean.F b/physics/sfc_ocean.F index e21ddb3a7..33a1d3082 100644 --- a/physics/sfc_ocean.F +++ b/physics/sfc_ocean.F @@ -29,7 +29,7 @@ end subroutine sfc_ocean_finalize !! subroutine sfc_ocean_run & & ( im, cp, rd, eps, epsm1, hvap, rvrdm1, ps, t1, q1, & ! --- inputs - & tskin, cm, ch, prsl1, prslki, wet, wind, & + & tskin, cm, ch, prsl1, prslki, wet, lake, wind, & & flag_iter, & & qsurf, cmm, chh, gflux, evap, hflx, ep, & ! --- outputs & errmsg, errflg & @@ -102,7 +102,7 @@ subroutine sfc_ocean_run & real (kind=kind_phys), dimension(im), intent(in) :: ps, & & t1, q1, tskin, cm, ch, prsl1, prslki, wind - logical, dimension(im), intent(in) :: flag_iter, wet + logical, dimension(im), intent(in) :: flag_iter, wet, lake ! --- outputs: real (kind=kind_phys), dimension(im), intent(inout) :: qsurf, & @@ -138,6 +138,7 @@ subroutine sfc_ocean_run & ! rho is density, qss is sat. hum. at surface if ( flag(i) ) then + if(.not.lake(i)) then q0 = max( q1(i), 1.0e-8 ) rho = prsl1(i) / (rd*t1(i)*(1.0 + rvrdm1*q0)) @@ -166,6 +167,7 @@ subroutine sfc_ocean_run & hflx(i) = hflx(i) * tem * cpinv evap(i) = evap(i) * tem * hvapi endif + endif !end of if not lake enddo ! return diff --git a/physics/sfc_ocean.meta b/physics/sfc_ocean.meta index d60c1ce2c..733e69f54 100644 --- a/physics/sfc_ocean.meta +++ b/physics/sfc_ocean.meta @@ -153,6 +153,14 @@ type = logical intent = in optional = F +[lake] + standard_name = flag_nonzero_lake_surface_fraction + long_name = flag indicating presence of some lake surface area fraction + units = flag + dimensions = (horizontal_dimension) + type = logical + intent = in + optional = F [wind] standard_name = wind_speed_at_lowest_model_layer long_name = wind speed at lowest model level