From 2d71e7ae8f32cefa919fa7d8b44a3e76fed60341 Mon Sep 17 00:00:00 2001 From: Fabian Bernhard Date: Tue, 9 Mar 2021 17:59:53 +0100 Subject: [PATCH] Update documentation (#8) * Update documentation * Fix some typos * Fix further typos * Increment version number --- Project.toml | 2 +- README.md | 6 +-- docs/src/index.md | 6 +-- docs/src/model.md | 3 +- docs/src/user-guide.md | 5 -- src/func_DiffEq_definition_cb.jl | 29 ++++------- src/func_DiffEq_definition_f.jl | 6 ++- src/func_DiffEq_definition_ode.jl | 17 ++++++ src/func_DiffEq_definition_p.jl | 21 ++++++-- src/func_DiffEq_definition_u0.jl | 7 ++- src/func_MSB_functions.jl | 10 +++- src/func_input_definition.jl | 42 ++++++++++----- src/module_EVP.jl | 84 +++++++++++++++++++++++------- src/module_GLBLDECL.jl | 32 ++++++++---- src/module_KPT.jl | 61 +++++++++++++--------- src/module_PET.jl | 81 +++++++++++++++++++++-------- src/module_SNO.jl | 47 ++++++++++++----- src/module_SUN.jl | 27 ++++++++-- src/module_WAT.jl | 86 ++++++++++++++++++++----------- 19 files changed, 396 insertions(+), 176 deletions(-) diff --git a/Project.toml b/Project.toml index 0eb9d694..f8653b2b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LWFBrook90" uuid = "100deb6c-4418-4b77-ad5c-a1fadb580743" authors = ["Fabian Bernhard"] -version = "0.1.2" +version = "0.1.3" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/README.md b/README.md index 55aaa6ff..df4470de 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ Implementation of the LWF-BROOK90 hydrological model in Julia

## What is LWFBrook90.jl? -The model LWF-BROOK90 is a 1D Soil Vegetation Atmosphere Transport (SVAT) model, calculating the soil water balance in forest soil. Modelled processes include vertical soil water movement, soil and plant evapotranspiration and temporary storages in snowpack or interception layer. +The model LWF-BROOK90 is a 1D Soil-Vegetation-Atmosphere Transfer (SVAT) model, calculating the soil water balance in forest soil. Modelled processes include vertical soil water movement, soil and plant evapotranspiration and temporary storages in snowpack or interception layer. LWFBrook90.jl is an implementation of the existing LWF-BROOK90 rewritten entirely in the Julia programming language. Processes and state variables of the model in LWF-BROOK90 are summarised visually in Figure 1 below. @@ -41,7 +41,7 @@ Development of LWFBrook90.jl builds on the following works: Furthermore, Matthias Häni, Katrin Meusburger, Peter Waldner, Lorenz Walthert, Stephan Zimmermann of [WSL](http://www.wsl.ch) and its Long-term Forest Ecosystem Research (LWF) project gratefully acknowledged for providing example data files located in `example/BEA2016*`. -LWFBrook90.jl is licensed under GPL-3.0 (see the file `LICENSE`). +LWFBrook90.jl is licensed under GPL-3.0 (for details see the file `LICENSE`). @@ -76,7 +76,7 @@ Tests are run to assert agreement with results from LWFBrook90R. Visualizations - Currently a part of LWFBrook90R activated when `Reset==1`, is not implemented in LWFBrook90.jl. (A test implementation of this is available in the branch `005-simplify-time-step-control`.) ### Improve agreement with LWFBrook90R -Note that some features of LWFBrook90R are not implemented in the main version of LWFBrook90.jl. The time step adaptivity and `Reset==1` of LWFBrook90R would require code refactoring, which goes slightly against the intended use of the library for ODEs DiffEqeuations.jl. Because of this, implementation of these features is left aways from the main version. However, an attempt at their implementation resides currently in a feature branch `005`. Below are some of the results of that branch: +Note that some features of LWFBrook90R are not implemented in the main version of LWFBrook90.jl. The time step adaptivity and `Reset==1` of LWFBrook90R would require code refactoring, which goes slightly against the intended use of the library for ODEs DifferentialEquations.jl. Because of this, implementation of these features is left away from the main version. However, an attempt at their implementation resides currently in a feature branch `005`. Below are some of the results generated with that branch:

Figure 6: Comparing daily outputs of LWFBrook90R and experimental LWFBrook90.jl:feature-005 for example data set over a year

diff --git a/docs/src/index.md b/docs/src/index.md index a2ae67dd..b1e1d789 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -14,7 +14,7 @@ Depth = 2 ## About LWFBrook90.jl -LWFBrook90.jl implements a 1D soil vegetation atmosphere transport model. +LWFBrook90.jl implements a 1D soil-vegetation-atmosphere transfer model. Intended use cases of this impelemntation are: - support of stable isotopes (δ¹⁸O and δ²H) - efficient model calibration for data analysis @@ -28,7 +28,7 @@ For further details read through the section [User Guide](@ref) or refer to sect ## Citing LWFBrook90.jl When using LWFBrook90.jl please cite ->TODO: generate citation: (Journal of Open Source Software? Zenodo? Journal of Open Research Software (DiffEq.jl)? Alternatives?) +>TODO: generate citation: (Journal of Open Source Software? Zenodo? Journal of Open Research Software (DifferentialEquations.jl)? Alternatives?) @@ -46,6 +46,4 @@ Federer, C. A. (2002). BROOK 90: A simulation model output for evaporation, soil Hammel, K., & Kennel, M. (2001). Charakterisierung und Analyse der Wasserverfügbarkeit und des Wasserhaushalts von Waldstandorten in Bayern mit dem Simulationsmodell BROOK90 (No. 185; *Forstliche Forschungsberichte München*, p. 135). Technische Uni München Wissenschaftszentrum Weihenstephan. ISBN 3-933506-16-6 -Hammel, K., & Kennel, M. (2001). Charakterisierung und Analyse der Wasserverfügbarkeit und des Wasserhaushalts von Waldstandorten in Bayern mit dem Simulationsmodell BROOK90 (No. 185; *Forstliche Forschungsberichte München*, p. 135). Technische Uni München Wissenschaftszentrum Weihenstephan. ISBN 3-933506-16-6 - Schmidt-Walter, P., Trotsiuk, V., Meusburger, K., Zacios, M., & Meesenburg, H. (2020). Advancing simulations of water fluxes, soil moisture and drought stress by using the LWF-Brook90 hydrological model in R. *Agricultural and Forest Meteorology*, 291, 108023. https://doi.org/10.1016/j.agrformet.2020.108023 \ No newline at end of file diff --git a/docs/src/model.md b/docs/src/model.md index 818c44cf..7a1e5a21 100644 --- a/docs/src/model.md +++ b/docs/src/model.md @@ -1,7 +1,7 @@ # SVAT Model ## Description -LWFBrook90.jl is a 1D Soil Vegetation Atmosphere Transport (SVAT) model, calculating the soil water balance in forest soil. Modelled processes include vertical soil water movement, soil and plant evapotranspiration and temporary storages in snowpack or interception layer. +LWFBrook90.jl is a 1D Soil-Vegetation-Atmosphere Transfer (SVAT) model, calculating the soil water balance in forest soil. Modelled processes include vertical soil water movement, soil and plant evapotranspiration and temporary storages in snowpack or interception layer. Vertical soil water movement is modelled using the Richards equations and preferential flow. Mass loss through evaporation from temporary storages (snowpack or interception by vegetation) is included. @@ -14,7 +14,6 @@ Processes and state variables in LWF-BROOK90 are summarised visually in Figure 1

``` - ## Implementation The model is implemented based code from the R package LWFBrook90R and its Fortran source code as well as the original implementation of BROOK90 (v4.8) (www.ecoshift.net/brook/b90doc.html). diff --git a/docs/src/user-guide.md b/docs/src/user-guide.md index 82a6d22a..e3dbf207 100644 --- a/docs/src/user-guide.md +++ b/docs/src/user-guide.md @@ -11,11 +11,6 @@ To install LWFBrook90.jl open a Julia REPL, enter the Pkg REPL by pressing `]` a ``` (@v1.5) pkg> add LWFBrook90 ``` -If the you want to use the develop version not in the package registry you can install directly from github using the command: -``` -(@v1.5) pkg> add https://github.com/fabern/LWFBrook90.jl -``` - Dependencies of LWFBrookJulia.jl should automatically be installed. ### Usage diff --git a/src/func_DiffEq_definition_cb.jl b/src/func_DiffEq_definition_cb.jl index a9fa921f..5c373893 100644 --- a/src/func_DiffEq_definition_cb.jl +++ b/src/func_DiffEq_definition_cb.jl @@ -1,8 +1,11 @@ -""" define_LWFB90_cb()\n - Generates callback function cb needed for ODE() probelm in DiffEq.jl package. - LWFBrook90 updates states INTS, INTR, SNOW, CC, SNOWLQ not continuously but only - once per day. This operator splitting (daily vs continuous update of ODEs) is - implemented by using this callback function which is called once per day. +""" + define_LWFB90_cb() + +Generate callback function cb needed for ODE() problem in DiffEq.jl package. + +LWFBrook90 updates states INTS, INTR, SNOW, CC, SNOWLQ not continuously but only +once per day. This operator splitting (daily vs continuous update of ODEs) is +implemented by using this callback function which is called once per day. """ function define_LWFB90_cb() # A) Define updating function @@ -271,18 +274,4 @@ function define_LWFB90_cb() cb_func = PeriodicCallback(LWFBrook90R_update_INTS_INTR_SNOW_CC_SNOWLQ!, 1.0; initial_affect = true); return cb_func -end - -""" define_LWFB90_ODE()\n - Generates an ODEProblem from DiffEq.jl.\n\n - # An ODE problem which consists of - # - definition of right-hand-side (RHS) function f - # - definition of callback function cb - # - initial condition of states - # - definition of simulation time span - # - parameters\n - Seperate updating of different states (INTS, INTR, SNOW, CC, SNOWLQ are updated once per - day while GWAT and SWATI are updated continuously) is implemented by means of operator - splitting using a callback function for the daily updates and a ODE RHS (right hand - side) for the continuous update. -""" \ No newline at end of file +end \ No newline at end of file diff --git a/src/func_DiffEq_definition_f.jl b/src/func_DiffEq_definition_f.jl index dca616a2..7c0e9cc2 100644 --- a/src/func_DiffEq_definition_f.jl +++ b/src/func_DiffEq_definition_f.jl @@ -1,5 +1,7 @@ -""" define_LWFB90_f()\n - Generates function f (right-hand-side of ODEs) needed for ODE() probelm in DiffEq.jl package. +""" + define_LWFB90_f() + +Generate function f (right-hand-side of ODEs) needed for ODE() problem in DiffEq.jl package. """ # function define_LWFB90_f() function f_LWFBrook90R(du,u,p,t) diff --git a/src/func_DiffEq_definition_ode.jl b/src/func_DiffEq_definition_ode.jl index 22d772a1..5640a77e 100644 --- a/src/func_DiffEq_definition_ode.jl +++ b/src/func_DiffEq_definition_ode.jl @@ -1,3 +1,20 @@ +""" + define_LWFB90_ODE() + +Generates an ODEProblem from DiffEq.jl + +An ODE problem which consists of + - definition of right-hand-side (RHS) function f + - definition of callback function cb + - initial condition of states + - definition of simulation time span + - parameters + +Seperate updating of different states (INTS, INTR, SNOW, CC, SNOWLQ are updated once per +day while GWAT and SWATI are updated continuously) is implemented by means of operator +splitting using a callback function for the daily updates and a ODE RHS (right hand +side) for the continuous update. +""" function define_LWFB90_ODE(u0, tspan, p) # Define callback functions diff --git a/src/func_DiffEq_definition_p.jl b/src/func_DiffEq_definition_p.jl index f6941d50..437ee593 100644 --- a/src/func_DiffEq_definition_p.jl +++ b/src/func_DiffEq_definition_p.jl @@ -1,6 +1,21 @@ -""" define_diff_eq_parameters(NLAYER, IMODEL, constant_dt_solver, NOOUTF, Reset, compute_intermediate_quantities, - pfile_meteo, pfile_siteparam, pfile_param, pfile_soil, pfile_pdur)\n - Generates vector p needed for ODE() problem in DiffEq.jl package.""" +""" + define_diff_eq_parameters() + +Generate vector p needed for ODE() problem in DiffEq.jl package. + +# Arguments +- `NLAYER::...`: TODO argument description. +- `IMODEL::...`: TODO argument description. +- `constant_dt_solver::...`: TODO argument description. +- `NOOUTF::...`: TODO argument description. +- `Reset::...`: TODO argument description. +- `compute_intermediate_quantities::...`: TODO argument description. +- `pfile_meteo::...`: TODO argument description. +- `pfile_siteparam::...`: TODO argument description. +- `pfile_param::...`: TODO argument description. +- `pfile_soil::...`: TODO argument description. +- `pfile_pdur::...`: TODO argument description. +""" function define_LWFB90_p(NLAYER, IMODEL, constant_dt_solver, NOOUTF, Reset, compute_intermediate_quantities, pfile_meteo, pfile_siteparam, pfile_param, pfile_soil, pfile_pdur) diff --git a/src/func_DiffEq_definition_u0.jl b/src/func_DiffEq_definition_u0.jl index d8a3af63..3b088dd5 100644 --- a/src/func_DiffEq_definition_u0.jl +++ b/src/func_DiffEq_definition_u0.jl @@ -1,5 +1,8 @@ -""" define_LWFB90_u0()\n - Generates vector u0 needed for ODE() problem in DiffEq.jl package.""" +""" + define_LWFB90_u0() + +Generate vector u0 needed for ODE() problem in DiffEq.jl package. +""" function define_LWFB90_u0(u_GWAT_init, u_INTS_init, u_INTR_init, diff --git a/src/func_MSB_functions.jl b/src/func_MSB_functions.jl index 09daa92c..a5271281 100644 --- a/src/func_MSB_functions.jl +++ b/src/func_MSB_functions.jl @@ -1,6 +1,12 @@ # TODO(bernhard): think about where to put function definition of MSBSETVARS(), MSBDAYNIGHT() -"""MSBSETVARS() function that computes state dependent parameters for -updating states INTS, INTR, SNOW, CC, SNOWLQ in callback function +""" + MSBSETVARS() + +Compute state dependent parameters for updating states INTS, INTR, SNOW, CC, SNOWLQ in +callback function. + +# Arguments +- many """ function MSBSETVARS(IDAY, #TODO(bernhard) just for debug... remove again! # arguments diff --git a/src/func_input_definition.jl b/src/func_input_definition.jl index 92ac69f6..a52c6658 100644 --- a/src/func_input_definition.jl +++ b/src/func_input_definition.jl @@ -6,17 +6,18 @@ using DataFramesMeta#: @linq, transform, DataFramesMeta using Dates: DateTime, Millisecond, Second, Day, month, value, dayofyear """ -read_LWFBrook90R_inputData(folder::String, prefix::String)\n -Intent: load different input files for LWFBrook90\n + read_LWFBrook90R_inputData(folder::String, prefix::String) + +Load different input files for LWFBrook90: - climveg - param - siteparam - pdur - soil_materials.csv - soil_nodes.csv -\n\n + These files were created with an R script `generate_LWFBrook90Julia_Input.R` that -takes the same arguements as the R funciton LWFBrook90R::run_LWFB90() and generates +takes the same arguements as the R funciton `LWFBrook90R::run_LWFB90()` and generates the corresponding Julia input functions. """ function read_LWFBrook90R_inputData(folder::String, prefix::String) @@ -117,31 +118,48 @@ end ###################### # Define functions to handle DateTimes and convert into Days as Floats -""" DateTime2RelativeDaysFloat(x, - reference_DateTime)\n Transforms DateTimes `x` to simulation time """ +""" + DateTime2RelativeDaysFloat(x,reference_DateTime) + +Transforms DateTimes `x` to simulation time +""" function DateTime2RelativeDaysFloat(x::DateTime, reference::DateTime) ms2days = 1.0/(24.0*3600.0*1000.0) # to convert milliseconds to days ms2days*value(convert(Millisecond, x-reference)) end -""" RelativeDaysFloat2DateTime(t, - reference_DateTime)\n Transforms simulation time t to DateTimes `dt` """ +""" + RelativeDaysFloat2DateTime(t, reference_DateTime) + +Transforms simulation time `t` to DateTimes +""" function RelativeDaysFloat2DateTime(t::Float64, reference::DateTime) # reference + Day(floor(t)) t_sec = 60*60*24*t # t is in days, t_sec in seconds reference + Second(floor(t_sec)) end -""" p_DOY(t::Float64, reference::DateTime)\n Get DOY (Day Of Year) from simulation time""" +""" + p_DOY(t::Float64, reference::DateTime) + +Get DOY (Day Of Year) from simulation time +""" function p_DOY(t::Float64, reference::DateTime) dayofyear(reference + Day(floor(t))) end -""" p_MONTHN(t::Float64, reference::DateTime)\n Get Month from simulation time""" +""" + p_MONTHN(t::Float64, reference::DateTime) + +Get Month from simulation time +""" function p_MONTHN(t::Float64, reference::DateTime) month(reference + Day(floor(t))) end # Subset input data and transform dates into floats relative to reference_date -# """ subset_standardize(data::DataFrame, start::DateTime, stop::DateTime, reference::DateTime)\n -# Returns DataFrame `data` that is subset between `start` and `stop` and colum `dates` transformed to simulation time.""" +# """ +# subset_standardize(data::DataFrame, start::DateTime, stop::DateTime, reference::DateTime) +# +# Returns DataFrame `data` that is subset between `start` and `stop` and colum `dates` transformed to simulation time. +# """ # function subset_standardize(data::DataFrame, start::DateTime, stop::DateTime, reference::DateTime) # @linq data |> # # Subset diff --git a/src/module_EVP.jl b/src/module_EVP.jl index b259407e..c3fb1b28 100644 --- a/src/module_EVP.jl +++ b/src/module_EVP.jl @@ -24,7 +24,11 @@ using ..CONSTANTS: p_RHOWG, p_PI # https://discourse.julialang.org/t/large-progr export PLNTRES, TBYLAYER, INTER, INTER24 -""" PLNTRES() allocates total plant resistance to xylem and root layers +""" + PLNTRES() + +Allocates total plant resistance to xylem and root layers. + Ecoshift: Subroutine PLNTRES is called at the beginning of each day to obtain resistivities to liquid water flow: rhizosphere resistivity for each soil layer, root resistivity in each soil @@ -172,7 +176,11 @@ function PLNTRES(NLAYER, p_THICK, p_STONEF, p_fu_RTLEN, p_fT_RELDEN, p_RTRAD, p_ end -"""TBYLAYER() allocates transporation among soil layers. +""" + + TBYLAYER() + +Allocate transporation among soil layers. Ecoshift: TBYLAYER determines the rate at which liquid water can be supplied to transpiring leaves, @@ -388,8 +396,26 @@ function TBYLAYER(J, p_fu_PTR, p_fu_DISPC, p_fu_ALPHA, p_fu_KK, p_fu_RROOTI, p_f end -"""Rain interception, used when p_NPINT > 1 +""" + INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_CINTRL, p_CINTRS, p_DTP, u_INTR) + +Compute rain catch rate (interception) and evaporation rate of intercepted rain in mm/d. +Rain interception, used when p_NPINT > 1. + +# Arguments +- `p_fT_RFAL`: rainfall rate, mm/d +- `p_fu_PINT`: potential interception rate, mm/d +- `p_fu_LAI`: projected leaf area index, m2/m2 +- `p_fu_SAI`: projected stem area index, m2/m2 +- `p_FRINTL`: intercepted fraction of p_fT_RFAL per unit p_fu_LAI +- `p_FRINTS`: intercepted fraction of p_fT_RFAL per unit p_fu_SAI +- `p_CINTRL`: maximum interception storage of rain per unit p_fu_LAI, mm +- `p_CINTRS`: maximum interception storage of rain per unit p_fu_SAI, mm +- `p_DTP`: precipitation interval time step, d +- `u_INTR`: intercepted rain storage, mm, + +# Ecoshift: Older studies of rain and snow interception regressed throughfall on precipitation, but such interpretation ignored the fact that energy supply rather than water supply may limit interception and also ignores storm duration/intensity and interstorm interval. Detailed @@ -482,7 +508,7 @@ function INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_C # p_FRINTS intercepted fraction of p_fT_RFAL per unit p_fu_SAI # p_CINTRL maximum interception storage of rain per unit p_fu_LAI, mm # p_CINTRS maximum interception storage of rain per unit p_fu_SAI, mm - # DTP precipitation interval time step, d + # p_DTP precipitation interval time step, d # u_INTR intercepted rain storage, mm, # maximum RINT, mm/d @@ -516,10 +542,29 @@ function INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_C return (aux_du_RINT, aux_du_IRVP) end -""" Rain interception with duration in hours, used when p_NPINT = 1. Same routine is used for -snow interception, with different calling variables - -Ecoshift: +""" + INTER(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_CINTRL, p_CINTRS, p_DTP, u_INTR) + +Compute rain catch rate (interception) and evaporation rate of intercepted rain in mm/d. + +Rain interception with duration in hours, used when p_NPINT = 1. Same routine is used for +snow interception, with different calling variables. + +# Arguments: +- `p_fT_RFAL`: 24-hour average rainfall rate, mm/d +- `p_fu_PINT`: potential interception rate, mm/d +- `p_fu_LAI`: projected leaf area index, m2/m2 +- `p_fu_SAI`: projected stem area index, m2/m2 +- `p_FRINTL`: intercepted fraction of p_fT_RFAL per unit p_fu_LAI +- `p_FRINTS`: intercepted fraction of p_fT_RFAL per unit p_fu_SAI +- `p_CINTRL`: maximum interception storage of rain per unit p_fu_LAI, mm +- `p_CINTRS`: maximum interception storage of rain per unit p_fu_SAI, mm +- `p_DURATN`: average storm duration, hr +- `u_INTR`: intercepted rain storage, mm, +- `MONTHN`: Month of the year + +# Ecoshift: +" Subroutine INTER24 - daily interception Proper representation and integration of the interception process is a problem for hydrologic models that use a daily interval for precipitation input (p_NPINT = 1), because the @@ -546,19 +591,20 @@ To determine appropriate values of DURATN I examined hourly precipitation data f of precipitation of 0.02 inch (0.5 mm) or greater over days with such precipitation gave the following results after a little smoothing - J F M A M J J A S O N D -San Juan PR 3 2 2 2 2 2 2 3 3 3 3 3 -Atlanta GA 5 5 5 5 4 4 3 3 4 4 5 6 -Caribou ME 4 4 5 5 4 4 4 4 4 6 6 5 -Madison WI 4 4 5 3 3 2 3 3 4 4 5 5 -Lake Charles LA 5 4 3 3 3 3 2 2 3 3 4 5 -Phoenix AZ 4 4 4 4 4 2 2 2 2 2 4 4 -Rapid City SD 3 3 3 4 4 3 2 2 2 2 4 4 -Tacoma WA 6 6 5 4 4 4 4 4 4 4 6 6 -Fairbanks AK 3 3 4 4 4 4 3 3 4 4 4 3 -Hubbard Brook NH 5 5 5 4 4 4 4 4 4 5 5 5 + J F M A M J J A S O N D + San Juan PR 3 2 2 2 2 2 2 3 3 3 3 3 + Atlanta GA 5 5 5 5 4 4 3 3 4 4 5 6 + Caribou ME 4 4 5 5 4 4 4 4 4 6 6 5 + Madison WI 4 4 5 3 3 2 3 3 4 4 5 5 + Lake Charles LA 5 4 3 3 3 3 2 2 3 3 4 5 + Phoenix AZ 4 4 4 4 4 2 2 2 2 2 4 4 + Rapid City SD 3 3 3 4 4 3 2 2 2 2 4 4 + Tacoma WA 6 6 5 4 4 4 4 4 4 4 6 6 + Fairbanks AK 3 3 4 4 4 4 3 3 4 4 4 3 + Hubbard Brook NH 5 5 5 4 4 4 4 4 4 5 5 5 Apparently a default DURATN of 4 hours is appropriate for all months anywhere in the U.S. +" """ function INTER24(p_fT_RFAL, p_fu_PINT, p_fu_LAI, p_fu_SAI, p_FRINTL, p_FRINTS, p_CINTRL, p_CINTRS, p_DURATN, u_INTR, MONTHN) # p_fT_RFAL 24-hour average rainfall rate, mm/d diff --git a/src/module_GLBLDECL.jl b/src/module_GLBLDECL.jl index 35f582ef..19955a53 100644 --- a/src/module_GLBLDECL.jl +++ b/src/module_GLBLDECL.jl @@ -15,9 +15,10 @@ export derive_params_from_inputData # Exported functions """ -derive_params_from_inputData(input_meteo::DataFrame, input_param::DataFrame, input_siteparam::DataFrame, -input_precdat::DataFrame, input_pdur::DataFrame, input_soil_materials::DataFrame, input_soil_nodes::DataFrame)\n -Takes input data sets and defines parameters needed for simulation. + derive_params_from_inputData(input_meteo::DataFrame, input_param::DataFrame, input_siteparam::DataFrame, + input_precdat::DataFrame, input_pdur::DataFrame, input_soil_materials::DataFrame, input_soil_nodes::DataFrame) + +Take input data sets and defines parameters needed for simulation. """ function derive_params_from_inputData(input_meteo::DataFrame, input_param::DataFrame, @@ -59,8 +60,9 @@ end ####################### # Unexported functions """ -derive_params_from_input_meteo(input_meteo::DataFrame)\n -Takes climate and vegetation parameters in `input_meteo` and generates continuous parameters. + derive_params_from_input_meteo(input_meteo::DataFrame) + +Take climate and vegetation parameters in `input_meteo` and generates continuous parameters. """ function derive_params_from_input_meteo(input_meteo::DataFrame, input_reference_date::DateTime) @@ -107,7 +109,10 @@ function derive_params_from_input_meteo(input_meteo::DataFrame, input_reference_ ("input_reference_date", input_reference_date)]) end -"""derive_params_from_input_pdur(input_pdur)\n Function that defines constant parameters from input_pdur. TO BE REDEFINED +""" + derive_params_from_input_pdur(input_pdur) + +Define constant parameters from input_pdur. TO BE REDEFINED """ function derive_params_from_input_pdur(input_pdur) # Interception parameters ------- @@ -115,7 +120,10 @@ function derive_params_from_input_pdur(input_pdur) return Dict([("DURATN",DURATN)]) end -"""derive_params_from_input_siteparam(input_siteparam)\n Function that defines constant parameters from input_siteparam. TO BE REDEFINED +""" + derive_params_from_input_siteparam(input_siteparam) + +Define constant parameters from input_siteparam. TO BE REDEFINED """ function derive_params_from_input_siteparam(input_siteparam) p_LAT = input_siteparam[1,3]/57.296 @@ -130,7 +138,10 @@ function derive_params_from_input_siteparam(input_siteparam) ("p_NPINT",p_NPINT)]) end -"""derive_params_from_input_soil(input_soil_materials, input_soil_nodes, IMODEL, ILAYER, QLAYER, NLAYER)\n Function that defines constant parameters from input_pdur. TO BE REDEFINED +""" + derive_params_from_input_soil(input_soil_materials,IMODEL, ILAYER, QLAYER, NLAYER) + +Define constant parameters from input_pdur. TO BE REDEFINED """ function derive_params_from_input_soil(input_soil_materials, input_soil_nodes, IMODEL, ILAYER, QLAYER, NLAYER, HEAT, nmat, inirdep, rgrorate) # Parse soil parameter according to material at different depth given by @@ -322,7 +333,10 @@ function derive_params_from_input_soil(input_soil_materials, input_soil_nodes, I ("TopInfT", TopInfT)]) end -"""derive_params_from_input_param(input_param)\nFunction that defines constant parameters from input_param. TO BE REDEFINED +""" + derive_params_from_input_param(input_param) + +unction that defines constant parameters from input_param. TO BE REDEFINED """ function derive_params_from_input_param(input_param) # from LWFBrook90R:PFILE.h diff --git a/src/module_KPT.jl b/src/module_KPT.jl index 45408531..20d3a4f0 100644 --- a/src/module_KPT.jl +++ b/src/module_KPT.jl @@ -78,7 +78,11 @@ using Roots: find_zero, Bisection # to find wetness for a given hydraulic conduc ### # from LWFBrook90R: real(kind=8) :: MvGl ! tortuosity parameter (default = 0.5) -"""SOILPAR()\n Function that defines constant parameters from input_siteparam. TO BE REDEFINED +""" + + SOILPAR() + +Define constant parameters from input_siteparam. TO BE REDEFINED """ function SOILPAR(p_RHOWG, p_THICK, # layer thicknesses (mm) @@ -241,12 +245,17 @@ function SOILPAR_MvG(p_RHOWG,p_THICK,p_THETAF,p_THSAT,p_STONEF,p_BEXP, return (p_ψg, p_SWATMX, p_WETfc, p_CHm, p_CHn, p_Ksat, p_PSIF, p_THETAF)#, u_aux_WETNES, u_aux_SWATI) end -""" derive_auxiliary_SOILVAR(u_SWATI, p_SWATMX, p_THSAT, +""" + + derive_auxiliary_SOILVAR(u_SWATI, p_SWATMX, p_THSAT, p_PSIF, p_BEXP, p_WETINF, p_WETF, p_CHM, p_CHN, p_KF, p_θr, p_MvGα, p_MvGn, p_MvGl, p_Ksat, - p_PSIG, NLAYER, IMODEL)\n -Derives alternative representations of soil water status. -I.e. based on the state u_SWATI it returns (u_aux_WETNES, u_aux_PSIM, u_aux_PSITI, p_fu_KK)""" + p_PSIG, NLAYER, IMODEL) + +Derive alternative representations of soil water status. + +Based on the state u_SWATI it returns (u_aux_WETNES, u_aux_PSIM, u_aux_PSITI, p_fu_KK) +""" function derive_auxiliary_SOILVAR(u_SWATI, p_SWATMX, p_THSAT, p_PSIF, p_BEXP, p_WETINF, p_WETF, p_CHM, p_CHN, p_KF, p_θr, p_MvGα, p_MvGn, p_MvGl, p_Ksat, @@ -330,10 +339,13 @@ end # TODO(bernhard): remove remark that I renamed LWFBrook90_MvG_FK() to FK_MvG() -""" FK_MvG(WETNES, KSAT, MvGl, MvGn)\n Computes hydraulic conductivity from -wetness for the Mualem van Genuchten parametrization.\n\n +""" + + FK_MvG(WETNES, KSAT, MvGl, MvGn) + +Compute hydraulic conductivity from wetness for the Mualem van Genuchten parametrization. -Computes unsaturated hydraulic conductivity: K(Se) a.k.a. K(W) using MvG equation 8) +Compute unsaturated hydraulic conductivity: K(Se) a.k.a. K(W) using MvG equation 8) K = Ks*W^l*[ 1 - (1-W^(1/m))^m ]^2 using m = 1-1/n yields: K = Ks*W^l*[ 1 - (1-W^(n/(n-1)))^(1-1/n) ]^2 """ function FK_MvG(WETNES, KSAT, MvGl, MvGn) @@ -344,18 +356,17 @@ function FK_MvG(WETNES, KSAT, MvGl, MvGn) return KSAT*AWET^MvGl*(1 - (1-AWET^(MvGn/(MvGn-1)) )^(1-1/MvGn) )^2 end -""" FPSIMF_CH(u_aux_WETNES,p_PSIF, p_BEXP, p_WET∞, p_WETF, p_CHM, p_CHN)\n Computes ψ(Se) = h(Se) a.k.a ψ(W) = h(W) """ + FPSIMF_CH(u_aux_WETNES,p_PSIF, p_BEXP, p_WET∞, p_WETF, p_CHM, p_CHN) -""" FPSIM_MvG(u_aux_WETNES, p_MvGα, p_MvGn)\n Computes ψ(Se) = h(Se) a.k.a ψ(W) = h(W) +Compute ψ(Se) = h(Se) a.k.a ψ(W) = h(W). +""" + +""" + FPSIM_MvG(u_aux_WETNES, p_MvGα, p_MvGn) + +Compute ψ(Se) = h(Se) a.k.a ψ(W) = h(W). """ -#TODO(bernhard): check correct usage in SOILVAR() i.e. split in FPSIM_MvG() and FPSIMF_CH -# old definition: TODO: remove function LWFBrook90_MvG_FPSIMF(u_aux_WETNES, -# old definition: TODO: remove p_PSIF, p_BEXP, p_WET∞, p_WETF, p_CHM, p_CHN, -# old definition: TODO: remove p_MvGα, p_MvGn, -# old definition: TODO: remove iModel) -# old definition: TODO: replace by (for iModel = 0): FPSIMF_CH(u_aux_WETNES,p_PSIF, p_BEXP, p_WET∞, p_WETF, p_CHM, p_CHN) -# old definition: TODO: replace by (for iModel = 1): FPSIM_MvG(u_aux_WETNES, p_MvGα, p_MvGn) function FPSIMF_CH(u_aux_WETNES,p_PSIF, p_BEXP, p_WET∞, p_WETF, p_CHM, p_CHN) # Computes ψ(Se) = h(Se) a.k.a ψ(W) = h(W) # @@ -441,9 +452,10 @@ function FDPSIDWF_MvG(u_aux_WETNES, α, n) # d PSI/d WETNES, kPa end -# TODO(bernhard): remove remark that I renamed LWFBrook90_MvG_FTheta() to FTheta_MvG() -"""FTheta_MvG(u_aux_WETNES, p_θs, p_θr, iModel) -Computes θ based on Se. +""" + FTheta_MvG(u_aux_WETNES, p_θs, p_θr, iModel) + +Compute θ based on Se. """ function FTheta_MvG(u_aux_WETNES, p_θs, p_θr) # Computes θ(Se) = Se*(θs-θr) + θr @@ -460,11 +472,9 @@ function FTheta_CH(u_aux_WETNES, p_θs) end # TODO(bernhard): remove remark that I renamed LWFBrook90_MvG_FWETNES() to FWETNES_MvG() -""" FWETNES_MvG(u_aux_PSIM, p_MvGα, p_MvGn)\n -Computes θ(ψ) = θ(h) by computing first Se(ψ)=Se(h) a.k.a W(ψ)=W(h) """ + FWETNES_MvG(u_aux_PSIM, p_MvGα, p_MvGn) -""" FWETNES_CH(u_aux_PSIM,p_WETF, p_WET∞, p_BEXP, p_PSIF, p_CHM, p_CHN)\n Computes θ(ψ) = θ(h) by computing first Se(ψ)=Se(h) a.k.a W(ψ)=W(h) """ function FWETNES_MvG(u_aux_PSIM, p_MvGα, p_MvGn) @@ -486,6 +496,11 @@ function FWETNES_MvG(u_aux_PSIM, p_MvGα, p_MvGn) return WETNEs end +""" + FWETNES_CH(u_aux_PSIM,p_WETF, p_WET∞, p_BEXP, p_PSIF, p_CHM, p_CHN) + +Computes θ(ψ) = θ(h) by computing first Se(ψ)=Se(h) a.k.a W(ψ)=W(h) +""" function FWETNES_CH(u_aux_PSIM,p_WETF, p_WET∞, p_BEXP, p_PSIF, p_CHM, p_CHN) # Computes θ(ψ) = θ(h) by computing first Se(ψ)=Se(h) a.k.a W(ψ)=W(h) # diff --git a/src/module_PET.jl b/src/module_PET.jl index 3fe07f5e..0a4ec2d3 100644 --- a/src/module_PET.jl +++ b/src/module_PET.jl @@ -167,9 +167,12 @@ export LWFBrook90_CANOPY, ROUGH, WEATHER, SWPE, SWGE, SWGRA, SRSC, ESAT using ..CONSTANTS # https://discourse.julialang.org/t/large-programs-structuring-modules-include-such-that-to-increase-performance-and-readability/29102/5 """ -LWFBrook90_CANOPY() computes evolution of plant parameters over the season. + LWFBrook90_CANOPY() -Ecoshift: +Compute evolution of plant parameters over the season. + +# Ecoshift +" Subroutine CANOPY calculates plant "parameters" that can vary with day of the year ( DOY). @@ -224,6 +227,7 @@ RTLEN are all reduced proportionally to DENSEF, and RPLANT is increased. However NOT reduce HEIGHT because the remaining canopy still has the same height. Therefore DENSEF should NOT be set to 0 to simulate a clearcut as HEIGHT is unchanged and the aerodynamic resistances will be wrong. Probably DENSEF should not be less than 0.05. +" """ function LWFBrook90_CANOPY(p_fT_HEIGHT, p_fT_LAI, # leaf area index, m2/m2, minimum of 0.00001 @@ -263,9 +267,13 @@ function LWFBrook90_CANOPY(p_fT_HEIGHT, return (p_fu_HEIGHTeff, p_fu_LAIeff, p_fT_SAIeff, p_fu_RTLEN, p_fu_RPLANT) end -"""ROUGH() computes canopy roughness height. +""" + + ROUGH() -Ecoshift: +Compute canopy roughness height. + +# Ecoshift: ROUGH obtains the roughness parameter, z0 , and the zero-plane displacement, d, based on canopy height, h, the projected leaf area index, Lp, and the projected stem area index, Sp. The methods used follow Shuttleworth and Gurney (1990) with some modifications. Shuttleworth @@ -356,9 +364,12 @@ end -"""WEATHER() computes solar radiation, temperature and wind speed. +""" + WEATHER() -Ecoshift: +Compute solar radiation, temperature and wind speed. + +# Ecoshift WEATHER includes all adjustments of input weather data, including separation into daytime and nighttime values. @@ -465,8 +476,11 @@ function WEATHER(p_fT_TMAX, p_fT_TMIN, p_fT_DAYLEN, p_fT_I0HDAY, p_fT_EA, p_fT_U p_fu_UANTM) # average wind speed for nighttime at ZA, m/s end -""" ESAT(p_fu_TA) calculates saturated vp (kPa) and DELTA=dES/dTA (kPa/K) from temperature based on -Murray J Applied Meteorol 6:203 using as input p_fu_TA (air temperature in °C) +""" + ESAT(p_fu_TA) + +Calculate saturated vp (kPa) and DELTA=dES/dTA (kPa/K) from temperature based on +Murray J Applied Meteorol 6:203 using as input p_fu_TA (air temperature in °C). """ function ESAT(p_fu_TA) # @@ -480,10 +494,13 @@ function ESAT(p_fu_TA) end @doc raw""" -WNDADJ(p_fu_ZA, p_fu_DISP, p_fu_Z0, p_FETCH, p_ZW, p_Z0W) returns ratio of wind speed at -reference height (above canopy) to wind speed at weather station + WNDADJ(p_fu_ZA, p_fu_DISP, p_fu_Z0, p_FETCH, p_ZW, p_Z0W) -Ecoshift: This function estimates the wind speed (UA) at reference height ZA above the +Compute ratio of wind speed at reference height (above canopy) to wind speed at weather station. + +# Ecoshift +" +This function estimates the wind speed (UA) at reference height ZA above the canopy from input wind speed at a remote weather station (UW). Assume that the weather station represents a new surface downwind that has a roughness of z0w (Z0W) and a fetch of F (FETCH). Brutsaert (1982) gives the height of the internal boundary layer, zb, as @@ -506,6 +523,7 @@ where zw (ZW) is the height of wind measurement at the weather station (Federer and d (DISP) is the zero-plane displacement of the canopy. This assumes that the weather station is over a smooth surface so its zero plane displacement can be ignored. If the parameter Z0W is set to zero, then no adjustment is made and ua = uw. +" """ function WNDADJ(p_fu_ZA, p_fu_DISP, p_fu_Z0, p_FETCH, p_ZW, p_Z0W) # Brutsaert (1982) equation 7-39 @@ -518,10 +536,12 @@ end """ -SWPE(AA, ASUBS, VPD, RAA, RAC, RAS, RSC, p_fu_RSS, DELTA)\n -computes Shuttleworth and Wallace (1985) transpiration and ground evaporation + SWPE(AA, ASUBS, VPD, RAA, RAC, RAS, RSC, p_fu_RSS, DELTA) -Ecoshift: +Compute Shuttleworth and Wallace (1985) transpiration and ground evaporation. + +# Ecoshift +" Shuttleworth and Wallace (1985) (SW) modified the Penman-Monteith method to account separately for the different water vapor and sensible heat pathways from the soil and from the leaves. Instead of the two resistances of equation (1), rc and ra, SW define five: rsc, @@ -577,6 +597,7 @@ less in reverse order. The outputs Ec (PRATE) and Es (ERATE) from SWPE are in units of mm/d whereas Lvρw E in (1) is output as W m-2 from function PM. The conversion is ETOM * WTOMJ . +" """ function SWPE(AA, ASUBS, VPD, RAA, RAC, RAS, RSC, p_fu_RSS, DELTA) # AA - net radiation at canopy top minus ground flux, W/m2 @@ -611,9 +632,12 @@ function SWPE(AA, ASUBS, VPD, RAA, RAC, RAS, RSC, p_fu_RSS, DELTA) end """ -SWGE() computes ground evaporation rate (mm/d) using Shuttleworth-Wallace with known transpiration + SWGE() -Ecoshift: +Compute ground evaporation rate (mm/d) using Shuttleworth-Wallace with known transpiration. + +# Ecoshift +" The Shuttleworth-Wallace approach incorporates the energy tradeoff between transpiration and soil evaporation. When transpiration is reduced by low availability of soil water or is zero, BROOK90 uses the new value of transpiration, Ec (ARATE), in subroutine SWGE to get a @@ -622,6 +646,7 @@ into (2) and solving for Lv ρw E gives (14) ... then (15) Lv ρw Es = Lv ρw E - Lv ρw Ec +" """ function SWGE(AA, ASUBS, VPD, RAA, RAS, p_fu_RSS, DELTA, ARATE) # AA - net radiation at canopy top minus ground flux, W/m2 @@ -645,9 +670,12 @@ end """ -SWGRA() + SWGRA() -Ecoshift: +Compute Shuttleworth - Wallace - Gurney Aerodynamic Resistances. + +# Ecoshift +" Shuttleworth - Wallace - Gurney Aerodynamic Resistances The three SW aerodynamic resistances, raa, ras, and rac are obtained in subroutine SWGRA by the methods of Shuttleworth and Gurney (1990). The derivations of their equations are not @@ -720,6 +748,7 @@ the leaf diffusion resistance are in series on each side of flat leaves (Jarvis McNaughton 1986). The two resistances should be summed over each side of the leaf before integrating over the canopy, but this is too complicated for practical application (Choudhury and Monteith 1988). +" """ function SWGRA(UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, p_fu_DISP, p_fu_Z0C, p_fu_DISPC, p_fu_Z0GS, p_LWIDTH, p_RHOTP, p_NN, p_fu_LAI, p_fu_SAI) # UA - wind speed at reference height, m/s @@ -806,9 +835,12 @@ function FRSS(p_RSSA, p_RSSB, p_PSIF_TopLayer, u_aux_PSIM_TopLayer, p_PsiCrit_To end """ -SRSC() computes canopy surface resistance +SRSC() -Ecoshift: +Compute canopy surface resistance. + +# Ecoshift +" This routine obtains the canopy surface resistance, rsc, which is the classic canopy resistance in the Penman-Monteith equation, using the Jarvis (1976) expression for the factors that control the individual leaf resistance, r, and its reciprocal the leaf @@ -930,6 +962,7 @@ and (33) is then which corresponds to Shuttleworth and Gurney (1990) and Saugier and Katerji (1991) when Sp = 0. The combination of (32) and (34) provides the value of rsc (RSC) +" """ function SRSC(RAD, p_fu_TA, VPD, p_fu_LAI, p_fu_SAI, p_GLMIN, p_GLMAX, p_R5, p_CVPD, p_RM, p_CR, p_TL, p_T1, p_T2, p_TH) # RAD - solar radiation on canopy, W/m2 @@ -981,7 +1014,12 @@ end """ -Ecoshift: + PM(AA, VPD, DELTA, RA, RC) + +Compute Penman-Monteith latent heat flux density, W/m2. + +# Ecoshift +" The Penman-Monteith equation is Lv ρw E = (Δ (Rn-S) + c_p ρ Da / ra ) / (Δ + γ + γ (rc/ra)) @@ -1002,6 +1040,7 @@ height of water vapor in the plant canopy. The aerodynamic resistance, ra, is a the turbulent transfer capability of the atmosphere between the effective source height and za. The Penman-Monteith equation is derived from the energy balance equation and the mass transfer equations for sensible and latent heat fluxes (e.g. Brutsaert 1982). +" """ function PM(AA, VPD, DELTA, RA, RC) Pm = (RA * DELTA * AA + p_CPRHO * VPD) / ((DELTA + p_GAMMA) * RA + p_GAMMA * RC) diff --git a/src/module_SNO.jl b/src/module_SNO.jl index 7d1ec078..844900b3 100644 --- a/src/module_SNO.jl +++ b/src/module_SNO.jl @@ -73,8 +73,13 @@ using ..PET: ESAT, SWGRA export SNOFRAC, SNOVAP, SNOENRGY -"""SNOFRAC(p_fT_TMAX, p_fT_TMIN, p_RSTEMP) separates RFAL from SFAL -Ecoshift: +""" + SNOFRAC(p_fT_TMAX, p_fT_TMIN, p_RSTEMP) + +Separate RFAL from SFAL. + +# Ecoshift +" Separation of daily precipitation into snow or rain is a major problem in hydrologic modeling. For instance, if the wrong precipitation form is chosen in December, simulated streamflow from that day's precipitation could be shifted from December to April or vice @@ -90,6 +95,7 @@ where RSTEMP is the "base" temperature for the rain-snow transition. When TMAX < SNOFRC = 1; when TMIN > RSTEMP, SNOFRC = 0. The default value of RSTEMP is -0.5°C because that seems to work best at Hubbard Brook. If precipitation is input more than once a day, the same SNOFRC is used for all precipitation intervals. +" """ function SNOFRAC(p_fT_TMAX, p_fT_TMIN, p_RSTEMP) if (p_fT_TMIN >= p_RSTEMP) @@ -102,13 +108,17 @@ function SNOFRAC(p_fT_TMAX, p_fT_TMIN, p_RSTEMP) return p_fT_SNOFRC end -""" SNOVAP(p_fu_TSNOW, p_fu_TA, p_fT_EA, UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, p_fu_DISP, p_fu_Z0C, p_fu_DISPC, p_fu_Z0GS, p_LWIDTH, p_RHOTP, p_NN, p_fu_LAI, p_fu_SAI, p_KSNVP) -Returns potential snow evaporation (mm/d) +""" + SNOVAP(p_fu_TSNOW, p_fu_TA, p_fT_EA, UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, p_fu_DISP, + p_fu_Z0C,p_fu_DISPC, p_fu_Z0GS, p_LWIDTH, p_RHOTP, p_NN, p_fu_LAI, p_fu_SAI, p_KSNVP) + +Compute potential snow evaporation (mm/d). -Ecoshift: -Evaporation rate from the snowpack, and its negative, condensation, are evaluated using the -aerodynamic flux equation +# Ecoshift +" +Evaporation rate from the snowpack, and its negative, condensation, are evaluated +using the aerodynamic flux equation E = (cp ρ / g Ls rw ) (e0 - ea) / (raa + ras) @@ -153,6 +163,7 @@ evaporation under forests. Note that although evaporation and condensation of water are simulated in SNOVAP, the accompanying latent transfer is not simulated. The snow energy balance in subroutine SNOENRGY is (unfortunately) decoupled from the snow evaporation-condensation process. +" """ function SNOVAP(p_fu_TSNOW, p_fu_TA, p_fT_EA, UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, p_fu_DISP, p_fu_Z0C, p_fu_DISPC, p_fu_Z0GS, p_LWIDTH, p_RHOTP, p_NN, p_fu_LAI, p_fu_SAI, p_KSNVP) @@ -173,10 +184,14 @@ function SNOVAP(p_fu_TSNOW, p_fu_TA, p_fT_EA, UA, p_fu_ZA, p_fu_HEIGHT, p_fu_Z0, return p_fu_PSNVP # potential snow evaporation (mm/d) end -""" SNOENRGY(p_fu_TSNOW, p_fu_TA, p_fT_DAYLEN, p_CCFAC, p_MELFAC, p_fT_SLFDAY, p_fu_LAI, -p_fu_SAI, p_LAIMLT, p_SAIMLT)\nreturns energy flux density to snow surface, MJ m-2 d-1 +""" + SNOENRGY(p_fu_TSNOW, p_fu_TA, p_fT_DAYLEN, p_CCFAC, p_MELFAC, p_fT_SLFDAY, p_fu_LAI, + p_fu_SAI, p_LAIMLT, p_SAIMLT) + +Compute energy flux density to snow surface, MJ m-2 d-1. -Ecoshift: +# Ecoshift +" The energy flux density to the snow surface (SNOEN, MJ m-2 d-1) is calculated in subroutine SNOENRGY independently of precipitation for the day. @@ -205,6 +220,7 @@ below this. Inclusion of SLFDAY in the MELT equation arises from an assumption t radiation melt plays an important role. If this is not so, then SLFDAY multiplier could be omitted, and snowmelt would not depend on slope-aspect. The functional forms of the SAI and LAI dependencies are based on the somewhat arbitrary curves used by Federer and Lash (1978). +" """ function SNOENRGY(p_fu_TSNOW, p_fu_TA, p_fT_DAYLEN, p_CCFAC, p_MELFAC, p_fT_SLFDAY, p_fu_LAI, p_fu_SAI, p_LAIMLT, p_SAIMLT) # snow surface energy balance @@ -220,11 +236,13 @@ function SNOENRGY(p_fu_TSNOW, p_fu_TA, p_fT_DAYLEN, p_CCFAC, p_MELFAC, p_fT_SLFD end """ -SNOWPACK(p_fu_RTHR, p_fu_STHR, p_fu_PSNVP, p_fu_SNOEN, u_CC, u_SNOW, - u_SNOWLQ, p_DTP, p_fu_TA, p_MAXLQF, p_GRDMLT)\n -Updates status of snowpack. + SNOWPACK(p_fu_RTHR, p_fu_STHR, p_fu_PSNVP, p_fu_SNOEN, u_CC, u_SNOW, + u_SNOWLQ, p_DTP, p_fu_TA, p_MAXLQF, p_GRDMLT) + +Update status of snowpack. -Ecoshift: +# Ecoshift +" In each precipitation interval, throughfall of rain (RTHR) and throughfall of snow (STHR) are calculated and subroutine SNOWPACK is entered if there is STHR or if there is SNOW. This subroutine adds throughfall to the snowpack, subtracts groundmelt, snow evaporation, and @@ -283,6 +301,7 @@ MSBPREINT then calculates the rain passing through the snow (RNET) as RNET = RTHR - RSNO. When SNOW exists at the beginning of the day, soil evaporation (SLVP) is zero. +" """ function SNOWPACK(p_fu_RTHR, p_fu_STHR, p_fu_PSNVP, p_fu_SNOEN, u_CC, u_SNOW, u_SNOWLQ, p_DTP, p_fu_TA, p_MAXLQF, p_GRDMLT) diff --git a/src/module_SUN.jl b/src/module_SUN.jl index 108176ff..a6f1d4eb 100644 --- a/src/module_SUN.jl +++ b/src/module_SUN.jl @@ -49,6 +49,12 @@ using ..CONSTANTS: p_SIGMA, p_PI # https://discourse.julialang.org/t/large-progr export EQUIVSLP, SUNDS, AVAILEN """ + EQUIVSLP(p_LAT, SLOPE, p_ASPECT) + +Correct solar radiation for slope and aspect and compute "equivalent slope" parameters. + +# Ecoshift +" Correction of solar radiation for slope and aspect requires calculation of an "equivalent slope" in subroutine EQUIVSLP. The equivalent slope is defined as the location on the earth's surface where a horizontal surface is parallel to the given sloping surface. @@ -64,6 +70,7 @@ L2 = ATAN { SIN(SLOPE) * SIN(ASPECT) /[ COS(SLOPE) * COS(LAT) - SIN(SLOPE) * SIN with fixes for negative or zero denominator in L2, where SLOPE (ESLOPE), ASPECT, and latitude (LAT) are input parameters describing the location. All angles in the subroutine are in radians. +" """ function EQUIVSLP(p_LAT, SLOPE, p_ASPECT) #Swift#s p_L1 and p_L2, Lee (3.31, 3.32) @@ -79,9 +86,13 @@ function EQUIVSLP(p_LAT, SLOPE, p_ASPECT) return (p_L1, p_L2) end -""" Function SUNDS() returns p_fT_DAYLEN, p_fT_I0HDAY, p_fT_SLFDAY +""" + SUNDS() + +Return p_fT_DAYLEN, p_fT_I0HDAY, p_fT_SLFDAY. -From ecoshift: +# Ecoshift +" Several radiation-related variables depend only on day of the year and location. These are calculated in SUNDS, which is called once a day. @@ -126,6 +137,7 @@ from the I0SDAY equation with L1 = LAT, L2 = 0, and T3 and T2 for a horizontal s LAT. The daylength (DAYLEN), which is the fraction of a day that the sun is above a horizontal horizon, is HAFDAY / π where function HAFDAY is used with L = LAT. SLFDAY is the ratio of I0SDAY to I0HDAY. SUNDS outputs DAYLEN, I0HDAY, and SLFDAY. +" """ function SUNDS(p_LAT, SLOPE, DOY, p_L1, p_L2, p_SC, p_PI, p_WTOMJ) SCD = p_SC / (1.0 - 0.0167 * cos(0.0172 * (DOY - 3))) ^ 2 @@ -198,10 +210,14 @@ end -""" AVAILEN(SLRAD, p_fu_ALBEDO, p_C1, p_C2, p_C3, p_fu_TA, p_fT_EA, RATIO, p_fu_SHEAT, p_CR, p_fu_LAI, p_fu_SAI)\n -estimates the available energy above and below the canopy. +""" + AVAILEN(SLRAD, p_fu_ALBEDO, p_C1, p_C2, p_C3, p_fu_TA, p_fT_EA, RATIO, p_fu_SHEAT, p_CR, + p_fu_LAI, p_fu_SAI) + +Estimate the available energy above and below the canopy. -Ecoshift: +# Ecoshift +" Estimates of available energy above and below the canopy are made in subroutine AVAILEN. Available energy is net radiation minus subsurface heat flux (SHEAT), and is the energy available for partitioning into heating the air and evaporating water. SHEAT is set to zero @@ -305,6 +321,7 @@ plays a major role in the value of PE estimates when net radiation is estimated. effort using worldwide longwave data (not estimates from models!) will be needed to improve this situation. Fortunately, the cloud cover correction, approximate as it is, brings the net longwave closer to zero and helps wash out the emissivity error. +" """ function AVAILEN(SLRAD, p_fu_ALBEDO, p_C1, p_C2, p_C3, p_fu_TA, p_fT_EA, RATIO, p_fu_SHEAT, p_CR, p_fu_LAI, p_fu_SAI) diff --git a/src/module_WAT.jl b/src/module_WAT.jl index e3bb5694..68a559e2 100644 --- a/src/module_WAT.jl +++ b/src/module_WAT.jl @@ -47,13 +47,16 @@ module WAT # WATER MOVEMENT IN SOIL export INFPAR, LWFRootGrowth, ITER -"""INFPAR(p_INFEXP, p_ILAYER, p_THICK, NLAYER)\n Computes fraction of infiltration to each -soil layer.\n\n -Arguments: -- p_INFEXP: infiltration exponent: 0 all to top, 1 uniform with depth, >1.0=more at bottom than at top -- p_ILAYER: number of layers over which infiltration is distributed -- p_THICK -- NLAYER +""" + INFPAR(p_INFEXP, p_ILAYER, p_THICK, NLAYER) + +Compute fraction of infiltration to each soil layer. + +# Arguments: +- `p_INFEXP`: infiltration exponent: 0 all to top, 1 uniform with depth, >1.0=more at bottom than at top +- `p_ILAYER`: number of layers over which infiltration is distributed +- `p_THICK` +- `NLAYER` """ function INFPAR(p_INFEXP, p_ILAYER, p_THICK, NLAYER) p_INFRAC = fill(NaN, NLAYER) # fraction of infiltration to each layer @@ -80,18 +83,21 @@ end -"""LWFRootGrowth(frelden, tini, age, rgroper, inirdep, inirlen, NLAYER)\n Computes root growth -according to LWF root growth model, (Hammel and Kennel 2000).\n\n -Arguments: - - frelden[] : final relative values of root length per unit volume - - tini[] : initial time for root growth in layer - - age : age of vegetation - - rgroper : period of root growth in layer, a - - inirdep : intial root depth, m - - inirlen : intial total root length, m m-2 - - NLAYER : number of soil layers +""" + LWFRootGrowth(frelden, tini, age, rgroper, inirdep, inirlen, NLAYER) + +Compute root growth according to LWF root growth model, (Hammel and Kennel 2000). + +# Arguments: + - `frelden[]`: final relative values of root length per unit volume + - `tini[]` : initial time for root growth in layer + - `age` : age of vegetation + - `rgroper` : period of root growth in layer, a + - `inirdep` : intial root depth, m + - `inirlen` : intial total root length, m m-2 + - `NLAYER` : number of soil layers Returns: - - RELDEN[] : current, age-dependent relative values of root length per unit volume + - `RELDEN[]` : current, age-dependent relative values of root length per unit volume """ function LWFRootGrowth(frelden, tini, age, rgroper, inirdep, inirlen, NLAYER) @@ -121,10 +127,12 @@ end """ -BYFLFR(NLAYER, p_BYPAR, p_QFPAR, p_QFFC, u_aux_WETNES, p_WETF)\n -TODO(bernhard): + BYFLFR(NLAYER, p_BYPAR, p_QFPAR, p_QFFC, u_aux_WETNES, p_WETF) + +Compute fraction of bypass flow. -Ecoshift: +# Ecoshift +" Bypass flow (BYFL) and surface or source area flow (SRFL) are the two stormflow or quickflow generating mechanisms in BROOK90. The conceptual difference is that SRFL is "new" water that has not infiltrated but has moved across the surface to a channel, whereas BYFL is "new" @@ -168,6 +176,7 @@ produces a constant BYFRAC of QFFC. Note that BYFRAC is calculated from soil water prior to the input of water for the time step. +" """ function BYFLFR(NLAYER, p_BYPAR, p_QFPAR, p_QFFC, u_aux_WETNES, p_WETF) # TODO(bernhard): could be optimized by not allocating each time new memory (versus in-place) @@ -211,7 +220,10 @@ function BYFLFR(NLAYER, p_BYPAR, p_QFPAR, p_QFFC, u_aux_WETNES, p_WETF) end -""" DSLOP() downslope flow rate from layer +""" + DSLOP() + +Compute downslope flow rate from layer. """ function DSLOP(p_DSLOPE, p_LENGTH, p_THICK_i, p_STONEF_i, u_aux_PSIM_i, p_RHOWG, p_fu_KK_i) @@ -228,7 +240,10 @@ function DSLOP(p_DSLOPE, p_LENGTH, p_THICK_i, p_STONEF_i, u_aux_PSIM_i, p_RHOWG, return aux_du_DSFLI end -""" VERT() vertical flow rate +""" + VERT() + +Compute vertical flow rate. """ function VERT(KK_i, KK_iplus1, KSAT_i, KSAT_iplus1, @@ -254,9 +269,13 @@ function VERT(KK_i, KK_iplus1, VRFLI = (GRAD * KK_mean / p_RHOWG) * (1 - (STONEF_i + STONEF_iplus1) / 2) return(VRFLI) end -"""KKMEAN(KK_i, KK_iplus1, THICK_i, THICK_iplus1)\n -computes average hydraulic conductivity. Note that between version 3.1 (where LWFBrook90 was -forked), and 4.8 of Brook90 there were different variants how to compute the average. +""" + KKMEAN(KK_i, KK_iplus1, THICK_i, THICK_iplus1) + +Compute average hydraulic conductivity. + +Note that between version 3.1 (where LWFBrook90 was forked), and 4.8 of Brook90 there were +different variants how to compute the average. """ function KKMEAN(KK_i, KK_iplus1, THICK_i, THICK_iplus1) # NOTE(bernhard): different averaging for KKMEAN and GRAD exist in different Brook90 versions: @@ -277,9 +296,14 @@ end -"""net inflow to soil layer +""" + INFLOW(NLAYER, DTI, p_INFRAC, p_fu_BYFRAC, p_fu_SLFL, + aux_du_DSFLI, aux_du_TRANI, aux_du_SLVP, p_SWATMX, u_SWATI, VRFLI_prior) + +Compute net inflow to soil layer. -Ecoshift: +# Ecoshift +" In this routine, infiltrating water (SLFL) is allocated to soil water in each layer (INFLIi ) and to bypass flow from each layer (BYFLIi ). The fraction of SLFL going to each layer (INFRACi ) is constant and is obtained in subroutine INFPAR. This fraction is separated into @@ -326,6 +350,7 @@ VRFLI(0) to reduce, excess water becomes negative INFLI(1) and increases BYFLI(1 The modified values of VRFLIi are output from the INFLOW routine as variable VV because the original VRFLIi are needed again if the iteration time step (DTI) is reduced. +" """ function INFLOW(NLAYER, DTI, p_INFRAC, p_fu_BYFRAC, p_fu_SLFL, aux_du_DSFLI, aux_du_TRANI, aux_du_SLVP, p_SWATMX, u_SWATI, VRFLI_prior) @@ -602,7 +627,10 @@ function ITER(IMODEL, NLAYER, DTI, DTIMIN, DPSIDW, return DTINEW # second estimate of DTI end -"""calculates groundwater flow and seepage loss +""" + GWATER(u_GWAT, p_GSC, p_GSP, p_DT, aux_du_VRFLIN) + +Calculate groundwater flow and seepage loss. """ function GWATER(u_GWAT, p_GSC, p_GSP, p_DT, aux_du_VRFLIN) if (p_GSC < 1.0e-8)