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/configs/config_lasam_Bushland.txt b/configs/config_lasam_Bushland.txt index a8c273d..6ccd653 100644 --- a/configs/config_lasam_Bushland.txt +++ b/configs/config_lasam_Bushland.txt @@ -14,4 +14,4 @@ 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 -TO_enabled=false \ No newline at end of file +TO_enabled=false 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 diff --git a/include/all.hxx b/include/all.hxx index bc6992a..2d0ed0d 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -36,6 +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 25 //changed back to 25 from 15, because the file that loads soil types for Bushland relies on entries 16, 17, and 18 in the .dat file. //and generally, although we might want to restrict soil types to the 12 recognized ones and a few invalid soil types, in theory what if a user wanted to use a custom .dat file with a large number of soils? #define MAX_SOIL_NAME_CHARS 25 @@ -135,9 +136,11 @@ struct lgar_bmi_parameters double root_zone_depth_cm; // maximum depth from which roots extract water 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 bool TO_enabled = false; // if set to true, model uses multilayer TO model for recharge double mbal_tol; // if a substep's mass balance error is larger than this number, the model will abort. By default it is set to a large value (10 cm). + 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) @@ -253,6 +256,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); extern bool GW_fronts_among_surf_WFs(struct wetting_front *head); @@ -260,6 +264,7 @@ extern bool GW_fronts_among_surf_WFs(struct wetting_front *h + /*########################################*/ /* van Genuchten function prototypes */ /*########################################*/ diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index 9657492..b29339a 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(); BmiLGAR():giuh_ordinates(nullptr), giuh_runoff_queue(nullptr) { diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 8461084..237318f 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -57,6 +57,7 @@ Initialize (std::string config_file) } + /** * @brief Allocate (or reallocate) storage for soil parameters * @@ -66,6 +67,7 @@ void BmiLGAR::realloc_soil(){ 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]; } @@ -85,9 +87,11 @@ 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; @@ -97,6 +101,9 @@ Update() 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; @@ -129,8 +136,9 @@ Update() state->lgar_bmi_params.calib_params_flag = false; } + // 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) @@ -174,6 +182,7 @@ Update() bool adaptive_timestep = state->lgar_bmi_params.adaptive_timestep; double mbal_tol = state->lgar_bmi_params.mbal_tol; + // constant value used in the AET function double AET_thresh_Theta = 0.85; // scaled soil moisture (0-1) above which AET=PET (fix later!) double AET_expon = 1.0; // exponent that allows curvature of the rising portion of the Budyko curve (fix later!) @@ -191,6 +200,7 @@ Update() // 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) { subtimestep_h = state->lgar_bmi_params.minimum_timestep_h; //case where precip > 1 cm/h } @@ -204,7 +214,8 @@ Update() 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) {//and adaptive time step is engaged? + if (verbosity.compare("high") == 0) { + printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h); } @@ -411,10 +422,12 @@ Update() temp_rch += temp_excess_water; temp_runoff += temp_excess_water; + 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; @@ -467,6 +480,7 @@ Update() state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties); state->state_previous = NULL; + state->state_previous = listCopy(state->head); volin_timestep_cm += volin_subtimestep_cm; @@ -573,6 +587,8 @@ Update() volon_timestep_cm = volon_subtimestep_cm; // surface ponded water at the end of the timestep + /*----------------------------------------------------------------------*/ + // increment runoff for the subtimestep surface_runoff_subtimestep_cm = volrunoff_subtimestep_cm; surface_runoff_timestep_cm += surface_runoff_subtimestep_cm ; @@ -654,7 +670,6 @@ Update() 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 a609460..bd3577d 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -210,7 +210,6 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) fp.clear(); fp.seekg(0, fp.beg); - if (verbosity.compare("none") != 0) { std::cerr<<"------------- Initialization from config file ---------------------- \n"; } @@ -222,6 +221,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.TO_enabled = false; // setting mass balance tolerance to be large by default; this can be specified in the config file state->lgar_bmi_params.mbal_tol = 1.E1; + bool is_layer_thickness_set = false; bool is_initial_psi_set = false; @@ -426,10 +426,12 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } else if (param_key == "adaptive_timestep") { - if (param_value == "false") { + + if ((param_value == "false") || (param_value == "0")) { state->lgar_bmi_params.adaptive_timestep = false; } - else if (param_value == "true") { + else if ( (param_value == "true") || (param_value == "1")) { + state->lgar_bmi_params.adaptive_timestep = true; } else { @@ -439,6 +441,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } + else if (param_key == "mbal_tol") { state->lgar_bmi_params.mbal_tol = stod(param_value); @@ -463,6 +466,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } + else if (param_key == "timestep") { state->lgar_bmi_params.timestep_h = stod(param_value); @@ -478,6 +482,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"; @@ -1027,7 +1033,6 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run 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"); @@ -1369,8 +1374,6 @@ extern double lgar_move_wetting_fronts(bool TO_enabled, double timestep_h, doubl 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); @@ -5194,7 +5197,6 @@ extern double adjust_new_theta(int new_wf_num, double target_mass, double *cum_l next = current->next; if ((current->thetatheta) && (listLength_surface(*head)>0)){ - soil_num = soil_type[current->layer_num]; theta_e_temp = soil_properties[soil_num].theta_e; @@ -5255,6 +5257,7 @@ extern void lgar_clean_redundant_fronts(struct wetting_front** head, int *soil_t break; } + current = next; next = current->next; } @@ -5311,5 +5314,5 @@ extern double lgarto_calc_largest_mag_TO_dzdt(struct wetting_front* head){ } +#endif -#endif \ No newline at end of file diff --git a/src/linked_list.cxx b/src/linked_list.cxx index 92f94b6..1f75a58 100755 --- a/src/linked_list.cxx +++ b/src/linked_list.cxx @@ -322,7 +322,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 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..c31f29f 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 +giuh_ordinates=0.0 diff --git a/tests/configs/config_lasam_synth_2.txt b/tests/configs/config_lasam_synth_2.txt index 2559ab0..22ba45b 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 +giuh_ordinates=0.0 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