diff --git a/configs/README.md b/configs/README.md index 4687c7e..60ea2c4 100644 --- a/configs/README.md +++ b/configs/README.md @@ -16,10 +16,11 @@ A detailed description of the parameters for model configuration (i.e., initiali | endtime | double (scalar)| >0 | sec, min, hr, d | simulation duration | - | time at which model simulation ends | | layer_soil_type | int (1D array) | - | - | state variable | - | layer soil type (read from the database file soil_params_file) | | max_soil_types | int | >1 | - | - | - | maximum number of soil types read from the file soil_params_file (default is set to 15) | -| wilting_point_psi | double (scalar) | - | cm | state variable | - | wilting point (the amount of water not available for plants) used in computing AET | +| wilting_point_psi | double (scalar) | - | cm | state variable | - | wilting point (the amount of water not available for plants) used in computing AET. Suggested value is 15495.0 cm, corresponding to 15 atm. | +| field_capacity_psi | double (scalar) | - | cm | state variable | - | capillary head corresponding to volumetric water content at which gravity drainage becomes slower, used in computing AET. Suggested value is 340.9 cm for most soils, corresponding to 1/3 atm, and 103.3 cm for sands, corresponding to 1/10 atm. | | use_closed_form_G | bool | true or false | - | - | - | determines whether the numeric integral or closed form for G is used; a value of true will use the closed form. This defaults to false. | | giuh_ordinates | double (1D array)| - | - | state parameter | - | GIUH ordinates (for giuh based surface runoff) | | verbosity | string | high, low, none | - | debugging | - | controls IO (screen outputs and writing to disk) | | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen| | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) | -| calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_m`, and `vg_alpha` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content | +| calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content | diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt index 82297eb..06b6e89 100644 --- a/configs/config_lasam_Bushland.txt +++ b/configs/config_lasam_Bushland.txt @@ -11,4 +11,5 @@ use_closed_form_G=false layer_soil_type=16,17,18 max_soil_types=25 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.06,0.51,0.28,0.12,0.03 diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt index 17b4e2f..1d230cf 100644 --- a/configs/config_lasam_Phillipsburg.txt +++ b/configs/config_lasam_Phillipsburg.txt @@ -11,4 +11,5 @@ use_closed_form_G=false layer_soil_type=13,14,15 max_soil_types=25 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.06,0.51,0.28,0.12,0.03 diff --git a/configs/config_lasam_sft_ngen.txt b/configs/config_lasam_sft_ngen.txt index f3992bc..e9ff131 100644 --- a/configs/config_lasam_sft_ngen.txt +++ b/configs/config_lasam_sft_ngen.txt @@ -11,6 +11,7 @@ use_closed_form_G=false layer_soil_type=13,14,15 max_soil_types=25 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.06,0.51,0.28,0.12,0.03 sft_coupled=true soil_z=10,20,30,40,50,60,70,80,90,100.0,110.,120,130.,140.,150.,160.,170.,180.,190.,200.0[cm] diff --git a/include/all.hxx b/include/all.hxx index 67b7ebc..6b7526c 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -3,7 +3,7 @@ #define _ALL_HXX /* - authors : Ahmad Jan and Fred Ogden + authors : Ahmad Jan and Fred Ogden and Peter La Follette year : 2022 email : ahmad.jan@noaa.gov - This header file constains functions' definitions used in the lgar.cxx, bmi_lasam.cxx and in other files. @@ -128,6 +128,7 @@ struct lgar_bmi_parameters depth from the surface in meters */ double *frozen_factor; // frozen factor added to the hydraulic conductivity due to coupling to soil freeze-thaw double wilting_point_psi_cm; // wilting point (the amount of water not available for plants or not accessible by plants) + double field_capacity_psi_cm; // field capacity represented as a capillary head. Note that both wilting point and field capacity are specified for the whole model domain with single values bool use_closed_form_G = false; /* true if closed form of capillary drive calculation is desired, false if numeric integral for capillary drive calculation is desired */ double ponded_depth_cm; // amount of water on the surface unavailable for surface runoff @@ -186,11 +187,14 @@ struct lgar_mass_balance_variables // the structure holds pointer bmi output variables struct lgar_calib_parameters { - double *theta_e; // theta_e = smcmax - double *theta_r; // theta_r = smcmin - double *vg_alpha; // Van Genuchton alpha - double *vg_m; // Van Genuchton m - double *Ksat; // Hydraulic conductivity [cm/hr] + double *theta_e; // theta_e = smcmax [-] + double *theta_r; // theta_r = smcmin [-] + double *vg_n; // Van Genuchten n [-] + double *vg_alpha; // Van Genuchten alpha [1/cm] + double *Ksat; // Hydraulic conductivity [cm/hr] + double field_capacity_psi; // field capacity in capillary head [cm] + double ponded_depth_max; // maximum ponded depth of surface water [cm] + }; // nested structure of structures; main structure for the use in bmi @@ -332,7 +336,7 @@ extern void InitializeWettingFronts(int num_layers, double initial_psi_cm, int * /*Other function prototypes for doing hydrology calculations, etc. */ /********************************************************************/ -extern double calc_aet(double PET_timestep_cm, double timestep_h, double wilting_point_psi_cm, int *soil_type, +extern double calc_aet(double PET_timestep_cm, double timestep_h, double wilting_point_psi_cm, double field_capacity_psi_cm, int *soil_type, double AET_thresh_Theta, double AET_expon, struct wetting_front* head, struct soil_properties_ *soil_props); /********************************************************************/ diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index 9f10a4f..a39fb7f 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -67,9 +67,11 @@ public: // calibratable parameters this->calib_var_names[0] = "smcmax"; this->calib_var_names[1] = "smcmin"; - this->calib_var_names[2] = "van_genuchten_m"; + this->calib_var_names[2] = "van_genuchten_n"; this->calib_var_names[3] = "van_genuchten_alpha"; this->calib_var_names[4] = "hydraulic_conductivity"; + this->calib_var_names[5] = "field_capacity"; + this->calib_var_names[6] = "ponded_depth_max"; }; void Initialize(std::string config_file); @@ -132,7 +134,7 @@ private: struct model_state* state; static const int input_var_name_count = 3; static const int output_var_name_count = 15; - static const int calib_var_name_count = 5; + static const int calib_var_name_count = 7; std::string input_var_names[input_var_name_count]; std::string output_var_names[output_var_name_count]; diff --git a/realizations/realization_config_lasam.json b/realizations/realization_config_lasam.json index 26cebb6..7066b47 100644 --- a/realizations/realization_config_lasam.json +++ b/realizations/realization_config_lasam.json @@ -53,6 +53,15 @@ "precipitation_rate" : "P", "potential_evapotranspiration_rate" : "PET" }, + "model_params": { + "smcmax": [0.50, 0.49, 0.48], + "smcmin": [0.09, 0.08, 0.07], + "van_genuchten_n": [2.0, 2.1, 2.2], + "van_genuchten_alpha": [0.009, 0.008, 0.007], + "hydraulic_conductivity": [10.0, 11.0, 12.0], + "field_capacity": 333.3, + "ponded_depth_max": 1.1 + }, "output_variables" : [ "precipitation", "potential_evapotranspiration", diff --git a/realizations/realization_config_lasam_sft.json b/realizations/realization_config_lasam_sft.json index 13d5f40..2248399 100644 --- a/realizations/realization_config_lasam_sft.json +++ b/realizations/realization_config_lasam_sft.json @@ -88,6 +88,15 @@ "precipitation_rate" : "P", "potential_evapotranspiration_rate" : "PET" }, + "model_params": { + "smcmax": [0.50, 0.49, 0.48], + "smcmin": [0.09, 0.08, 0.07], + "van_genuchten_n": [2.0, 2.1, 2.2], + "van_genuchten_alpha": [0.009, 0.008, 0.007], + "hydraulic_conductivity": [10.0, 11.0, 12.0], + "field_capacity": 333.3, + "ponded_depth_max": 1.1 + }, "output_variables" : [ "precipitation", "potential_evapotranspiration", diff --git a/realizations/realization_config_lasam_smp.json b/realizations/realization_config_lasam_smp.json index c4c1e6f..3203240 100644 --- a/realizations/realization_config_lasam_smp.json +++ b/realizations/realization_config_lasam_smp.json @@ -74,6 +74,15 @@ "precipitation_rate" : "P", "potential_evapotranspiration_rate" : "PET" }, + "model_params": { + "smcmax": [0.50, 0.49, 0.48], + "smcmin": [0.09, 0.08, 0.07], + "van_genuchten_n": [2.0, 2.1, 2.2], + "van_genuchten_alpha": [0.009, 0.008, 0.007], + "hydraulic_conductivity": [10.0, 11.0, 12.0], + "field_capacity": 333.3, + "ponded_depth_max": 1.1 + }, "output_variables" : [ "precipitation", "potential_evapotranspiration", diff --git a/src/aet.cxx b/src/aet.cxx index 60cb8ab..03e4bfa 100644 --- a/src/aet.cxx +++ b/src/aet.cxx @@ -4,7 +4,7 @@ #include "../include/all.hxx" //################################################################################ -/* authors : Fred Ogden and Ahmad Jan +/* authors : Fred Ogden and Ahmad Jan and Peter La Follette year : 2022 the code computes actual evapotranspiration given PET. It uses an S-shaped function used in HYDRUS-1D (Simunek & Sejna, 2018). @@ -14,7 +14,7 @@ //################################################################################ -extern double calc_aet(double PET_timestep_cm, double time_step_h, double wilting_point_psi_cm, +extern double calc_aet(double PET_timestep_cm, double time_step_h, double wilting_point_psi_cm, double field_capacity_psi_cm, int *soil_type, double AET_thresh_Theta, double AET_expon, struct wetting_front* head, struct soil_properties_ *soil_properties) { @@ -44,8 +44,8 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin vg_n = soil_properties[soil_num].vg_n; // compute theta field capacity - double head_at_which_PET_equals_AET_cm = 340.9;//*10/33; //340.9 is 0.33 atm, expressed in water depth, which is a good field capacity for most soils. - //Coarser soils like sand will have a field capacity of 0.1 atm or so. + double head_at_which_PET_equals_AET_cm = field_capacity_psi_cm; //340.9 is 0.33 atm, expressed in water depth, which is a good field capacity for most soils. + //Coarser soils like sand will have a field capacity of 0.1 atm or so, which would be 103.3 cm. double theta_fc = calc_theta_from_h(head_at_which_PET_equals_AET_cm, vg_a,vg_m, vg_n, theta_e, theta_r); double wp_head_theta = calc_theta_from_h(wilting_point_psi_cm, vg_a,vg_m, vg_n, theta_e, theta_r); diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 7d8002b..c11a8aa 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -115,6 +115,7 @@ Update() double subtimestep_h = state->lgar_bmi_params.timestep_h; int nint = state->lgar_bmi_params.nint; double wilting_point_psi_cm = state->lgar_bmi_params.wilting_point_psi_cm; + double field_capacity_psi_cm = state->lgar_bmi_params.field_capacity_psi_cm; bool use_closed_form_G = state->lgar_bmi_params.use_closed_form_G; // constant value used in the AET function @@ -192,7 +193,7 @@ Update() // Calculate AET from PET if PET is non-zero if (PET_subtimestep_cm_per_h > 0.0) { - AET_subtimestep_cm = calc_aet(PET_subtimestep_cm_per_h, subtimestep_h, wilting_point_psi_cm, + AET_subtimestep_cm = calc_aet(PET_subtimestep_cm_per_h, subtimestep_h, wilting_point_psi_cm, field_capacity_psi_cm, state->lgar_bmi_params.layer_soil_type, AET_thresh_Theta, AET_expon, state->head, state->soil_properties); } @@ -509,27 +510,27 @@ update_calibratable_parameters() double volstart_before = lgar_calc_mass_bal(state->lgar_bmi_params.cum_layer_thickness_cm, state->head); - for (int i=0; ilgar_bmi_params.num_wetting_fronts; i++) { + for (int i=0; ilgar_bmi_params.num_wetting_fronts; i++) {//first we update the parameters that depend on soil layer, for each layer layer_num = current->layer_num; soil = state->lgar_bmi_params.layer_soil_type[layer_num]; assert (current != NULL); if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters (initial values) ----------- \n"; + std::cerr<<"----------- Calibratable parameters depending on soil layer (initial values) ----------- \n"; std::cerr<<"| soil_type = "<< soil <<", layer = "<soil_properties[soil].theta_e = state->lgar_calib_params.theta_e[layer_num-1]; state->soil_properties[soil].theta_r = state->lgar_calib_params.theta_r[layer_num-1]; - state->soil_properties[soil].vg_m = state->lgar_calib_params.vg_m[layer_num-1]; - state->soil_properties[soil].vg_n = 1.0/(1.0 - state->soil_properties[soil].vg_m); + state->soil_properties[soil].vg_n = state->lgar_calib_params.vg_n[layer_num-1]; + state->soil_properties[soil].vg_m = 1.0 - 1.0/state->soil_properties[soil].vg_n; state->soil_properties[soil].vg_alpha_per_cm = state->lgar_calib_params.vg_alpha[layer_num-1]; state->soil_properties[soil].Ksat_cm_per_h = state->lgar_calib_params.Ksat[layer_num-1]; @@ -538,18 +539,34 @@ update_calibratable_parameters() state->soil_properties[soil].theta_e, state->soil_properties[soil].theta_r); if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters (updated values) ----------- \n"; + std::cerr<<"----------- Calibratable parameters depending on soil layer (updated values) ----------- \n"; std::cerr<<"| soil_type = "<< soil <<", layer = "<next; } + + //next we update the parameters that apply to the whole model domain and do not depend on soil layer + if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { + std::cerr<<"----------- Calibratable parameters independent of soil layer (initial values) ----------- \n"; + std::cerr<<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm + <<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm <<"\n"; + } + + state->lgar_bmi_params.field_capacity_psi_cm = state->lgar_calib_params.field_capacity_psi; + state->lgar_bmi_params.ponded_depth_max_cm = state->lgar_calib_params.ponded_depth_max; + + if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { + std::cerr<<"----------- Calibratable parameters independent of soil layer (updated values) ----------- \n"; + std::cerr<<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm + <<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm <<"\n"; + } if (verbosity.compare("high") == 0) listPrint(state->head); @@ -581,7 +598,7 @@ GetVarGrid(std::string name) || name.compare("actual_evapotranspiration") == 0) // double return 1; else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 - || name.compare("soil_storage") == 0) // double + || name.compare("soil_storage") == 0 || name.compare("field_capacity") == 0 || name.compare("ponded_depth_max") == 0)// double return 1; else if (name.compare("total_discharge") == 0 || name.compare("infiltration") == 0 || name.compare("percolation") == 0 || name.compare("groundwater_to_stream_recharge") == 0) // double @@ -589,7 +606,7 @@ GetVarGrid(std::string name) else if (name.compare("mass_balance") == 0) return 1; else if (name.compare("soil_depth_layers") == 0 || name.compare("smcmax") == 0 || name.compare("smcmin") == 0 - || name.compare("van_genuchten_m") == 0 || name.compare("van_genuchten_alpha") == 0 + || name.compare("van_genuchten_m") == 0 || name.compare("van_genuchten_alpha") == 0 || name.compare("van_genuchten_n") == 0 || name.compare("hydraulic_conductivity") == 0) // array of doubles (fixed length) return 2; else if (name.compare("soil_moisture_wetting_fronts") == 0 || name.compare("soil_depth_wetting_fronts") == 0) // array of doubles (dynamic length) @@ -806,12 +823,16 @@ GetValuePtr (std::string name) return (void*)this->state->lgar_calib_params.theta_e; else if (name.compare("smcmin") == 0) return (void*)this->state->lgar_calib_params.theta_r; - else if (name.compare("van_genuchten_m") == 0) - return (void*)this->state->lgar_calib_params.vg_m; + else if (name.compare("van_genuchten_n") == 0) + return (void*)this->state->lgar_calib_params.vg_n; else if (name.compare("van_genuchten_alpha") == 0) return (void*)this->state->lgar_calib_params.vg_alpha; else if (name.compare("hydraulic_conductivity") == 0) return (void*)this->state->lgar_calib_params.Ksat; + else if (name.compare("ponded_depth_max") == 0) + return (void*)&this->state->lgar_calib_params.ponded_depth_max; + else if (name.compare("field_capacity") == 0) + return (void*)&this->state->lgar_calib_params.field_capacity_psi; else { std::stringstream errMsg; errMsg << "variable "<< name << " does not exist"; diff --git a/src/lgar.cxx b/src/lgar.cxx index bfcf3b9..505b010 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -86,12 +86,15 @@ extern void lgar_initialize(string config_file, struct model_state *state) // initialize array for holding calibratable parameters state->lgar_calib_params.theta_e = new double[state->lgar_bmi_params.num_layers]; state->lgar_calib_params.theta_r = new double[state->lgar_bmi_params.num_layers]; - state->lgar_calib_params.vg_m = new double[state->lgar_bmi_params.num_layers]; + state->lgar_calib_params.vg_n = new double[state->lgar_bmi_params.num_layers]; state->lgar_calib_params.vg_alpha = new double[state->lgar_bmi_params.num_layers]; state->lgar_calib_params.Ksat = new double[state->lgar_bmi_params.num_layers]; // initialize thickness/depth and soil moisture of wetting fronts (used for model coupling) // also initialize calibratable parameters + state->lgar_calib_params.field_capacity_psi = state->lgar_bmi_params.field_capacity_psi_cm; + state->lgar_calib_params.ponded_depth_max = state->lgar_bmi_params.ponded_depth_max_cm; + struct wetting_front *current = state->head; for (int i=0; ilgar_bmi_params.num_wetting_fronts; i++) { assert (current != NULL); @@ -103,7 +106,7 @@ extern void lgar_initialize(string config_file, struct model_state *state) state->lgar_calib_params.theta_e[i] = state->soil_properties[soil].theta_e; state->lgar_calib_params.theta_r[i] = state->soil_properties[soil].theta_r; - state->lgar_calib_params.vg_m[i] = state->soil_properties[soil].vg_m; + state->lgar_calib_params.vg_n[i] = state->soil_properties[soil].vg_n; state->lgar_calib_params.vg_alpha[i] = state->soil_properties[soil].vg_alpha_per_cm; state->lgar_calib_params.Ksat[i] = state->soil_properties[soil].Ksat_cm_per_h; @@ -159,6 +162,7 @@ extern void lgar_initialize(string config_file, struct model_state *state) @param frozen_factor : frozen factor causing the hydraulic conductivity to decrease due to frozen soil (when coupled to soil freeze thaw model) @param wilting_point_psi_cm : wilting point (the amount of water not available for plants or not accessible by plants) + @param field_capacity_psi_cm : field capacity, represented with a capillary head (head above which drainage is much faster) @param ponded_depth_cm : amount of water on the surface not available for surface drainage (initialized to zero) @param ponded_depth_max cm : maximum amount of water on the surface not available for surface drainage (default is zero) @param nint : number of trapezoids used in integrating the Geff function (set to 120) @@ -210,18 +214,19 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.sft_coupled = false; state->lgar_bmi_params.use_closed_form_G = false; - bool is_layer_thickness_set = false; - bool is_initial_psi_set = false; - bool is_timestep_set = false; - bool is_endtime_set = false; - bool is_forcing_resolution_set = false; - bool is_layer_soil_type_set = false; - bool is_wilting_point_psi_cm_set = false; - bool is_soil_params_file_set = false; - bool is_max_soil_types_set = false; - bool is_giuh_ordinates_set = false; - bool is_soil_z_set = false; - bool is_ponded_depth_max_cm_set = false; + bool is_layer_thickness_set = false; + bool is_initial_psi_set = false; + bool is_timestep_set = false; + bool is_endtime_set = false; + bool is_forcing_resolution_set = false; + bool is_layer_soil_type_set = false; + bool is_wilting_point_psi_cm_set = false; + bool is_field_capacity_psi_cm_set = false; + bool is_soil_params_file_set = false; + bool is_max_soil_types_set = false; + bool is_giuh_ordinates_set = false; + bool is_soil_z_set = false; + bool is_ponded_depth_max_cm_set = false; string soil_params_file; @@ -368,6 +373,17 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } + else if (param_key == "field_capacity_psi") { + state->lgar_bmi_params.field_capacity_psi_cm = stod(param_value); + is_field_capacity_psi_cm_set = true; + + if (verbosity.compare("high") == 0) { + std::cerr<<"Field capacity Psi [cm] : "<lgar_bmi_params.field_capacity_psi_cm<<"\n"; + std::cerr<<" ***** \n"; + } + + continue; + } else if (param_key == "use_closed_form_G") { if (param_value == "false") { state->lgar_bmi_params.use_closed_form_G = false; @@ -522,6 +538,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->soil_properties = new soil_properties_[state->lgar_bmi_params.num_soil_types+1]; int num_soil_types = state->lgar_bmi_params.num_soil_types; double wilting_point_psi_cm = state->lgar_bmi_params.wilting_point_psi_cm; + double field_capacity_psi_cm = state->lgar_bmi_params.field_capacity_psi_cm; int max_num_soil_in_file = lgar_read_vG_param_file(soil_params_file.c_str(), num_soil_types, wilting_point_psi_cm, state->soil_properties); @@ -572,7 +589,13 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) if(!is_wilting_point_psi_cm_set) { stringstream errMsg; - errMsg << "The configuration file \'" << config_file <<"\' does not set wilting_point_psi. \n"; + errMsg << "The configuration file \'" << config_file <<"\' does not set wilting_point_psi. \n Recommended value of 15495.0[cm], corresponding to 15 atm. \n"; + throw runtime_error(errMsg.str()); + } + + if(!is_field_capacity_psi_cm_set) { + stringstream errMsg; + errMsg << "The configuration file \'" << config_file <<"\' does not set field_capacity_psi. \n Recommended value of 340.9[cm] for most soils, corresponding to 1/3 atm, or 103.3[cm] for sands, corresponding to 1/10 atm. \n"; throw runtime_error(errMsg.str()); } diff --git a/tests/configs/config_lasam_synth_0.txt b/tests/configs/config_lasam_synth_0.txt index 7b900f5..c6c27ad 100644 --- a/tests/configs/config_lasam_synth_0.txt +++ b/tests/configs/config_lasam_synth_0.txt @@ -9,5 +9,6 @@ forcing_resolution=3600[sec] layer_soil_type=13,14,15 max_soil_types=15 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.06,0.51,0.28,0.12,0.03 ponded_depth_max=1.0[cm] diff --git a/tests/configs/config_lasam_synth_1.txt b/tests/configs/config_lasam_synth_1.txt index a7bce12..f6547f5 100644 --- a/tests/configs/config_lasam_synth_1.txt +++ b/tests/configs/config_lasam_synth_1.txt @@ -9,4 +9,5 @@ forcing_resolution=300[sec] layer_soil_type=13,14,15 max_soil_types=25 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.0,0.0,0.0,0.0,0.0 diff --git a/tests/configs/config_lasam_synth_2.txt b/tests/configs/config_lasam_synth_2.txt index 0765a82..2559ab0 100644 --- a/tests/configs/config_lasam_synth_2.txt +++ b/tests/configs/config_lasam_synth_2.txt @@ -9,4 +9,5 @@ forcing_resolution=300[sec] layer_soil_type=13,14,15 max_soil_types=25 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.0,0.0,0.0,0.0,0.0 diff --git a/tests/configs/unittest.txt b/tests/configs/unittest.txt index 09e3554..bb46d12 100644 --- a/tests/configs/unittest.txt +++ b/tests/configs/unittest.txt @@ -6,8 +6,10 @@ initial_psi=2000.0[cm] timestep=300[sec] endtime=1.0[hr] forcing_resolution=3600[sec] +ponded_depth_max=0[cm] layer_soil_type=13,14,15 max_soil_types=15 wilting_point_psi=15495.0[cm] +field_capacity_psi=340.9[cm] giuh_ordinates=0.06,0.51,0.28,0.12,0.03 calib_params=true diff --git a/tests/main_unit_test_bmi.cxx b/tests/main_unit_test_bmi.cxx index bc2caae..84a4c56 100644 --- a/tests/main_unit_test_bmi.cxx +++ b/tests/main_unit_test_bmi.cxx @@ -453,8 +453,8 @@ int main(int argc, char *argv[]) // Benchmark values of wetting fronts depth and moisture (b is for benchmark) //std::vector depth_wf_b = {1.873813, 44.00,175.0, 200.0}; // in cm - std::vector depth_wf_b = {4.55355, 44.00,175.0, 200.0}; // in cm - std::vector theta_wf_b = {0.213716, 0.172703948143618, 0.252113867764474, 0.179593529195751}; + std::vector depth_wf_b = {4.55355239489608365, 44.00,175.0, 200.0}; // in cm + std::vector theta_wf_b = {0.21371581122514613, 0.17270389607163267, 0.25211383152603861, 0.17959348005962811}; int m_to_cm = 100; int m_to_mm = 1000; @@ -518,7 +518,7 @@ int main(int argc, char *argv[]) // check total infiltration, AET, and PET. double infiltration_check_mm = 1.896; // in mm - double AET_check_mm = 0.029801; // in mm + double AET_check_mm = 0.02980092620558239; // in mm double PET_check_mm = 0.104; // in mm double infiltration_computed = 0.0; double PET_computed = 0.0; @@ -528,7 +528,6 @@ int main(int argc, char *argv[]) model.GetValue("potential_evapotranspiration", &PET_computed); model.GetValue("actual_evapotranspiration", &AET_computed); - std::cout< 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between ponded_depth_max calibrated values set and get "<< ponded_depth_max_set<<" "<< ponded_depth_max + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(field_capacity - field_capacity_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between field_capacity calibrated values set and get "<< field_capacity_set<<" "<< field_capacity + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + for (int i=0; i < num_layers; i++) { if (fabs(smcmax[i] - smcmax_set[i]) > 1.E-5) { @@ -621,9 +645,9 @@ int main(int argc, char *argv[]) throw std::runtime_error(errMsg.str()); } - if (fabs(vg_m[i] - vg_m_set[i]) > 1.E-5) { + if (fabs(vg_n[i] - vg_n_set[i]) > 1.E-5) { std::stringstream errMsg; - errMsg << "Mismatch between vg_m calibrated values set and get "<< vg_m_set[i]<<" "<< vg_m[i] + errMsg << "Mismatch between vg_n calibrated values set and get "<< vg_n_set[i]<<" "<< vg_n[i] << " which is unexpected. \n"; throw std::runtime_error(errMsg.str()); } @@ -647,8 +671,10 @@ int main(int argc, char *argv[]) std::cout<<"| \n"; for (int i=0; i < num_layers; i++) std::cout<<"| Calib. values: layer = "<< i+1 <<", smcmax = "<< smcmax[i] - <<", vg_m = "<< vg_m[i] <<", vg_alpha = " << vg_alpha[i] + <<", vg_n = "<< vg_n[i] <<", vg_alpha = " << vg_alpha[i] <<", Ksat = "<< Ksat[i] <<"\n"; + printf("field_capacity = %lf \n", field_capacity); + printf("ponded_depth_max = %lf \n", ponded_depth_max); std::cout<<"| *************************************** \n"; std::cout<<"| LASAM Calibration test passed? YES \n"; std::cout<