Merge branch 'master' into LGARTO

bmi/bmi.hxx

@@ -9,10 +9,10 @@
 #include <string>
 #include <vector>
 
-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;

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
+TO_enabled=false

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

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,13 +256,15 @@ 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);
 
 
 
 
 
+
 /*########################################*/
 /*   van Genuchten function prototypes    */
 /*########################################*/

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) {

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
 

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]: "<<state->lgar_bmi_params.timestep_h<<" , "
 		 <<state->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->theta<next->theta) && (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

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

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]
+ponded_depth_max=1.0[cm]

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

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

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