Merge pull request #5 from NOAA-OWP/ajk/calib_params

calibratable_parameters

README.md

@@ -52,7 +52,6 @@ _________________________________________________________________
 | soil_storage_model | string | conceptual or layered or topmodel | - | - | if `conceptual`, conceptual models are used for computing the soil moisture profile (e.g., CFE). If `layered`, layered-based soil moisture models are used (e.g., LGAR). If `topmodel`, topmodel's variables are used
 | soil_moisture_profile_option | string | constant or linear | - | - | `constant` for layered-constant profile. `linear`  for linearly interpolated values between two consecutive layers. Needed if `soil_storage_model = layered`.
 | soil_depth_layers | double (1D array) | - | - | - | Absolute depth of soil layers. Needed if `soil_storage_model = layered`.
-| max_num_wetting_fronts | int | - | - | - | Maximum number of wetting fronts. Default is set to 30. Needed if `soil_storage_model = layered`.
 | soil_moisture_fraction_depth | double | (0, domain_depth] | m | - | *user specified depth for the soil moisture fraction (default is 40 cm)
 | water_table_based_method | string | flux-based or deficit-based | - | - | Needed if `soil_storage_model = topmodel`. `flux-based` uses an iterative scheme, and `deficit-based` uses catchment deficit to compute soil moisture profile
 

include/bmi_soil_moisture_profile.hxx

@@ -19,19 +19,24 @@ namespace coupler {
 class BmiSoilMoistureProfile : public bmixx::Bmi {
 public:
   BmiSoilMoistureProfile() {
-    this->input_var_names[0] = "soil_storage";
-    this->input_var_names[1] = "soil_storage_change";
-    this->input_var_names[2] = "soil_moisture_wetting_fronts";
-    this->input_var_names[3] = "soil_depth_wetting_fronts";
-    this->input_var_names[4] = "num_wetting_fronts";
-    this->input_var_names[5] = "Qb_topmodel";        // baseflow in the topmodel
-    this->input_var_names[6] = "Qv_topmodel";        // recharge rate to the saturated zone to the un saturated zone
+    this->input_var_names[0]  = "soil_storage";
+    this->input_var_names[1]  = "soil_storage_change";
+    this->input_var_names[2]  = "num_wetting_fronts";
+    this->input_var_names[3]  = "soil_moisture_wetting_fronts";
+    this->input_var_names[4]  = "soil_depth_wetting_fronts";
+    this->input_var_names[5]  = "Qb_topmodel";        // baseflow in the topmodel
+    this->input_var_names[6]  = "Qv_topmodel";        // recharge rate to the saturated zone to the un saturated zone
                                                      // in the topmodel
-    this->input_var_names[7] = "global_deficit";     // global soil deficit in the topmodel
+    this->input_var_names[7]  = "global_deficit";    // global soil deficit in the topmodel
 
     this->output_var_names[0] = "soil_moisture_profile";  // entire profile of the soil column (1D array)
     this->output_var_names[1] = "soil_water_table";       // depth of the water table from the surface in meters
     this->output_var_names[2] = "soil_moisture_fraction"; // fraction of soil moisture, top 0.4 m (or user specified depth)
+
+    // add calibratable parameters
+    this->calib_var_names[0]  = "smcmax";
+    this->calib_var_names[1]  = "b";
+    this->calib_var_names[2]  = "satpsi";
   };
 
   void Initialize(std::string config_file);
@@ -87,13 +92,17 @@ public:
   void GetGridFaceEdges(const int grid, int *face_edges);
   void GetGridFaceNodes(const int grid, int *face_nodes);
   void GetGridNodesPerFace(const int grid, int *nodes_per_face);
+  void ResetSize (std::string name);
+
 private:
   soil_moisture_profile::soil_profile_parameters* state;
   static const int input_var_name_count  = 8;
   static const int output_var_name_count = 3;
+  static const int calib_var_name_count  = 3;
 
   std::string input_var_names[input_var_name_count];
   std::string output_var_names[output_var_name_count];
+  std::string calib_var_names[calib_var_name_count];
   std::string verbosity;
 };
 

include/soil_moisture_profile.hxx

@@ -18,7 +18,7 @@
   @param soil_moisture_wetting_fronts  [-] : soil moisture content of wetting fronts (bmi input to layered the model)
   @param soil_depth_wetting_fronts     [m] : absolute depth of the wetting fronts (bmi input to layered the model)
   @param smcmax                    [-] : maximum soil moisture content (porosity)
-  @param bb                        [-] : pore size distribution, beta exponent in Clapp-Hornberger (1978) function
+  @param b                         [-] : pore size distribution, beta exponent in Clapp-Hornberger (1978) function
   @param satpsi                    [m] : saturated capillary head (saturated moisture potential)
   @param ncells                    [-] : number of cells of the discretized soil column
   @param nlayers                   [-] : number of soil moisture layers, typically different than the ncells
@@ -75,7 +75,7 @@ namespace soil_moisture_profile {
     double *soil_moisture_profile;
 
     double *smcmax;
-    double bb;
+    double b;
     double satpsi;
     int    ncells;
     double soil_depth;

src/bmi_soil_moisture_profile.cxx

@@ -62,10 +62,14 @@ GetVarGrid(std::string name)
     return 1;
   else if (name.compare("Qb_topmodel") == 0 || name.compare("Qv_topmodel") == 0 || name.compare("global_deficit") == 0)
     return 1;
+  else if (name.compare("b") == 0 || name.compare("satpsi") == 0)
+    return 1;
   else if (name.compare("soil_moisture_profile") == 0) // array of doubles (conceptual model)
     return 2;
   else if (name.compare("soil_moisture_wetting_fronts") == 0 || name.compare("soil_depth_wetting_fronts") == 0) // array of doubles (layered model)
-    return 3; 
+    return 3;
+  else if (name.compare("smcmax") == 0) // fixed number of layers for calibratable params
+    return 4;
   else
     return -1;
 }
@@ -78,7 +82,7 @@ GetVarType(std::string name)
 
   if (var_grid == 0)
     return "int";
-  else if (var_grid == 1 || var_grid == 2 || var_grid == 3)
+  else if (var_grid == 1 || var_grid == 2 || var_grid == 3 || var_grid == 4)
     return "double";
   else
     return "";
@@ -152,6 +156,9 @@ GetGridShape(const int grid, int *shape)
   else if (grid == 3) {
     shape[0] = this->state->shape[1];
   }
+  else if (grid == 4) {
+    shape[0] = this->state->num_layers;
+  }
 }
 
 
@@ -192,6 +199,8 @@ GetGridSize(const int grid)
     return this->state->shape[0];
   else if (grid == 3)
     return this->state->shape[1];
+  else if (grid == 4)
+    return this->state->num_layers;
   else
     return -1;
 }
@@ -322,6 +331,12 @@ GetValuePtr (std::string name)
     return (void*)(&this->state->Qv_topmodel);
   else if (name.compare("global_deficit") == 0)
     return (void*)(&this->state->global_deficit);
+  else if (name.compare("smcmax") == 0)
+    return (void*)this->state->smcmax;
+  else if (name.compare("b") == 0)
+    return (void*)(&this->state->b);
+  else if (name.compare("satpsi") == 0)
+    return (void*)(&this->state->satpsi);
   else {
     std::stringstream errMsg;
     errMsg << "variable "<< name << " does not exist";
@@ -345,27 +360,45 @@ GetValueAtIndices (std::string name, void *dest, int *inds, int len)
     char *ptr;
 
     itemsize = this->GetVarItemsize(name);
-
     for (i=0, ptr=(char *)dest; i<len; i++, ptr+=itemsize) {
       offset = inds[i] * itemsize;
       memcpy(ptr, (char *)src + offset, itemsize);
     }
   }
 }
 
-
+void BmiSoilMoistureProfile::
+ResetSize (std::string name)
+{
+// reset the size of wetting fronts array to the number of wetting fronts at the timestep
+  if (name.compare("soil_moisture_wetting_fronts") == 0) {
+    assert (this->state->num_wetting_fronts > 0);
+    state->soil_moisture_wetting_fronts = new double[this->state->num_wetting_fronts]();
+  }
+  else if (name.compare("soil_depth_wetting_fronts") == 0) {
+    assert (this->state->num_wetting_fronts > 0);
+    state->soil_depth_wetting_fronts = new double[this->state->num_wetting_fronts]();
+  }
+}
+
 void BmiSoilMoistureProfile::
 SetValue (std::string name, void *src)
 {
   void * dest = NULL;
+  ResetSize(name);
+
   dest = this->GetValuePtr(name);
-
+  
   if (dest) {
     int nbytes = 0;
     nbytes = this->GetVarNbytes(name);
     memcpy(dest, src, nbytes);
+
+    if (name.compare("num_wetting_fronts") == 0)
+      this->state->shape[1] = this->state->num_wetting_fronts;
+
   }
-
+  
 }
 
 
@@ -374,19 +407,24 @@ SetValueAtIndices (std::string name, int * inds, int len, void *src)
 {
   void * dest = NULL;
 
+  ResetSize(name);
+
   dest = this->GetValuePtr(name);
-
+  
   if (dest) {
     int i;
     int itemsize = 0;
     int offset;
     char *ptr;
 
     itemsize = this->GetVarItemsize(name);
-
+    
     for (i=0, ptr=(char *)src; i<len; i++, ptr+=itemsize) {
       offset = inds[i] * itemsize;
       memcpy((char *)dest + offset, ptr, itemsize);
+
+      if (name.compare("num_wetting_fronts") == 0)
+	this->state->shape[1] = this->state->num_wetting_fronts;
     }
   }
 }
@@ -395,7 +433,7 @@ SetValueAtIndices (std::string name, int * inds, int len, void *src)
 std::string BmiSoilMoistureProfile::
 GetComponentName()
 {
-  return "SoilMoistureProfile BMI";
+  return "SoilMoistureProfiles BMI";
 }
 
 

src/soil_moisture_profile.cxx

@@ -27,12 +27,12 @@ SoilMoistureProfile(string config_file, struct soil_profile_parameters* paramete
   if (parameters->soil_storage_model == Conceptual || parameters->soil_storage_model == Topmodel)
     parameters->shape[1] = 1;
   else if (parameters->soil_storage_model == Layered)
-    parameters->shape[1] = parameters->max_num_wetting_fronts;
+    parameters->shape[1] = parameters->num_wetting_fronts;
 
   parameters->soil_moisture_profile = new double[parameters->ncells];
 
-  parameters->soil_moisture_wetting_fronts = new double[parameters->shape[1]];
-  parameters->soil_depth_wetting_fronts    = new double[parameters->shape[1]];
+  parameters->soil_moisture_wetting_fronts = new double[parameters->shape[1]]();
+  parameters->soil_depth_wetting_fronts    = new double[parameters->shape[1]]();
 
   // For water_table_based_method
   parameters->cat_area     = 1.0;        // catchment area used in the topmodel (normalized)
@@ -43,16 +43,16 @@ SoilMoistureProfile(string config_file, struct soil_profile_parameters* paramete
   parameters->origin[1]    = 0.;
   parameters->soil_storage = 0.0;
   parameters->init_profile = true;
-  parameters->num_wetting_fronts = 1;
   parameters->soil_depth_NWM = 2.0;
+  parameters->num_wetting_fronts = parameters->shape[1];
   parameters->soil_storage_change_per_timestep = 0.0;
 }
 
 /*
   Read and initialize values from configuration file
   @input soil_z   (1D)        :   soil discretization; array of depths from the surface [m]
   @input layers_z  (1D)       : depth of each layer from the surface [m]
-  @input bb  (double)         : pore size distribution [-], beta exponent on Clapp-Hornberger (1978)
+  @input b  (double)          : pore size distribution [-], beta exponent on Clapp-Hornberger (1978)
   @input satpsi  (double)     : saturated capillary head (saturated moisture potential) [m]
   @input ncells  (int)        : number of cells of the discretized soil column
   @input ncells_layered       : number of layers (wetting fronts) for the layered model
@@ -84,11 +84,10 @@ InitFromConfigFile(string config_file, struct soil_profile_parameters* parameter
   bool is_soil_z_set                    = false;
   bool is_soil_depth_layers_set         = false;
   bool is_smcmax_set                    = false;
-  bool is_bb_set                        = false;
+  bool is_b_set                        = false;
   bool is_satpsi_set                    = false;
   bool is_soil_storage_model_set        = false;
   bool is_soil_storage_model_depth_set  = false;
-  bool is_max_num_wetting_fronts_set    = false;
   bool is_water_table_depth_set         = false;
   bool is_water_table_based_method_set  = false;
   bool is_soil_moisture_fraction_depth_set        = false;
@@ -166,9 +165,9 @@ InitFromConfigFile(string config_file, struct soil_profile_parameters* parameter
       continue;
     }
     else if (param_key == "soil_params.b") {
-      parameters->bb = stod(param_value);
-      assert (parameters->bb > 0);
-      is_bb_set = true;
+      parameters->b = stod(param_value);
+      assert (parameters->b > 0);
+      is_b_set = true;
       continue;
     }
     else if (param_key == "soil_params.satpsi") {
@@ -201,12 +200,6 @@ InitFromConfigFile(string config_file, struct soil_profile_parameters* parameter
       is_soil_storage_model_depth_set = true;
       continue;
     }
-    else if (param_key == "max_num_wetting_fronts") {
-      parameters->max_num_wetting_fronts =  stod(param_value);
-      assert (parameters->max_num_wetting_fronts > 0);
-      is_max_num_wetting_fronts_set = true;
-      continue;
-    }
     else if (param_key == "water_table_depth") {
       parameters->water_table_depth = stod(param_value);
       is_water_table_depth_set = true;
@@ -261,9 +254,9 @@ InitFromConfigFile(string config_file, struct soil_profile_parameters* parameter
     throw runtime_error(errMsg.str());
   }
 
-  if (!is_bb_set) {
+  if (!is_b_set) {
     stringstream errMsg;
-    errMsg << "bb (Clapp-Hornberger's parameter) not set in the config file "<< config_file << "\n";
+    errMsg << "b (Clapp-Hornberger's parameter) not set in the config file "<< config_file << "\n";
     throw runtime_error(errMsg.str());
   }
 
@@ -296,16 +289,12 @@ InitFromConfigFile(string config_file, struct soil_profile_parameters* parameter
       errMsg << "soil_moisture_profile_option_set key is not set in the config file "<< config_file << ", options = constant or linear \n";
       throw runtime_error(errMsg.str());
     }
-
-    if (!is_max_num_wetting_fronts_set) {
-      parameters->max_num_wetting_fronts = 30;
-    }
 
     if (!is_water_table_depth_set) {
       parameters->water_table_depth = 6.0;
     }
 
-    assert (parameters->max_num_wetting_fronts > 0);
+    assert (parameters->num_wetting_fronts > 0);
     assert (parameters->water_table_depth >= 0);
   }
 
@@ -339,6 +328,7 @@ SoilMoistureProfileUpdate(struct soil_profile_parameters* parameters)
   double soil_moisture_fraction_depth = parameters->soil_moisture_fraction_depth;
   parameters->soil_moisture_fraction = 0.0; // reset to 0 at each timestep
 
+
   if (parameters->soil_storage_model == Conceptual) {
     SoilMoistureProfileFromConceptualReservoir(parameters);
   }
@@ -403,7 +393,7 @@ SoilMoistureProfileUpdate(struct soil_profile_parameters* parameters)
   Computes 1D soil moisture profile for conceptual reservoir using Newton-Raphson iterative method
   For detailed decription of the model implemented here, please see README.md on the github repo
   local_variables:
-  @param lam  [-] : 1/bb (bb: pore size distribution)
+  @param lam  [-] : 1/b (b: pore size distribution)
   @param satpsi_cm [cm] : saturated moisture potential
   @param soil_storage_model_depth [m] : depth of the soil reservoir model (e.g., CFE)
   @param zb  [cm] : bottom of the computational domain
@@ -429,7 +419,7 @@ SoilMoistureProfileFromConceptualReservoir(struct soil_profile_parameters* param
   double zb          = 0.0;  // bottom of the computational domain
   double z0          = 0.0;  // bottom of the fictitious domain (to track fictitious water table location)
   double zi          = 0.01; // initial guess for the water table location, use Newton-Raphson to find new zi
-  double lam         = 1.0/parameters->bb;
+  double lam         = 1.0/parameters->b;
   double beta        = 1.0 - lam;
   double alpha       = pow(satpsi_cm,lam)/beta; // a constant term obtained in the integration of the soil moisture function
   double tolerance   = 1.0E-6;
@@ -575,7 +565,7 @@ SoilMoistureProfileFromConceptualReservoir(struct soil_profile_parameters* param
 		     value between consecutive layers.
   - Note: the water table location is the thickness of the water table plus saturated capillary head (satpsi)
   - local_variables:
-  @param lam                [-]  : 1/bb (bb: pore size distribution)
+  @param lam                [-]  : 1/b (b: pore size distribution)
   @param soil_depth         [cm] : depth of the soil column
   @param last_layer_depth   [cm] : depth of the last layer from the surface
   @param tolerance          [-]  : tolerance to find location of the water table
@@ -597,7 +587,7 @@ SoilMoistureProfileFromLayeredReservoir(struct soil_profile_parameters* paramete
 
   parameters->last_layer_depth = parameters->soil_depth_wetting_fronts[num_wf-1];
 
-  double lam = 1.0/parameters->bb; // pore distribution index
+  double lam = 1.0/parameters->b; // pore distribution index
 
   vector<double> z_layers_n(1,0.0);
 
@@ -718,7 +708,7 @@ FindWaterTableLayeredReservoir(struct soil_profile_parameters* parameters)
   int num_wf            = parameters->num_wetting_fronts; //number of wetting fronts
   int num_layers        = parameters->num_layers;
   double tolerance      = 1.0e-3;
-  double lam = 1.0/parameters->bb; // pore distribution index
+  double lam = 1.0/parameters->b; // pore distribution index
   bool is_water_table_found = false;
 
   // find and update water table location if it is within the model domain
@@ -818,7 +808,7 @@ ExtendedProfileBelowDomainDepth(struct soil_profile_parameters* parameters) {
   Computes 1D soil moisture profile for models (e.g. Topmodel) using the water table depth
   For detailed decription of the model implemented here, please see README.md on the github repo
   local_variables:
-  @param lam  [-] : 1/bb (bb: pore size distribution)
+  @param lam  [-] : 1/b (b: pore size distribution)
   @param satpsi_cm [cm] : saturated moisture potential
   @param soil_moisture_profile [-] : OUTPUT (soil moisture content vertical profile [-])
   @param dt                    [h] : topmodel's timestep
@@ -834,7 +824,7 @@ SoilMoistureProfileFromWaterTableDepth(struct soil_profile_parameters* parameter
 {
   // converting variables to cm for numerical reasons only
   double satpsi_cm   = parameters->satpsi * 100.;
-  double lam         = 1.0/parameters->bb;
+  double lam         = 1.0/parameters->b;
   double model_depth = 600;
   double dt = 1.0;
   double to_cm = 100;