1. Background diffusivity for shear-instability from KPP added 2. Alsta lake config refined

setup/alsta_driver.dat

@@ -218,11 +218,15 @@ lakeform  1
 trib_inflow    -9999.
 effl_outflow 1
 0. 0. 0.
-morphometry    5
-0	1.13E+6
-1	0.8E+6
-2       0.5E+6
-3       0.2E+6
+morphometry    9
+0	1.14E+6
+0.5 0.99E+6
+1	0.85E+6
+1.5 0.70E+6
+2   0.56E+6
+2.5 0.40E+6
+3   0.24E+6
+3.5 0.11E+6
 4       100.
 #
 #-----------------------------------------------------------------------------------------

setup/alsta_driver_new.dat

+# The parameters for driver of Lake model
+#-----------------------------------------------------------------------------------------
+#                     INFORMATION ON THE FILE WITH ATMOSPHERIC DATA
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# dataname*   --- name of file with atmospheric data (must be in data directory)
+# height_T_q* --- height of temperature and humidity measurements, m
+# height_u*   --- height of wind measurements, m
+# interval*   --- time interval of measurements, hours
+# rad*        --- defines, atmospheric radiation (1) or net radiation (2) is
+#                 in the appropriate coloumn of datafile, relevant if input file is ASCII (forc_format = 0)
+# forc_format*--- defines the input file format: 0 - ASCII(text), 1 - netcdf
+# npoints*    --- the number of points of the forcing, must be 1, if input file is ASCII (forc_format = 0)
+# select_call --- the length of the set of numbers of forcing points (maximal 20), for which the Lake model
+#                 will be launched, other points will be omitted
+# form*       --- defines the input file format, relevant if input file is ASCII (forc_format = 0):
+#                                                     0 - "free" (adjustable) format (see below)
+#                                                     other options are disabled
+#
+# The parameters for adjustable format of input text file, relevant if forc_format = 0, form = 0 :
+#
+# N_header_lines* --- the number of lines, occupied by file header 
+# N_coloumns*     --- the total number of coloumns in the file
+# N_Year*         --- the number of coloumn with the number of year (not used in the model)
+# N_Month*        --- the number of coloumn with the number of month (not used in the model)
+# N_Day*          --- the number of coloumn with the number of day (not used in the model)
+# N_Hour*         --- the number of coloumn with the number of hour (not used in the model)
+# N_Uspeed*       --- the number of coloumn with x-component speed values,      (m/s)
+# N_Vspeed*       --- the number of coloumn with y-component speed values,      (m/s)
+# N_Temp*         --- the number of coloumn with air temperature values,        (K)
+# N_Hum*          --- the number of coloumn with air humidity values,           (kg/kg)
+# N_Pres*         --- the number of coloumn with atmospheric pressure value,s   (Pa)
+# N_SWdown*       --- the number of coloumn with net solar radiation values,    (W/m**2)
+# N_LWdown*       --- the number of coloumn with net longwave radiation values, (W/m**2)
+# N_Precip*       --- the number of coloumn with precipitation intensity,       (m/s)
+#-----------------------------------------------------------------------------------------
+#
+dataname    'alsta'
+forc_format 0
+npoints     1
+# 8
+#select_call 3
+#1
+#2
+#3
+#
+lakinterac  1
+form        0
+height_T_q  2.5
+height_u    2.5
+interval    1
+rad         1
+#
+N_header_lines  1
+N_coloumns      8
+#
+N_Year          -1
+N_Month         -1
+N_Day           -1
+N_Hour          -1
+N_Uspeed        4
+N_Vspeed        5
+N_Temp          1
+N_Hum           2
+N_Pres          3
+N_SWdown        6
+N_LWdown        7
+N_Precip        8
+N_SensFlux      -1
+N_LatentFlux    -1
+N_Ustar         -1
+N_surfrad       -1
+N_cloud         -1
+N_NetRad        -1
+N_SurfTemp      -1
+#
+#-----------------------------------------------------------------------------------------
+#                              TIME INTEGRATION PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# year0*       --- julian year of start of integration
+# month0*      --- julian month of start of integration
+# day0*        --- julian day of start of integration
+# hour0*       --- hour of start of integration (is real value)
+# dt*          --- timestep,         s
+# tinteg*      --- integration time (including spinup period!), days
+# spinup_times*  --- number for spinup periods
+# spinup_period--- the duration of spinup period, s
+# call_Flake   --- the switch for integrating Flake model (1 - on, 0 - off)
+#-----------------------------------------------------------------------------------------
+#
+year0       2020
+month0      1
+day0        1
+hour0       0
+#
+tinteg      731
+spinup_times  0
+spinup_period 0
+cp_period     0
+control_point 0
+dt          10
+call_Flake  0
+#
+#----------------------------------------------------------------------------------------
+#                                   PHYSICAL PARAMETERS
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# extwat       --- coefficient of solar radiation extinction in water body, m**(-1)
+# extice       --- coefficient of solar radiation extinction in ice layer,  m**(-1)
+# alphax       --- slope angle of water surface in the x-direction,         deg
+# alphay       --- slope angle of water surface in the y-direction,         deg
+# c_veg        --- friction coefficient of vegetation in water,             n/d
+# a_veg        --- effective cross-section of vegetation,                   m**2/m**2
+# h_veg        --- the height of vegetation in the lake,                    m
+# kor          --- Coriolis parameter,                                      s**(-1)
+#                  if -999, it is calculated from latitude
+# phi*         --- latitude (positive to North), deg (required only for water albedo calculation)
+# lam*         --- longitude (positive to East), dag (required only for water albedo calculation)
+# fetch        --- the wind fetch,                                          m
+#----------------------------------------------------------------------------------------
+#
+extwat      3.4 #
+# 3.
+#select_extwat 9
+#2   0.02
+#20  0.12
+#108 0.14
+#86  0.17
+#176 0.12
+#105 0.24
+#39  0.17
+#170 0.14
+#163 1.2
+#
+extice      1.E+7
+alphax      0.0
+alphay      0.0
+a_veg       1.
+c_veg       1.e-3
+h_veg       0.
+#kor         0.
+kor         -999.
+# 1.e-4
+phi         59.74204 
+lam         17.252910
+fetch       1500.
+#
+#----------------------------------------------------------------------------------------
+#                                   INITIAL CONDITIONS
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# l10   --- thickness of ice,                          m
+# h10   --- thickness of liquid water,                 m
+# select_h10 --- thickness of liquid water in selected points (maximal 20), m
+#                (must be specified AFTER select_call)
+# ls10  --- thickness of ice at the bottom,            m
+# hs10  --- thickness of snow cover,                   m
+# Ts0   --- temperature of mixed layer,                C
+# Tb0   --- temperature at the bottom,                 C
+# Sals0 --- salinity in mixed layer,                   kg/kg
+# Salb0 --- salinity at the bottom,                    kg/kg
+# us0   --- x-component of speed at the surface,       m/s
+# vs0   --- y-component of speed at the surface,       m/s
+# Tbb0  --- temperature at the lower boundary of soil, C
+# h_ML0 --- thickness of mixed layer,                  m
+# init_T--- the type of temperature profile initialization:
+#           1 - using h_ML0, Ts0 and Tb0
+#           2 - using Tm,    Ts0 and Tb0
+#           3 - using the temperature profile, specified in *_setup.dat file after keyword T_profile
+#-----------------------------------------------------------------------------------------
+#
+l10         0.
+h10         4.
+#h10         3.
+#select_h10  3
+#1  1.2
+#2  2.
+#3  5.
+#select_h10  8
+#1  3.
+#2  4.
+#3  5.
+#4  6.
+#5  7.
+#6  8.
+#7  9.
+#8  10.
+#
+ls10        0.
+hs10        0.
+Ts0         15.
+Tb0         5.
+Tbb0        5.
+Tm          3.
+h_ML0       1.5
+Sals0       1.E-3
+Salb0       30.E-3
+us0         1.E-3
+vs0         1.E-3
+init_T      3
+#
+#-----------------------------------------------------------------------------------------
+#                              SOME LAKE PARAMETERS
+#-----------------------------------------------------------------------------------------
+# area_lake     --- the area of the lake,                              m**2
+# trib_inflow*  --- total tributaries' inflow,                         m**3/s
+# morphometry   --- depth - lake cross-section area table,             m**2
+# effl_outflow  --- effluent discharge parameters group, specifying polynomial dependence
+#                   of discharge on water level; the last value is relative altitude of
+#                   effluent bottom over lake bottom (at deepest points, respectively);
+#                   the first N values are polynomial coefficients, where N stands after
+#                   'effl_outflow' keyword
+#
+area_lake      1.13E+6
+cellipt 2.5
+lakeform  1
+trib_inflow    -9999.
+effl_outflow 1
+0. 0. 0.
+morphometry    9
+0	1.14E+6
+0.5 0.99E+6
+1	0.85E+6
+1.5 0.70E+6
+2   0.56E+6
+2.5 0.40E+6
+3   0.24E+6
+3.5 0.11E+6
+4       100.
+#
+#-----------------------------------------------------------------------------------------
+#                               NETCDF OUTPUT PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# nstep_ncout  --- the interval of netcdf output from driver, timesteps (if -1 no netcdf output from driver)
+#-----------------------------------------------------------------------------------------
+nstep_ncout -1
+#15
+#
+#-----------------------------------------------------------------------------------------
+#                               FLAKE MODEL OUTPUT PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# nstep_out_Flake  --- the interval of output of Flake variables from driver, timesteps 
+#                   (if -1 the output of Flake variables from driver is not implemented)
+#                   relevant if call_Flake = 1
+#-----------------------------------------------------------------------------------------
+#
+nstep_out_Flake 3
+#
+#-----------------------------------------------------------------------------------------
+#                              POSTPROCESSING OPTIONS
+#-----------------------------------------------------------------------------------------
+# moving_average_window --- the ineterval of moving average, netcdf output steps (intervals)
+# mean_cycle_period     --- the length of mean cycle, netcdf output steps (intervals)
+#-----------------------------------------------------------------------------------------
+#
+moving_average_window -1
+mean_cycle_period     -1
+#
+end

setup/alsta_setup.dat

@@ -149,7 +149,7 @@ dyn_pgrad 0
 pgrad    0.
 nManning 5.E-2
 horvisc  0.
-backdiff 0
+backdiff 2
 backdiff0 -999.
 botfric  1
 zero_model 0

source/model/lake.f90

@@ -222,7 +222,7 @@ use PHYS_FUNC, only: &
 & WL_MAX, &
 & MIXED_LAYER_CALC, &
 & HC_CORR_CARBEQUIL, &
-& DIFFMIN_HS, &
+& DIFFMIN_HS, DIFFMIN_KPP, &
 & LONGWAVE_RADIATION, &
 & SHORTWAVE_RADIATION
 
@@ -1072,7 +1072,14 @@ select case(backdiff%par)
 case(0)
   lamw_back(:) = 0.
 case(1)
-  call DIFFMIN_HS(area_lake,wst,gsp,ls,M) !Calculating background heat conductance
+  !Background heat conductivity following (Fang and Stefan, 1996)
+  call DIFFMIN_HS(area_lake,wst,gsp,ls,M) 
+case(2)
+  !Background heat conductivity from KPP model (e.g., Zhang et al., HESS, 2019)
+  call DIFFMIN_KPP(wst,gsp,ls,M) 
+case(3)
+  !Constant background vertical diffusivity 
+  lamw_back(:) = max(backdiff0%par,0._ireals) 
 end select
 
 ! Heat conductance at the next timestep

source/model/phys_func.f90

@@ -2082,6 +2082,57 @@ END FUNCTION FP_THETA
   if (firstcall) firstcall = .false.
   END SUBROUTINE DIFFMIN_HS
 
+
+  !> Subroutine calculates background eddy diffusivity from KPP model 
+  !! (Zhang et al., HESS, 2019 and references therein)
+  SUBROUTINE DIFFMIN_KPP(wst,gsp,ls,M) 
+  use PHYS_CONSTANTS, only : row0, g, cw_m_row0
+  use ARRAYS_GRID, only : gridspacing_type
+  use ARRAYS_WATERSTATE, only : waterstate_type
+  use ARRAYS_BATHYM, only : layers_type
+  use DRIVING_PARAMS, only : backdiff0
+  implicit none
+  !Input/output variables
+  type(waterstate_type) , intent(inout) :: wst
+  type(gridspacing_type), intent(in) :: gsp
+  type(layers_type)     , intent(in) :: ls
+  integer(kind=iintegers), intent(in) :: M
+  !Local variables
+  logical, save :: firstcall = .true.
+  real(kind=ireals), parameter :: k0_default = 1.E-5, Ri0 = 7.E-1, p = 3.
+  real(kind=ireals), save :: k0
+  real(kind=ireals) :: x, y, Ri, dz
+  integer(kind=iintegers) :: i !Loop index
+
+  if (firstcall) then
+    if (backdiff0%par >= 0.) then
+      k0 = backdiff0%par 
+    else
+      k0 = k0_default !default value
+    endif
+  endif
+
+  do i = 1, M
+    dz = ls%h1*gsp%ddz(i)
+    !Brunt-Vaisala frequency squared
+    x = g/row0*(wst%row(i+1) - wst%row(i)) / dz
+    !Shear frequency squared
+    y = ( (wst%u1(i+1) - wst%u1(i))**2 + (wst%v1(i+1) - wst%v1(i))**2 ) / dz**2
+    !Gradient Richardson number
+    Ri = x / y
+    if     (Ri <= 0.) then
+      wst%lamw_back(i) = k0
+    elseif (Ri < Ri0) then
+      wst%lamw_back(i) = k0*(1. - (Ri/Ri0)**2)**p
+    else
+      wst%lamw_back(i) = 0.
+    endif
+  enddo
+  wst%lamw_back(:) = wst%lamw_back(:) * cw_m_row0
+
+  if (firstcall) firstcall = .false.
+  END SUBROUTINE DIFFMIN_KPP
+
   !> Function HORIZCONC returns a value of dissolved concentration in the mixed layer
   !! at distance r from center of circular lake; according to 
   !! (DelSontro et al., Ecosystems, 2018) extended by inclusion of linear sink term