diff --git a/setup/alsta_driver.dat b/setup/alsta_driver.dat index d45bce80c0a3a05c3d2c8c4088422e3379801bd8..479e75bbf20a14b2a01f1ed68380d4219ad54be2 100644 --- a/setup/alsta_driver.dat +++ b/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. # #----------------------------------------------------------------------------------------- diff --git a/setup/alsta_driver_new.dat b/setup/alsta_driver_new.dat new file mode 100644 index 0000000000000000000000000000000000000000..479e75bbf20a14b2a01f1ed68380d4219ad54be2 --- /dev/null +++ b/setup/alsta_driver_new.dat @@ -0,0 +1,262 @@ +# 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 diff --git a/setup/alsta_setup.dat b/setup/alsta_setup.dat index cf02db3c10464dc833b2e7b07348906d8d78de83..d90832ce98b80e51720a9aa58e7b44e4e70346a1 100644 --- a/setup/alsta_setup.dat +++ b/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 diff --git a/source/model/lake.f90 b/source/model/lake.f90 index 2e5e9d8aa0ab81b03588e4fc1a78d55d7d22a1da..c86a38409aff9cacafe80b16ef3832a98bc95590 100644 --- a/source/model/lake.f90 +++ b/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 diff --git a/source/model/phys_func.f90 b/source/model/phys_func.f90 index eab6a6e5a0c3ca02f22cf7da95395adf69f74c54..58d6cd48f48fb853b660af5ee941a3edc36b2b3b 100644 --- a/source/model/phys_func.f90 +++ b/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