From 5b64baa29e57668e181d9483e79ffa8569be5435 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 10 Jun 2024 17:07:27 -0700 Subject: [PATCH 01/32] Update BMI C++ header to match upstream with virtual dtor and correct namespace --- bmi/bmi.hxx | 8 +++++--- include/bmi_lgar.hxx | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/bmi/bmi.hxx b/bmi/bmi.hxx index 87b290b..b25fd8c 100644 --- a/bmi/bmi.hxx +++ b/bmi/bmi.hxx @@ -9,10 +9,10 @@ #include #include -namespace bmixx { +namespace bmi { - //const int BMI_SUCCESS = 0; - // const int BMI_FAILURE = 1; + const int BMI_SUCCESS = 0; + const int BMI_FAILURE = 1; const int MAX_COMPONENT_NAME = 2048; const int MAX_VAR_NAME = 2048; @@ -21,6 +21,8 @@ namespace bmixx { class Bmi { public: + virtual ~Bmi() { } + // Model control functions. virtual void Initialize(std::string config_file) = 0; virtual void Update() = 0; diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index a39fb7f..606f59d 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -27,7 +27,7 @@ class NotImplemented : public std::logic_error { } -class BmiLGAR : public bmixx::Bmi { +class BmiLGAR : public bmi::Bmi { public: BmiLGAR() { this->input_var_names[0] = "precipitation_rate"; From 01dc874996c4fcaa748a43b49d5a8cb041698a96 Mon Sep 17 00:00:00 2001 From: PeterLaFollette-NOAA <59893312+peterlafollette@users.noreply.github.com> Date: Wed, 26 Jun 2024 11:31:29 -0700 Subject: [PATCH 02/32] Ptl add adaptive timestep (#25) * adding the capability for an adaptive time step. If not provided in the config file, it defaults to false. * updating readme * updating readme * updating giuh fxn with adaptive time step * giuh now uses a fixed time step that is hardcoded to be the the same as the smallest adaptive time step for the rest of the model. If the adaptive time step is off, the giuh time step will be the same as the timestep input in the config file. * updating readme in configs * cleaned up code that was used for debugging * updating Phillipsburg config file * cleaned up code that was used for debugging * After chatting with Ahmad, we determined that it is best if the adaptive time step strategy does not alter the function in giuh.c because this function is used in multiple models. So, the adaptive time step method has been altered such that the giuh is computed at the hourly time step rather than at the sub timestep. This allows for the giuh fxn to be unaltered, and has the added bonus of not reqiring resampling for the giuh based on time step. * updating in response to comments, mostly focused on removing hardcoding from min and max time step values * more changes in response to comments, including formatting of giuh.c, removing hardcoding for internal time steps and replacing with user specified minimum time step, and establishing a maximum time step equal to the forcing resolution * removing code that was commented out --------- Co-authored-by: Peter La Follette Co-authored-by: Peter La Follette --- configs/README.md | 3 +- configs/config_lasam_Bushland.txt | 1 + configs/config_lasam_Phillipsburg.txt | 1 + configs/config_lasam_sft_ngen.txt | 1 + include/all.hxx | 2 ++ src/bmi_lgar.cxx | 49 +++++++++++++++++++-------- src/lgar.cxx | 37 +++++++++++++++----- src/soil_funcs.cxx | 2 +- 8 files changed, 72 insertions(+), 24 deletions(-) diff --git a/configs/README.md b/configs/README.md index 60ea2c4..8980fd6 100644 --- a/configs/README.md +++ b/configs/README.md @@ -12,7 +12,7 @@ A detailed description of the parameters for model configuration (i.e., initiali | initial_psi | double (scalar)| >=0 | cm | capillary head | - | used to initialize layers with a constant head | | ponded_depth_max | double (scalar)| >=0 | cm | maximum surface ponding | - | the maximum amount of water unavailable for surface drainage, default is set to zero | | timestep | double (scalar)| >0 | sec/min/hr | temporal resolution | - | timestep of the model | -| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data | +| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. | | 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) | @@ -24,3 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali | 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_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content | +| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is used as a minimum timestep (recommended value of 300 seconds). The adaptive timestep will never be larger than the forcing resolution. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. | diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt index 06b6e89..cd141ef 100644 --- a/configs/config_lasam_Bushland.txt +++ b/configs/config_lasam_Bushland.txt @@ -13,3 +13,4 @@ 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 +adaptive_timestep=true \ No newline at end of file diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt index 1d230cf..d9f3041 100644 --- a/configs/config_lasam_Phillipsburg.txt +++ b/configs/config_lasam_Phillipsburg.txt @@ -13,3 +13,4 @@ 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 +adaptive_timestep=true diff --git a/configs/config_lasam_sft_ngen.txt b/configs/config_lasam_sft_ngen.txt index e9ff131..6c22a0f 100644 --- a/configs/config_lasam_sft_ngen.txt +++ b/configs/config_lasam_sft_ngen.txt @@ -15,3 +15,4 @@ 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] +adaptive_timestep=true diff --git a/include/all.hxx b/include/all.hxx index 6b7526c..4708964 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -114,6 +114,7 @@ struct lgar_bmi_parameters double initial_psi_cm; // model initial (psi) condition double timestep_h; // model timestep in hours double forcing_resolution_h; // forcing resolution in hours + double minimum_timestep_h; // minimum time step in hours, only used if adaptive_timestep is true int forcing_interval; // = forcing_resolution_h/timestep_h int num_soil_types; // number of soil types; must be less than or equal to MAX_NUM_SOIL_TYPES double AET_cm; // actual evapotranspiration in cm @@ -131,6 +132,7 @@ struct lgar_bmi_parameters 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 */ + bool adaptive_timestep = false; // if set to true, model uses adaptive timestep. In this case, the minimum timestep is the timestep specified in the config file. The maximum time step will be equal to the forcing resolution. double ponded_depth_cm; // amount of water on the surface unavailable for surface runoff double ponded_depth_max_cm; // maximum amount of water on the surface unavailable for surface runoff double precip_previous_timestep_cm; // amount of rainfall (previous time step) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index c11a8aa..b9f92ba 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -33,18 +33,19 @@ Initialize (std::string config_file) num_giuh_ordinates = state->lgar_bmi_params.num_giuh_ordinates; - /* giuh ordinates are static and read in the lgar.cxx, and we need to have a copy of it to pass to giuh.cxx, so allocating/copying here*/ - + giuh_ordinates = new double[num_giuh_ordinates]; giuh_runoff_queue = new double[num_giuh_ordinates+1]; - for (int i=0; ilgar_bmi_params.giuh_ordinates[i+1]; // note lgar uses 1-indexing + } - for (int i=0; i<=num_giuh_ordinates;i++) + for (int i=0; i<=num_giuh_ordinates;i++){ giuh_runoff_queue[i] = 0.0; + } } @@ -76,7 +77,7 @@ Update() double mm_to_cm = 0.1; // unit conversion // local variables for readibility - int subcycles = state->lgar_bmi_params.forcing_interval; + int subcycles; int num_layers = state->lgar_bmi_params.num_layers; // local variables for a full timestep (i.e., timestep of the forcing data) @@ -109,7 +110,6 @@ Update() double volrech_subtimestep_cm; double surface_runoff_subtimestep_cm; // direct surface runoff double precip_previous_subtimestep_cm; - double volrunoff_giuh_subtimestep_cm; double volQ_gw_subtimestep_cm = 0.0; // fix it for non-zero values after adding groundwater reservoir double subtimestep_h = state->lgar_bmi_params.timestep_h; @@ -117,6 +117,7 @@ Update() 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; + bool adaptive_timestep = state->lgar_bmi_params.adaptive_timestep; // constant value used in the AET function double AET_thresh_Theta = 0.85; // scaled soil moisture (0-1) above which AET=PET (fix later!) @@ -131,6 +132,26 @@ Update() assert (state->lgar_bmi_input_params->precipitation_mm_per_h >= 0.0); assert(state->lgar_bmi_input_params->PET_mm_per_h >=0.0); + + // adaptive time step is set + if (adaptive_timestep) { + subtimestep_h = state->lgar_bmi_params.forcing_resolution_h; + if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || volon_timestep_cm > 0.0 ) { + subtimestep_h = state->lgar_bmi_params.minimum_timestep_h; //case where precip > 1 cm/h, or there is ponded head from the last time step + } + else if (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) { + subtimestep_h = state->lgar_bmi_params.minimum_timestep_h * 2.0; //case where precip is less than 1 cm/h but greater than 0, and there is no ponded head + } + subtimestep_h = fmin(subtimestep_h, state->lgar_bmi_params.forcing_resolution_h); //just in case the user has specified a minimum time step that would make the subtimestep_h greater than the forcing resolution + state->lgar_bmi_params.timestep_h = subtimestep_h; + } + + state->lgar_bmi_params.forcing_interval = int(state->lgar_bmi_params.forcing_resolution_h/state->lgar_bmi_params.timestep_h+1.0e-08); // add 1.0e-08 to prevent truncation error + subcycles = state->lgar_bmi_params.forcing_interval; + + if (verbosity.compare("high") == 0) { + printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h); + } // subcycling loop (loop over model's timestep) for (int cycle=1; cycle <= subcycles; cycle++) { @@ -358,16 +379,9 @@ Update() /*----------------------------------------------------------------------*/ - // compute giuh runoff for the subtimestep + // increment runoff for the subtimestep surface_runoff_subtimestep_cm = volrunoff_subtimestep_cm; - volrunoff_giuh_subtimestep_cm = giuh_convolution_integral(volrunoff_subtimestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue); - surface_runoff_timestep_cm += surface_runoff_subtimestep_cm ; - volrunoff_giuh_timestep_cm += volrunoff_giuh_subtimestep_cm; - - // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well. - - volQ_timestep_cm += volrunoff_giuh_subtimestep_cm; // adding groundwater flux to stream channel (note: this will be updated/corrected after adding the groundwater reservoir) volQ_gw_timestep_cm += volQ_gw_subtimestep_cm; @@ -415,6 +429,13 @@ Update() } // end of subcycling + //update giuh at the time step level (was previously updated at the sub time step level) + volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue); + + // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well. + // when groundwater component is added, it should probably happen inside of the subcycling loop. + volQ_timestep_cm = volrunoff_giuh_timestep_cm; + /*----------------------------------------------------------------------*/ // Everything related to lgar state is done at this point, now time to update some dynamic variables diff --git a/src/lgar.cxx b/src/lgar.cxx index fd412a3..815f37f 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -213,6 +213,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) // setting these options to false (defualt) state->lgar_bmi_params.sft_coupled = false; state->lgar_bmi_params.use_closed_form_G = false; + state->lgar_bmi_params.adaptive_timestep = false; bool is_layer_thickness_set = false; bool is_initial_psi_set = false; @@ -398,6 +399,20 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } + else if (param_key == "adaptive_timestep") { + if ((param_value == "false") || (param_value == "0")) { + state->lgar_bmi_params.adaptive_timestep = false; + } + else if ( (param_value == "true") || (param_value == "1")) { + state->lgar_bmi_params.adaptive_timestep = true; + } + else { + std::cerr<<"Invalid option: adaptive_timestep must be true or false. \n"; + abort(); + } + + continue; + } else if (param_key == "timestep") { state->lgar_bmi_params.timestep_h = stod(param_value); @@ -411,6 +426,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) assert (state->lgar_bmi_params.timestep_h > 0); is_timestep_set = true; + state->lgar_bmi_params.minimum_timestep_h = state->lgar_bmi_params.timestep_h; + if (verbosity.compare("high") == 0) { std::cerr<<"Model timestep [hours,seconds]: "<lgar_bmi_params.timestep_h<<" , " <lgar_bmi_params.timestep_h*3600<<"\n"; @@ -605,29 +622,28 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) throw runtime_error(errMsg.str()); } - if (is_giuh_ordinates_set) { - int factor = int(1.0/state->lgar_bmi_params.timestep_h); + int factor = int(1.0/state->lgar_bmi_params.forcing_resolution_h); state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1); state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1]; for (int i=0; ilgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor); + int index = j + i * factor + 1; + state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor); } } if (verbosity.compare("high") == 0) { for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++) - std::cerr<<"GIUH ordinates (scaled) : "<lgar_bmi_params.giuh_ordinates[i]<<"\n"; + std::cerr<<"GIUH ordinates (scaled) : "<lgar_bmi_params.giuh_ordinates[i]<<"\n"; std::cerr<<" ***** \n"; } - giuh_ordinates_temp.clear(); } + else if (!is_giuh_ordinates_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set giuh_ordinates. \n"; @@ -1272,7 +1288,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (wf_free_drainage_demand == wf) prior_mass += precip_mass_to_add - (free_drainage_demand + actual_ET_demand); // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front - double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm, delta_thetas, delta_thickness, soil_type, soil_properties); actual_ET_demand = *AET_demand_cm; @@ -2573,7 +2588,8 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm // stop the loop if the error between the current and previous psi is less than 10^-15 // 1. enough accuracy, 2. the algorithm can't improve the error further, - // 3. avoid infinite loop, 4. handles a corner case when prior mass is tiny (e.g., <1.E-5) + // 3. avoid infinite loop, 4. handles case where theta is very close to theta_r and convergence might be possible but would be extremely slow + // 5. handles a corner case when prior mass is tiny (e.g., <1.E-5) // printf("A1 = %.20f, %.18f %.18f %.18f %.18f \n ",fabs(psi_cm_loc - psi_cm_loc_prev) , psi_cm_loc, psi_cm_loc_prev, factor, delta_mass); if (fabs(psi_cm_loc - psi_cm_loc_prev) < 1E-15 && factor < 1E-13) break; @@ -2588,6 +2604,11 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm if (count_no_mass_change == break_no_mass_change) break; + if (psi_cm_loc > 1e7){//there are rare cases where theta is very close to theta_r, and delta_mass - delta_mass_prev will change extremely slowly. Convergence might be possible but the model will take hours to converge rather than seconds. + //an alternative solution was to change the threshold in if (fabs(delta_mass - delta_mass_prev) < 1e-15) to 1e-11, but that solution is somewhat slow. + break; + } + // -ve pressure will return NAN, so terminate the loop if previous psi is way small and current psi is zero // the wetting front is almost saturated if (psi_cm_loc <= 0 && psi_cm_loc_prev < 0) break; diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx index bfedb11..7da1adc 100755 --- a/src/soil_funcs.cxx +++ b/src/soil_funcs.cxx @@ -154,7 +154,7 @@ double calc_theta_from_h(double h,double alpha, double m, double n, double theta /***********************************/ double calc_Se_from_h(double h,double alpha, double m, double n) { - if(is_epsilon_less_than(h,1.0e-01)) return 1.0; // this function doesn't work well ffor tiny h + if(is_epsilon_less_than(h,1.0e-10)) return 1.0; // this function doesn't work well ffor tiny h else return(1.0/(pow(1.0+pow(alpha*h,n),m))); } From eb7292cf37be77f9a20f3d0a74cddcbe2f8c4540 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:15:24 -0600 Subject: [PATCH 03/32] feat: add listDelete for linked list --- include/all.hxx | 2 +- src/linked_list.cxx | 14 ++++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/include/all.hxx b/include/all.hxx index 4708964..e829422 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -241,7 +241,7 @@ extern void listReverseOrder(struct wetting_front** head_ref extern bool listFindLayer(struct wetting_front* link, int num_layers, double *cum_layer_thickness_cm, int *lives_in_layer, bool *extends_to_bottom_flag); extern struct wetting_front* listCopy(struct wetting_front* current, struct wetting_front* state_previous=NULL); - +extern void listDelete(struct wetting_front* head); diff --git a/src/linked_list.cxx b/src/linked_list.cxx index 926f68c..c63dbdd 100755 --- a/src/linked_list.cxx +++ b/src/linked_list.cxx @@ -26,6 +26,20 @@ //___________________________________________________________ +/*###########################################################*/ +/* listDelete() - deletes memory allocated to a linked list */ +/* This function must be called on any list to deallocate */ +/* the dynamic memory used in creating and manipultating the */ +/* list. (added by NJF) */ +/*###########################################################*/ +extern void listDelete(struct wetting_front* head) +{ + if( head != NULL ){ + struct wetting_front *next = head->next; + free( head ); + listDelete(next); + } +} /*#########################################################*/ /* listPrint() - prints a linked list to screen */ From edb6dce0025911c0ea8f4ddf44642d1373facca8 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:22:36 -0600 Subject: [PATCH 04/32] fix: deallocate memory with listDelete before overwriting with listCopy --- src/bmi_lgar.cxx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index b9f92ba..397291e 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -164,7 +164,10 @@ Update() std::cerr<<"BMI Update |Timesteps = "<< state->lgar_bmi_params.timesteps<<", Time [h] = "<state->lgar_bmi_params.time_s / 3600.<<", Subcycle = "<< cycle <<" of "<state_previous = NULL; + if( state->state_previous != NULL ){ + listDelete(state->state_previous); + state->state_previous = NULL; + } state->state_previous = listCopy(state->head); // ensure precip and PET are non-negative @@ -291,7 +294,10 @@ Update() listPrint(state->head); } - state->state_previous = NULL; + if(state->state_previous != NULL ){ + listDelete(state->state_previous); + state->state_previous = NULL; + } state->state_previous = listCopy(state->head); volin_timestep_cm += volin_subtimestep_cm; From ce5d7c5b642c53307781bbb77b515ca1daa7a5f6 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:18:09 -0600 Subject: [PATCH 05/32] fix: add destructor to deallocate class held pointers --- include/bmi_lgar.hxx | 1 + src/bmi_lgar.cxx | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index 606f59d..abcec7c 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -29,6 +29,7 @@ class NotImplemented : public std::logic_error { class BmiLGAR : public bmi::Bmi { public: + ~BmiLGAR(); BmiLGAR() { this->input_var_names[0] = "precipitation_rate"; this->input_var_names[1] = "potential_evapotranspiration_rate"; diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 397291e..e9b488e 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -17,6 +17,14 @@ // default verbosity is set to 'none' other option 'high' or 'low' needs to be specified in the config file string verbosity="none"; +/** + * @brief Delete dynamic arrays allocated in Initialize() and held by this object + * + */ +BmiLGAR::~BmiLGAR(){ + if( giuh_ordinates != nullptr ) delete [] giuh_ordinates; + if( giuh_runoff_queue != nullptr ) delete [] giuh_runoff_queue; +} /* The `head` pointer stores the address in memory of the first member of the linked list containing all the wetting fronts. The contents of struct wetting_front are defined in "all.h" */ From 0dc6980c70d7268a57dfaee237ab06b33c02eed9 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:20:20 -0600 Subject: [PATCH 06/32] fix: ensure cleanup of old memory when allocating soil params --- include/bmi_lgar.hxx | 1 + src/bmi_lgar.cxx | 17 +++++++++++++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index abcec7c..741b419 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -132,6 +132,7 @@ public: struct model_state* get_model(); private: + void realloc_soil(); struct model_state* state; static const int input_var_name_count = 3; static const int output_var_name_count = 15; diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index e9b488e..b5bec52 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -57,6 +57,20 @@ Initialize (std::string config_file) } +/** + * @brief Allocate (or reallocate) storage for soil parameters + * + */ +void BmiLGAR::realloc_soil(){ + if(state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr) + delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; + if(state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) + delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; + + state->lgar_bmi_params.soil_depth_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; + state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; +} + /* This is the main function calling lgar subroutines for creating, moving, and merging wetting fronts. Calls to AET and mass balance module are also happening here @@ -457,8 +471,7 @@ Update() state->lgar_bmi_params.num_wetting_fronts = listLength(state->head); // allocate new memory based on updated wetting fronts; we could make it conditional i.e. create only if no. of wf are changed - state->lgar_bmi_params.soil_depth_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; - state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; + realloc_soil(); // update thickness/depth and soil moisture of wetting fronts (used for state coupling) struct wetting_front *current = state->head; From 8a76135fd6d0be107043f655a42ca14fafa60010 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:25:09 -0600 Subject: [PATCH 07/32] fix: delete memory when deleting front from list --- src/linked_list.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linked_list.cxx b/src/linked_list.cxx index c63dbdd..3782d12 100755 --- a/src/linked_list.cxx +++ b/src/linked_list.cxx @@ -254,7 +254,7 @@ extern struct wetting_front* listDeleteFront(int front_num, struct wetting_front previous = current->next; } - + if( current != NULL ) free( current ); current = previous; while(previous != NULL) { // decrement all front numbers From b74f03084407a35860aaf06155da059b5f34e51a Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:27:04 -0600 Subject: [PATCH 08/32] fix: delete tempory heap allocated variables in bmi_main --- src/bmi_main_lgar.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx index 1c64c93..37b0ac9 100644 --- a/src/bmi_main_lgar.cxx +++ b/src/bmi_main_lgar.cxx @@ -155,6 +155,8 @@ int main(int argc, char *argv[]) // write layers data to file fprintf(outlayer_fptr,"# Timestep = %d, %s \n", i, time[i].c_str()); write_state(outlayer_fptr, model_state.get_model()->head); + delete [] soil_moisture_wetting_front; + delete [] soil_thickness_wetting_front; } } From 7fd7c3b4f5cb439d04c4727821fe2205dbcb238a Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:33:07 -0600 Subject: [PATCH 09/32] fix: delete temporary dynamic allocated arrays --- src/lgar.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/lgar.cxx b/src/lgar.cxx index 815f37f..f57de87 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -1150,7 +1150,9 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm, delta_thetas, delta_thickness, soil_type, soil_properties); actual_ET_demand = *AET_demand_cm; - + //done with delta_thetas and delta_thickness, cleanup memory + delete delta_thetas; + delete delta_thickness; current->theta = fmax(theta_r, fmin(theta_new, theta_e)); double Se = calc_Se_from_theta(current->theta,theta_e,theta_r); From f1ab3a41c92fb566938c1b126a0b2f912fbbcecb Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:28:59 -0600 Subject: [PATCH 10/32] fix: clean up list memory prior to listCopy in main.cxx --- src/main.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/main.cxx b/src/main.cxx index 5c9e586..8f15113 100755 --- a/src/main.cxx +++ b/src/main.cxx @@ -435,6 +435,7 @@ for(time_step_num=0;time_step_num Date: Wed, 12 Jun 2024 19:29:55 -0600 Subject: [PATCH 11/32] fix: cleanup dynamic memory in bmi state during Finalize --- src/bmi_lgar.cxx | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index b5bec52..e3fae68 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -631,6 +631,30 @@ void BmiLGAR:: Finalize() { global_mass_balance(); + if( state->head != NULL ) listDelete(state->head); + if( state->state_previous != NULL ) listDelete(state->state_previous); + + if( state->soil_properties != nullptr ) delete [] state->soil_properties; + + if( state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr ) delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; + if( state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; + + if( state->lgar_bmi_params.soil_temperature != nullptr ) delete [] state->lgar_bmi_params.soil_temperature; + if( state->lgar_bmi_params.soil_temperature_z != nullptr) delete [] state->lgar_bmi_params.soil_temperature_z; + if( state->lgar_bmi_params.layer_soil_type != nullptr ) delete [] state->lgar_bmi_params.layer_soil_type; + + if( state->lgar_calib_params.theta_e != nullptr ) delete [] state->lgar_calib_params.theta_e; + if( state->lgar_calib_params.theta_r != nullptr ) delete [] state->lgar_calib_params.theta_r; + if( state->lgar_calib_params.vg_n != nullptr ) delete [] state->lgar_calib_params.vg_n; + if( state->lgar_calib_params.vg_alpha != nullptr ) delete [] state->lgar_calib_params.vg_alpha; + if( state->lgar_calib_params.Ksat != nullptr ) delete [] state->lgar_calib_params.Ksat; + + if( state->lgar_bmi_params.layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.layer_thickness_cm; + if( state->lgar_bmi_params.cum_layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.cum_layer_thickness_cm; + if( state->lgar_bmi_params.giuh_ordinates != nullptr ) delete [] state->lgar_bmi_params.giuh_ordinates; + if( state->lgar_bmi_params.frozen_factor != nullptr ) delete [] state->lgar_bmi_params.frozen_factor; + if( state->lgar_bmi_input_params != nullptr ) delete state->lgar_bmi_input_params; + if( state != nullptr ) delete state; } From c4a942f6e68c733cf648c41885c4ea678f6e62bc Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:27:36 -0600 Subject: [PATCH 12/32] fix: call Finalize on bmi model --- src/bmi_main_lgar.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx index 37b0ac9..26c0ed8 100644 --- a/src/bmi_main_lgar.cxx +++ b/src/bmi_main_lgar.cxx @@ -163,7 +163,7 @@ int main(int argc, char *argv[]) // do final mass balance model_state.global_mass_balance(); - + model_state.Finalize(); if (outdata_fptr) { fclose(outdata_fptr); fclose(outlayer_fptr); From 6032f640f0a53efffc01ce96a2ba9eae20cc96a4 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 19:32:03 -0600 Subject: [PATCH 13/32] chore: inline comments on things that don't make sense --- src/lgar.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/lgar.cxx b/src/lgar.cxx index f57de87..ad2ff63 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -659,6 +659,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } } else { + //NJF FIXME these arrays should be allocated based on num_cells_temp... state->lgar_bmi_params.soil_temperature = new double[1](); state->lgar_bmi_params.soil_temperature_z = new double[1](); state->lgar_bmi_params.num_cells_temp = 1; @@ -1133,6 +1134,8 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf double theta_below = 0.0; new_mass += layer_thickness * (theta - theta_below); + //NJF theta_below is always 0, so all delta_thetas are always 0... + //does this really need a dynamic array in this case??? delta_thetas[k] = theta_below; delta_thickness[k] = layer_thickness; } From 6e9f0a65abd7205a682202f0c61cbe21f2f80c49 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 12 Jun 2024 20:51:20 -0600 Subject: [PATCH 14/32] fix: free instead of delete for malloc'd mem --- src/lgar.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lgar.cxx b/src/lgar.cxx index ad2ff63..bf59dcf 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -1154,8 +1154,8 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf delta_thetas, delta_thickness, soil_type, soil_properties); actual_ET_demand = *AET_demand_cm; //done with delta_thetas and delta_thickness, cleanup memory - delete delta_thetas; - delete delta_thickness; + free(delta_thetas); + free(delta_thickness); current->theta = fmax(theta_r, fmin(theta_new, theta_e)); double Se = calc_Se_from_theta(current->theta,theta_e,theta_r); From 15a0842912ab1ddeeb32443e40c33e6d85c55736 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Thu, 13 Jun 2024 08:39:34 -0600 Subject: [PATCH 15/32] fix: free another set of tempory allocated pointers --- src/lgar.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/lgar.cxx b/src/lgar.cxx index bf59dcf..f5a497d 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -1296,7 +1296,9 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm, delta_thetas, delta_thickness, soil_type, soil_properties); actual_ET_demand = *AET_demand_cm; - + //done with delta_thetas and delta_thickness, cleanup memory + free(delta_thetas); + free(delta_thickness); current->theta = fmax(theta_r, fmin(theta_new, theta_e)); } From a15801fb4d1c35dc0daccf3ec19afb563af858ca Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Thu, 13 Jun 2024 12:27:33 -0600 Subject: [PATCH 16/32] fix: explicitly intialize class pointers to nullptr --- include/bmi_lgar.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index 741b419..b29339a 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -30,7 +30,7 @@ class NotImplemented : public std::logic_error { class BmiLGAR : public bmi::Bmi { public: ~BmiLGAR(); - BmiLGAR() { + BmiLGAR():giuh_ordinates(nullptr), giuh_runoff_queue(nullptr) { this->input_var_names[0] = "precipitation_rate"; this->input_var_names[1] = "potential_evapotranspiration_rate"; this->input_var_names[2] = "soil_temperature_profile"; From 00420a50f4dc0f933a4e8a9ae86712e5606934f5 Mon Sep 17 00:00:00 2001 From: Nels Date: Wed, 10 Jul 2024 13:16:21 -0700 Subject: [PATCH 17/32] chore: change recursive call to loop in listDelete Co-authored-by: Phil Miller - NOAA --- src/linked_list.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linked_list.cxx b/src/linked_list.cxx index 3782d12..0176d41 100755 --- a/src/linked_list.cxx +++ b/src/linked_list.cxx @@ -34,10 +34,10 @@ /*###########################################################*/ extern void listDelete(struct wetting_front* head) { - if( head != NULL ){ + while (head != NULL) { struct wetting_front *next = head->next; free( head ); - listDelete(next); + head = next; } } From 9a2c0418b119dccb4c33c7bc1fe412023872b5d1 Mon Sep 17 00:00:00 2001 From: Nels Date: Wed, 10 Jul 2024 13:20:13 -0700 Subject: [PATCH 18/32] chore: remove redundant NULL check Co-authored-by: Phil Miller - NOAA --- src/main.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/main.cxx b/src/main.cxx index 8f15113..c3e51dd 100755 --- a/src/main.cxx +++ b/src/main.cxx @@ -435,8 +435,7 @@ for(time_step_num=0;time_step_num Date: Wed, 10 Jul 2024 14:22:46 -0600 Subject: [PATCH 19/32] fix: remove extra mass balance check at end of bmi_main test --- src/bmi_main_lgar.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx index 26c0ed8..a106c8b 100644 --- a/src/bmi_main_lgar.cxx +++ b/src/bmi_main_lgar.cxx @@ -161,8 +161,7 @@ int main(int argc, char *argv[]) } - // do final mass balance - model_state.global_mass_balance(); + // do final mass balance ( inside Finalize() ) and finish the simulation model_state.Finalize(); if (outdata_fptr) { fclose(outdata_fptr); From a00443160570a2f907f9744a1739eb8173912574 Mon Sep 17 00:00:00 2001 From: hellkite500 Date: Wed, 10 Jul 2024 14:31:35 -0600 Subject: [PATCH 20/32] chore: remove defensive checks for nullptr --- src/bmi_lgar.cxx | 59 ++++++++++++++++++++++++------------------------ 1 file changed, 29 insertions(+), 30 deletions(-) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index e3fae68..f1b744c 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -22,8 +22,8 @@ string verbosity="none"; * */ BmiLGAR::~BmiLGAR(){ - if( giuh_ordinates != nullptr ) delete [] giuh_ordinates; - if( giuh_runoff_queue != nullptr ) delete [] giuh_runoff_queue; + delete [] giuh_ordinates; + delete [] giuh_runoff_queue; } /* The `head` pointer stores the address in memory of the first member of the linked list containing @@ -62,10 +62,9 @@ Initialize (std::string config_file) * */ void BmiLGAR::realloc_soil(){ - if(state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr) - delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; - if(state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) - delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; + + delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; + delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; state->lgar_bmi_params.soil_depth_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts]; @@ -631,30 +630,30 @@ void BmiLGAR:: Finalize() { global_mass_balance(); - if( state->head != NULL ) listDelete(state->head); - if( state->state_previous != NULL ) listDelete(state->state_previous); - - if( state->soil_properties != nullptr ) delete [] state->soil_properties; - - if( state->lgar_bmi_params.soil_depth_wetting_fronts != nullptr ) delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; - if( state->lgar_bmi_params.soil_moisture_wetting_fronts != nullptr) delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; - - if( state->lgar_bmi_params.soil_temperature != nullptr ) delete [] state->lgar_bmi_params.soil_temperature; - if( state->lgar_bmi_params.soil_temperature_z != nullptr) delete [] state->lgar_bmi_params.soil_temperature_z; - if( state->lgar_bmi_params.layer_soil_type != nullptr ) delete [] state->lgar_bmi_params.layer_soil_type; - - if( state->lgar_calib_params.theta_e != nullptr ) delete [] state->lgar_calib_params.theta_e; - if( state->lgar_calib_params.theta_r != nullptr ) delete [] state->lgar_calib_params.theta_r; - if( state->lgar_calib_params.vg_n != nullptr ) delete [] state->lgar_calib_params.vg_n; - if( state->lgar_calib_params.vg_alpha != nullptr ) delete [] state->lgar_calib_params.vg_alpha; - if( state->lgar_calib_params.Ksat != nullptr ) delete [] state->lgar_calib_params.Ksat; - - if( state->lgar_bmi_params.layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.layer_thickness_cm; - if( state->lgar_bmi_params.cum_layer_thickness_cm != nullptr ) delete [] state->lgar_bmi_params.cum_layer_thickness_cm; - if( state->lgar_bmi_params.giuh_ordinates != nullptr ) delete [] state->lgar_bmi_params.giuh_ordinates; - if( state->lgar_bmi_params.frozen_factor != nullptr ) delete [] state->lgar_bmi_params.frozen_factor; - if( state->lgar_bmi_input_params != nullptr ) delete state->lgar_bmi_input_params; - if( state != nullptr ) delete state; + listDelete(state->head); + listDelete(state->state_previous); + + delete [] state->soil_properties; + + delete [] state->lgar_bmi_params.soil_depth_wetting_fronts; + delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts; + + delete [] state->lgar_bmi_params.soil_temperature; + delete [] state->lgar_bmi_params.soil_temperature_z; + delete [] state->lgar_bmi_params.layer_soil_type; + + delete [] state->lgar_calib_params.theta_e; + delete [] state->lgar_calib_params.theta_r; + delete [] state->lgar_calib_params.vg_n; + delete [] state->lgar_calib_params.vg_alpha; + delete [] state->lgar_calib_params.Ksat; + + delete [] state->lgar_bmi_params.layer_thickness_cm; + delete [] state->lgar_bmi_params.cum_layer_thickness_cm; + delete [] state->lgar_bmi_params.giuh_ordinates; + delete [] state->lgar_bmi_params.frozen_factor; + delete state->lgar_bmi_input_params; + delete state; } From 567adc75bea3ed3580427592f198638b0a8a32ef Mon Sep 17 00:00:00 2001 From: Austin Raney Date: Wed, 11 Sep 2024 13:59:46 -0400 Subject: [PATCH 21/32] chore: add soil_params_file ordered by stat categories in nom soil table (#30) vG_default_params.dat re-ordered using stat categories like in nom https://github.com/NOAA-OWP/noah-owp-modular/blob/0abb891b48b043cc626c4e4bbd0efe54ad357fe1/parameters/SOILPARM.TBL --- data/vG_params_stat_nom_ordered.dat | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 data/vG_params_stat_nom_ordered.dat diff --git a/data/vG_params_stat_nom_ordered.dat b/data/vG_params_stat_nom_ordered.dat new file mode 100644 index 0000000..80c8932 --- /dev/null +++ b/data/vG_params_stat_nom_ordered.dat @@ -0,0 +1,13 @@ +"Texture theta_r theta_e alpha (cm^-1) n Ks (cm/h)" +"Sand" 0.05 0.38 4.00E-02 3.18 26.64 +"Loamy-sand" 0.05 0.39 3.00E-02 1.75 4.32 +"Sandy-loam" 0.04 0.39 3.00E-02 1.45 1.584 +"Silt-loam" 0.07 0.44 1.00E-02 1.66 0.756 +"Silt" 0.05 0.49 1.00E-02 1.68 1.836 +"Loam" 0.06 0.4 1.00E-02 1.47 0.504 +"Sandy-clay loam" 0.06 0.38 2.00E-02 1.33 0.54 +"Silty-clay loam" 0.09 0.48 1.00E-02 1.52 0.468 +"Clay-loam" 0.08 0.44 2.00E-02 1.42 0.3348 +"Sandy-clay" 0.12 0.39 3.00E-02 1.21 0.468 +"Silty-clay" 0.11 0.48 2.00E-02 1.32 0.432 +"Clay" 0.1 0.46 1.00E-02 1.25 0.612 From 5120df42f6eb0dea32aef4659daaaa7abb1e53d8 Mon Sep 17 00:00:00 2001 From: Austin Raney Date: Thu, 12 Sep 2024 10:01:15 -0400 Subject: [PATCH 22/32] fix: seek to beginning of config file after looking for verbose key (#33) * fix: seek to beginning of config file after looking for verbose key; fixes #32 * fix: clear ifstream flags before seeking Co-authored-by: Nels --------- Co-authored-by: Nels --- src/lgar.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/lgar.cxx b/src/lgar.cxx index f5a497d..5c02def 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -205,6 +205,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } } + // seek to beginning of input after searching for 'verbosity' + fp.clear(); + fp.seekg(0, fp.beg); if (verbosity.compare("none") != 0) { std::cerr<<"------------- Initialization from config file ---------------------- \n"; @@ -2646,4 +2649,4 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm } -#endif \ No newline at end of file +#endif From 1f8d4f8b6f83b7636df7541b0623d352cd999274 Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Wed, 11 Sep 2024 13:29:05 -0400 Subject: [PATCH 23/32] modification to handle invalid soil_type such as waterbody, glacier, etc., the model returns input precip as discharge. --- configs/README.md | 2 +- configs/config_lasam_Bushland.txt | 4 +-- configs/config_lasam_Phillipsburg.txt | 2 +- configs/config_lasam_sft_ngen.txt | 2 +- include/all.hxx | 3 ++- src/bmi_lgar.cxx | 36 +++++++++++++++++++++++++-- src/lgar.cxx | 16 +++++++++--- 7 files changed, 53 insertions(+), 12 deletions(-) diff --git a/configs/README.md b/configs/README.md index 8980fd6..d783909 100644 --- a/configs/README.md +++ b/configs/README.md @@ -15,7 +15,7 @@ A detailed description of the parameters for model configuration (i.e., initiali | forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. | | 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) | +| max_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | | 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. | diff --git a/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt index cd141ef..4f3877c 100644 --- a/configs/config_lasam_Bushland.txt +++ b/configs/config_lasam_Bushland.txt @@ -9,8 +9,8 @@ forcing_resolution=3600[sec] ponded_depth_max=0[cm] use_closed_form_G=false layer_soil_type=16,17,18 -max_soil_types=25 +max_valid_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 -adaptive_timestep=true \ No newline at end of file +adaptive_timestep=true diff --git a/configs/config_lasam_Phillipsburg.txt b/configs/config_lasam_Phillipsburg.txt index d9f3041..4044ed2 100644 --- a/configs/config_lasam_Phillipsburg.txt +++ b/configs/config_lasam_Phillipsburg.txt @@ -9,7 +9,7 @@ forcing_resolution=3600[sec] ponded_depth_max=2[cm] use_closed_form_G=false layer_soil_type=13,14,15 -max_soil_types=25 +max_valid_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 6c22a0f..1ce014d 100644 --- a/configs/config_lasam_sft_ngen.txt +++ b/configs/config_lasam_sft_ngen.txt @@ -9,7 +9,7 @@ forcing_resolution=3600[sec] ponded_depth_max=2[cm] use_closed_form_G=false layer_soil_type=13,14,15 -max_soil_types=25 +max_valid_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/include/all.hxx b/include/all.hxx index e829422..52750f0 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -146,7 +146,8 @@ struct lgar_bmi_parameters double *giuh_ordinates; // geomorphological instantaneous unit hydrograph int num_giuh_ordinates; // number of giuh ordinates - int calib_params_flag = 0; // flag for calibratable parameters; if true, then calibratable params are updated otherwise not + int calib_params_flag = 0; // flag for calibratable parameters; if true, then calibratable params are updated otherwise not + bool is_invalid_soil_type = true; // checks if the provided soil type is valid for the model, if not then return Q_out = precip }; // Define a data structure for local (timestep) and global mass balance parameters diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index f1b744c..7e3d1bd 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -83,7 +83,40 @@ Update() std::cerr<<"|****************** LASAM BMI Update... ******************|\n"; std::cerr<<"---------------------------------------------------------\n"; } - + + double mm_to_cm = 0.1; // unit conversion + + if (state->lgar_bmi_params.is_invalid_soil_type) { + // add to mass balance accumulated variables + state->lgar_mass_balance.volprecip_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + state->lgar_mass_balance.volin_cm += 0.0; + state->lgar_mass_balance.volon_cm = 0.0; + state->lgar_mass_balance.volend_cm = state->lgar_mass_balance.volstart_cm;; + state->lgar_mass_balance.volAET_cm += 0.0; + state->lgar_mass_balance.volrech_cm += 0.0; + state->lgar_mass_balance.volrunoff_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + state->lgar_mass_balance.volQ_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + state->lgar_mass_balance.volQ_gw_cm += 0.0; + state->lgar_mass_balance.volPET_cm += 0.0; + state->lgar_mass_balance.volrunoff_giuh_cm += 0.0; + state->lgar_mass_balance.volchange_calib_cm += 0.0; + + // converted values, a struct local to the BMI and has bmi output variables + bmi_unit_conv.mass_balance_m = state->lgar_mass_balance.local_mass_balance * state->units.cm_to_m; + bmi_unit_conv.volprecip_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + bmi_unit_conv.volin_timestep_m = 0.0; + bmi_unit_conv.volend_timestep_m = 0.0; + bmi_unit_conv.volAET_timestep_m = 0.0; + bmi_unit_conv.volrech_timestep_m = 0.0; + bmi_unit_conv.volrunoff_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + bmi_unit_conv.volQ_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + bmi_unit_conv.volQ_gw_timestep_m = 0.0; + bmi_unit_conv.volPET_timestep_m = 0.0; + bmi_unit_conv.volrunoff_giuh_timestep_m = 0.0; + + return; + } + // if lasam is coupled to soil freeze-thaw, frozen fraction module is called if (state->lgar_bmi_params.sft_coupled) frozen_factor_hydraulic_conductivity(state->lgar_bmi_params); @@ -95,7 +128,6 @@ Update() state->lgar_bmi_params.calib_params_flag = false; } - double mm_to_cm = 0.1; // unit conversion // local variables for readibility int subcycles; diff --git a/src/lgar.cxx b/src/lgar.cxx index 5c02def..0df5315 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -350,7 +350,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } - else if (param_key == "max_soil_types") { + else if (param_key == "max_valid_soil_types") { state->lgar_bmi_params.num_soil_types = stod(param_value); is_max_soil_types_set = true; continue; @@ -537,7 +537,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } if(!is_max_soil_types_set) - state->lgar_bmi_params.num_soil_types = 15; // maximum number of soil types defaults to 15 + state->lgar_bmi_params.num_soil_types = 12; // maximum number of soil types defaults to 12 if (verbosity.compare("high") == 0) { std::cerr<<"Maximum number of soil types: "<lgar_bmi_params.num_soil_types<<"\n"; @@ -563,9 +563,17 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) wilting_point_psi_cm, state->soil_properties); // check if soil layers provided are within the range + state->lgar_bmi_params.is_invalid_soil_type = false; // model not valid for soil types = waterbody, glacier, lava, etc. for (int layer=1; layer <= state->lgar_bmi_params.num_layers; layer++) { - assert (state->lgar_bmi_params.layer_soil_type[layer] <= state->lgar_bmi_params.num_soil_types); - assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file); + //assert (state->lgar_bmi_params.layer_soil_type[layer] <= state->lgar_bmi_params.num_soil_types); + //assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file); + if (state->lgar_bmi_params.layer_soil_type[layer] > state->lgar_bmi_params.num_soil_types) { + state->lgar_bmi_params.is_invalid_soil_type = true; + std::cerr << "Invalid soil type: " + << state->lgar_bmi_params.layer_soil_type[layer] + <<". Model returns input_precip = ouput_Qout. \n"; + break; + } } if (verbosity.compare("high") == 0) { From 5f2269dc79dd05ed21617d0f1a334aacbc1d33ea Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Wed, 11 Sep 2024 13:41:41 -0400 Subject: [PATCH 24/32] update tests/config files, and replace max_soil_type with max_valid_soil_type. --- src/lgar.cxx | 6 +++--- tests/configs/config_lasam_synth_0.txt | 4 ++-- tests/configs/config_lasam_synth_1.txt | 2 +- tests/configs/config_lasam_synth_2.txt | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/lgar.cxx b/src/lgar.cxx index 0df5315..7bcd9d2 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -227,7 +227,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) 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_max_valid_soil_types_set = false; bool is_giuh_ordinates_set = false; bool is_soil_z_set = false; bool is_ponded_depth_max_cm_set = false; @@ -352,7 +352,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } else if (param_key == "max_valid_soil_types") { state->lgar_bmi_params.num_soil_types = stod(param_value); - is_max_soil_types_set = true; + is_max_valid_soil_types_set = true; continue; } else if (param_key == "soil_params_file") { @@ -536,7 +536,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) std::cerr<<" ***** \n"; } - if(!is_max_soil_types_set) + if(!is_max_valid_soil_types_set) state->lgar_bmi_params.num_soil_types = 12; // maximum number of soil types defaults to 12 if (verbosity.compare("high") == 0) { diff --git a/tests/configs/config_lasam_synth_0.txt b/tests/configs/config_lasam_synth_0.txt index a359360..801de13 100644 --- a/tests/configs/config_lasam_synth_0.txt +++ b/tests/configs/config_lasam_synth_0.txt @@ -7,8 +7,8 @@ timestep=300[sec] endtime=12[h] forcing_resolution=3600[sec] layer_soil_type=13,14,15 -max_soil_types=15 +max_valid_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] \ No newline at end of file +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 f6547f5..761d5c0 100644 --- a/tests/configs/config_lasam_synth_1.txt +++ b/tests/configs/config_lasam_synth_1.txt @@ -7,7 +7,7 @@ timestep=300[sec] endtime=12[hr] forcing_resolution=300[sec] layer_soil_type=13,14,15 -max_soil_types=25 +max_valid_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 2559ab0..e3b1da8 100644 --- a/tests/configs/config_lasam_synth_2.txt +++ b/tests/configs/config_lasam_synth_2.txt @@ -7,7 +7,7 @@ timestep=300[sec] endtime=12[hr] forcing_resolution=300[sec] layer_soil_type=13,14,15 -max_soil_types=25 +max_valid_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 From 634b5764e4a7a2eb536fc60ca83caf3b3ba0ac5e Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Wed, 11 Sep 2024 13:47:28 -0400 Subject: [PATCH 25/32] updated unittest config file. --- tests/configs/unittest.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/configs/unittest.txt b/tests/configs/unittest.txt index bb46d12..118d554 100644 --- a/tests/configs/unittest.txt +++ b/tests/configs/unittest.txt @@ -8,7 +8,7 @@ endtime=1.0[hr] forcing_resolution=3600[sec] ponded_depth_max=0[cm] layer_soil_type=13,14,15 -max_soil_types=15 +max_valid_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 From 80661bd133461a3585e3d45152ab9d846e396f22 Mon Sep 17 00:00:00 2001 From: AJKhattak-NOAA Date: Wed, 11 Sep 2024 18:17:27 -0400 Subject: [PATCH 26/32] Update src/bmi_lgar.cxx Co-authored-by: Austin Raney --- src/bmi_lgar.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 7e3d1bd..c76235d 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -91,7 +91,7 @@ Update() state->lgar_mass_balance.volprecip_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; state->lgar_mass_balance.volin_cm += 0.0; state->lgar_mass_balance.volon_cm = 0.0; - state->lgar_mass_balance.volend_cm = state->lgar_mass_balance.volstart_cm;; + state->lgar_mass_balance.volend_cm = state->lgar_mass_balance.volstart_cm; state->lgar_mass_balance.volAET_cm += 0.0; state->lgar_mass_balance.volrech_cm += 0.0; state->lgar_mass_balance.volrunoff_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; From 6aab34624daacfeeb41f82503fd211e41f492d63 Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Wed, 11 Sep 2024 18:54:43 -0400 Subject: [PATCH 27/32] set max limit on the max_valid_soil_types (25 soils), removed +=0. --- configs/README.md | 2 +- src/bmi_lgar.cxx | 14 +++++++------- src/lgar.cxx | 2 +- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/configs/README.md b/configs/README.md index d783909..616e764 100644 --- a/configs/README.md +++ b/configs/README.md @@ -15,7 +15,7 @@ A detailed description of the parameters for model configuration (i.e., initiali | forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. | | 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_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | +| max_valid_soil_types | int | >= 1 and <= 25 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | | 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. | diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index c76235d..6580b51 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -89,17 +89,17 @@ Update() if (state->lgar_bmi_params.is_invalid_soil_type) { // add to mass balance accumulated variables state->lgar_mass_balance.volprecip_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; - state->lgar_mass_balance.volin_cm += 0.0; + state->lgar_mass_balance.volin_cm = 0.0; state->lgar_mass_balance.volon_cm = 0.0; state->lgar_mass_balance.volend_cm = state->lgar_mass_balance.volstart_cm; - state->lgar_mass_balance.volAET_cm += 0.0; - state->lgar_mass_balance.volrech_cm += 0.0; + state->lgar_mass_balance.volAET_cm = 0.0; + state->lgar_mass_balance.volrech_cm = 0.0; state->lgar_mass_balance.volrunoff_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; state->lgar_mass_balance.volQ_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; - state->lgar_mass_balance.volQ_gw_cm += 0.0; - state->lgar_mass_balance.volPET_cm += 0.0; - state->lgar_mass_balance.volrunoff_giuh_cm += 0.0; - state->lgar_mass_balance.volchange_calib_cm += 0.0; + state->lgar_mass_balance.volQ_gw_cm = 0.0; + state->lgar_mass_balance.volPET_cm = 0.0; + state->lgar_mass_balance.volrunoff_giuh_cm = 0.0; + state->lgar_mass_balance.volchange_calib_cm = 0.0; // converted values, a struct local to the BMI and has bmi output variables bmi_unit_conv.mass_balance_m = state->lgar_mass_balance.local_mass_balance * state->units.cm_to_m; diff --git a/src/lgar.cxx b/src/lgar.cxx index 7bcd9d2..85543cf 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -351,7 +351,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } else if (param_key == "max_valid_soil_types") { - state->lgar_bmi_params.num_soil_types = stod(param_value); + state->lgar_bmi_params.num_soil_types = fmax(stod(param_value), 25.0); is_max_valid_soil_types_set = true; continue; } From d768ac3cdb9ec9e39b9573b9c0fbdce50e0186b6 Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Thu, 12 Sep 2024 10:34:41 -0400 Subject: [PATCH 28/32] changed stod -> stoi for num_soil_types, removed 25 that was added as an upper limit on num_soil_type as the default is already there. --- configs/README.md | 2 +- src/lgar.cxx | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/configs/README.md b/configs/README.md index 616e764..d783909 100644 --- a/configs/README.md +++ b/configs/README.md @@ -15,7 +15,7 @@ A detailed description of the parameters for model configuration (i.e., initiali | forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. | | 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_valid_soil_types | int | >= 1 and <= 25 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | +| max_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | | 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. | diff --git a/src/lgar.cxx b/src/lgar.cxx index 85543cf..a4b293c 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -351,7 +351,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } else if (param_key == "max_valid_soil_types") { - state->lgar_bmi_params.num_soil_types = fmax(stod(param_value), 25.0); + state->lgar_bmi_params.num_soil_types = stoi(param_value); is_max_valid_soil_types_set = true; continue; } @@ -569,9 +569,11 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) //assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file); if (state->lgar_bmi_params.layer_soil_type[layer] > state->lgar_bmi_params.num_soil_types) { state->lgar_bmi_params.is_invalid_soil_type = true; - std::cerr << "Invalid soil type: " - << state->lgar_bmi_params.layer_soil_type[layer] - <<". Model returns input_precip = ouput_Qout. \n"; + if (verbosity.compare("high") == 0) { + std::cerr << "Invalid soil type: " + << state->lgar_bmi_params.layer_soil_type[layer] + <<". Model returns input_precip = ouput_Qout. \n"; + } break; } } From 2eddeea3544e18ba6696ba24ae4e00aa31af6902 Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Thu, 12 Sep 2024 11:58:16 -0400 Subject: [PATCH 29/32] unit correction (for invalid soiltype outputs). --- src/bmi_lgar.cxx | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 6580b51..38e42a1 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -85,7 +85,8 @@ Update() } double mm_to_cm = 0.1; // unit conversion - + double mm_to_m = 0.001; + if (state->lgar_bmi_params.is_invalid_soil_type) { // add to mass balance accumulated variables state->lgar_mass_balance.volprecip_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; @@ -102,14 +103,14 @@ Update() state->lgar_mass_balance.volchange_calib_cm = 0.0; // converted values, a struct local to the BMI and has bmi output variables - bmi_unit_conv.mass_balance_m = state->lgar_mass_balance.local_mass_balance * state->units.cm_to_m; - bmi_unit_conv.volprecip_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + bmi_unit_conv.mass_balance_m = 0.0; + bmi_unit_conv.volprecip_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_m; bmi_unit_conv.volin_timestep_m = 0.0; bmi_unit_conv.volend_timestep_m = 0.0; bmi_unit_conv.volAET_timestep_m = 0.0; bmi_unit_conv.volrech_timestep_m = 0.0; - bmi_unit_conv.volrunoff_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; - bmi_unit_conv.volQ_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; + bmi_unit_conv.volrunoff_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_m; + bmi_unit_conv.volQ_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_m; bmi_unit_conv.volQ_gw_timestep_m = 0.0; bmi_unit_conv.volPET_timestep_m = 0.0; bmi_unit_conv.volrunoff_giuh_timestep_m = 0.0; From 46aa6236541c3ff837296016c706c9e1a1053436 Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Thu, 12 Sep 2024 12:53:44 -0400 Subject: [PATCH 30/32] set MAX_NUM_SOIL_TYPE as both the upper limit the default value for num_soil_types. --- configs/README.md | 2 +- src/lgar.cxx | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/configs/README.md b/configs/README.md index d783909..4b839ee 100644 --- a/configs/README.md +++ b/configs/README.md @@ -15,7 +15,7 @@ A detailed description of the parameters for model configuration (i.e., initiali | forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. | | 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_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) | +| max_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 15; model not valid for soil_type = waterbody, lava, glacier, etc.) | | 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. | diff --git a/src/lgar.cxx b/src/lgar.cxx index a4b293c..5d84872 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -351,7 +351,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } else if (param_key == "max_valid_soil_types") { - state->lgar_bmi_params.num_soil_types = stoi(param_value); + state->lgar_bmi_params.num_soil_types = std::min(stoi(param_value), MAX_NUM_SOIL_TYPES); is_max_valid_soil_types_set = true; continue; } @@ -537,7 +537,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } if(!is_max_valid_soil_types_set) - state->lgar_bmi_params.num_soil_types = 12; // maximum number of soil types defaults to 12 + state->lgar_bmi_params.num_soil_types = MAX_NUM_SOIL_TYPES; // maximum number of valid soil types defaults to 15 if (verbosity.compare("high") == 0) { std::cerr<<"Maximum number of soil types: "<lgar_bmi_params.num_soil_types<<"\n"; From 0c150043cea22301329cdd1505522ad56d0bc822 Mon Sep 17 00:00:00 2001 From: AJKhattak-NOAA Date: Fri, 13 Sep 2024 16:39:57 -0400 Subject: [PATCH 31/32] fix_giuh_array_water (#34) * bug fix in computing leftover water in the giuh array. * bug fix in the giuh ordinates reader. --- include/all.hxx | 2 +- src/lgar.cxx | 46 +++++++++++--------------- tests/configs/config_lasam_synth_1.txt | 2 +- tests/configs/config_lasam_synth_2.txt | 2 +- 4 files changed, 23 insertions(+), 29 deletions(-) diff --git a/include/all.hxx b/include/all.hxx index 52750f0..a3c3fd2 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -36,7 +36,7 @@ extern string verbosity; #define use_bmi_flag FALSE // TODO set to TRUE to run in BMI environment #define MAX_NUM_SOIL_LAYERS 4 -#define MAX_NUM_SOIL_TYPES 16 +#define MAX_NUM_SOIL_TYPES 15 #define MAX_SOIL_NAME_CHARS 25 #define MAX_NUM_WETTING_FRONTS 300 diff --git a/src/lgar.cxx b/src/lgar.cxx index 5d84872..4a09fd0 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -787,12 +787,6 @@ ReadVectorData(string key) if (z1.find(delimiter) == string::npos) { double v = stod(z1); - if (v == 0.0) { - stringstream errMsg; - errMsg << "soil_z (depth of soil reservior) should be greater than zero. It it set to "<< v << " in the config file "<< "\n"; - throw runtime_error(errMsg.str()); - } - value.push_back(v); } @@ -890,29 +884,29 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run //check if the giuh queue have some water left at the end of simulaiton; needs to be included in the global mass balance // hold on; this is probably not needed as we have volrunoff in the balance; revist AJK - for(int i=1; i <= state->lgar_bmi_params.num_giuh_ordinates; i++) + for(int i=0; i <= state->lgar_bmi_params.num_giuh_ordinates; i++) volend_giuh_cm += giuh_runoff_queue_cm[i]; - double global_error_cm = volstart + volprecip - volrunoff - volAET - volon - volrech - volend + volchange_calib_cm; - - printf("\n********************************************************* \n"); - printf("-------------------- Simulation Summary ----------------- \n"); - //printf("Time (sec) = %6.10f \n", elapsed); - printf("------------------------ Mass balance ------------------- \n"); - printf("Initial water in soil = %14.10f cm\n", volstart); - printf("Total precipitation = %14.10f cm\n", volprecip); - printf("Total infiltration = %14.10f cm\n", volin); - printf("Final water in soil = %14.10f cm\n", volend); - printf("Surface ponded water = %14.10f cm\n", volon); - printf("Surface runoff = %14.10f cm\n", volrunoff); - printf("GIUH runoff = %14.10f cm\n", volrunoff_giuh); - printf("Total percolation = %14.10f cm\n", volrech); - printf("Total AET = %14.10f cm\n", volAET); - printf("Total PET = %14.10f cm\n", volPET); - printf("Total discharge (Q) = %14.10f cm\n", total_Q_cm); - printf("Vol change (calibration) = %14.10f cm\n", volchange_calib_cm); - printf("Global balance = %.6e cm\n", global_error_cm); + + printf("\n********************************************************* \n"); + printf("-------------------- Simulation Summary ----------------- \n"); + //printf("Time (sec) = %6.10f \n", elapsed); + printf("------------------------ Mass balance ------------------- \n"); + printf("Initial water in soil = %14.10f cm\n", volstart); + printf("Total precipitation = %14.10f cm\n", volprecip); + printf("Total infiltration = %14.10f cm\n", volin); + printf("Final water in soil = %14.10f cm\n", volend); + printf("Surface ponded water = %14.10f cm\n", volon); + printf("Surface runoff = %14.10f cm\n", volrunoff); + printf("GIUH runoff = %14.10f cm\n", volrunoff_giuh); + printf("GIUH water (in array) = %14.10f cm\n", volend_giuh_cm); + printf("Total percolation = %14.10f cm\n", volrech); + printf("Total AET = %14.10f cm\n", volAET); + printf("Total PET = %14.10f cm\n", volPET); + printf("Total discharge (Q) = %14.10f cm\n", total_Q_cm); + printf("Vol change (calibration) = %14.10f cm\n", volchange_calib_cm); + printf("Global balance = %.6e cm\n", global_error_cm); } diff --git a/tests/configs/config_lasam_synth_1.txt b/tests/configs/config_lasam_synth_1.txt index 761d5c0..c31f29f 100644 --- a/tests/configs/config_lasam_synth_1.txt +++ b/tests/configs/config_lasam_synth_1.txt @@ -10,4 +10,4 @@ layer_soil_type=13,14,15 max_valid_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 +giuh_ordinates=0.0 diff --git a/tests/configs/config_lasam_synth_2.txt b/tests/configs/config_lasam_synth_2.txt index e3b1da8..22ba45b 100644 --- a/tests/configs/config_lasam_synth_2.txt +++ b/tests/configs/config_lasam_synth_2.txt @@ -10,4 +10,4 @@ layer_soil_type=13,14,15 max_valid_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 +giuh_ordinates=0.0 From fe0f2245a283c884f81da5d1c4be5b14de514d4d Mon Sep 17 00:00:00 2001 From: PeterLaFollette-NOAA <59893312+peterlafollette@users.noreply.github.com> Date: Fri, 13 Sep 2024 13:40:37 -0700 Subject: [PATCH 32/32] adding HYDRUS soils params (#35) Co-authored-by: Peter La Follette --- data/vG_default_params_HYDRUS.dat | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100755 data/vG_default_params_HYDRUS.dat diff --git a/data/vG_default_params_HYDRUS.dat b/data/vG_default_params_HYDRUS.dat new file mode 100755 index 0000000..ee19052 --- /dev/null +++ b/data/vG_default_params_HYDRUS.dat @@ -0,0 +1,13 @@ +"Texture theta_r theta_e alpha (cm^-1) n Ks (cm/h)" +"Sand" 0.045 0.43 0.145 2.68 29.7 +"Loamy-sand" 0.057 0.41 0.124 2.28 14.5917 +"Sandy-loam" 0.065 0.41 0.075 1.89 4.42083 +"Silt-loam" 0.067 0.45 0.02 1.41 0.45 +"Silt" 0.034 0.46 0.016 1.37 0.25 +"Loam" 0.078 0.43 0.036 1.56 1.04 +"Sandy-clay loam" 0.1 0.39 0.059 1.48 1.31 +"Silty-clay loam" 0.089 0.43 0.01 1.23 0.07 +"Clay-loam" 0.095 0.41 0.019 1.31 0.26 +"Sandy-clay" 0.1 0.38 0.027 1.23 0.12 +"Silty-clay" 0.07 0.36 0.005 1.09 0.02 +"Clay" 0.068 0.38 0.008 1.09 0.2