Merge pull request #711 from rarutter/baseline_climate_config

This moves the calendar year range selection for EQ averaged climate to the
config file for easier modification without recompiling. It includes basic
checks for the specified years being out of bounds for the input file or the
end year being prior to the start year.

Author: rarutter

config/config.js

@@ -121,7 +121,9 @@
   },
 
   "model_settings": {
-    "dynamic_lai": 1                    // from model (1) or from input (0)
+    "dynamic_lai": 1,                   // from model (1) or from input (0)
+    "baseline_start": 1901,  //start year for baseline EQ climate
+    "baseline_end": 1931     //end year for baseline EQ climate
 //    //"dynamic_climate": 0,
 //    //"varied_co2": 0,
 //    //"fire_severity_as_input": 0,    // fire sev. as input or ??

include/Climate.h

@@ -10,10 +10,16 @@ class Climate {
   Climate();
   Climate(const std::string& fname, const std::string& co2fname, int y, int x);
 
-  // misc climate variables
-  // This value will change during a run when switching from
+  // Misc. climate variables
+  // The following two values determine the span of historic climate
+  //  years used to produce a baseline average climate for EQ
+  //  TODO: SP repeated climate is not currently affected by these.
+  int baseline_start;
+  int baseline_end;
+  // These values will change during a run when switching from
   //  historic to projected climate data
   int tseries_start_year;
+  int tseries_end_year;
 
   // driving variables
   std::vector<float> co2;
@@ -75,6 +81,8 @@ class Climate {
   void load_proj_climate(const std::string&, int, int);
   void load_proj_co2(const std::string&);
 
+  void prep_avg_climate();
+
 private:
 
   void load_from_file(const std::string& fname, int y, int x);

include/ModelData.h

@@ -86,6 +86,10 @@ class ModelData {
   bool nc_tr;
   bool nc_sc;
   int output_interval; //How many years to store for output
+  //The following two config values are temporarily stored in
+  // ModelData, to be transferred to Climate.
+  int baseline_start;//Start year for baseline EQ climate
+  int baseline_end;//End year for baseline EQ climate
 
   // Maps holding data about variables to be output at specific timesteps
   // C++11 would allow the use of unordered_maps, which have a faster

include/TEMUtilityFunctions.h

@@ -145,6 +145,7 @@ namespace temutil {
                                     const int y, const int x);
 
   int get_timeseries_start_year(const std::string &filename);
+  int get_timeseries_end_year(const std::string &filename);
 
   // draft - reads all timesteps co2 data
   std::vector<float> get_timeseries(const std::string &filename,

src/Climate.cpp

@@ -405,6 +405,7 @@ void Climate::load_from_file(const std::string& fname, int y, int x) {
     nirr = temutil::get_timeseries<float>(fname, "nirr", y, x);
 
     tseries_start_year = temutil::get_timeseries_start_year(fname);
+    tseries_end_year = temutil::get_timeseries_end_year(fname);
 
   }//End critical(load_climate)
 
@@ -454,14 +455,15 @@ void Climate::load_from_file(const std::string& fname, int y, int x) {
   BOOST_LOG_SEV(glg, debug) << "par = [" << temutil::vec2csv(par) << "]";
 
   // Create a simplified historic climate for EQ by averaging input data
-  // over a year range specified here. The year choices should be exposed
-  // in the config file eventually. TODO
-  if(fname.find("historic") != std::string::npos){
-    avgX_tair = avg_over(tair, 1901, 1931);
-    avgX_prec = avg_over(prec, 1901, 1931);
-    avgX_nirr = avg_over(nirr, 1901, 1931);
-    avgX_vapo = avg_over(vapo, 1901, 1931);
-  }
+  // over a year range specified here.
+//Commenting out but leaving temporarily for easier comparison with
+// the currently in-progress daily data ingestion branch. 20240429
+//  if(fname.find("historic") != std::string::npos){
+//    avgX_tair = avg_over(tair, baseline_start, baseline_end);
+//    avgX_prec = avg_over(prec, 1901, 1931);
+//    avgX_nirr = avg_over(nirr, 1901, 1931);
+//    avgX_vapo = avg_over(vapo, 1901, 1931);
+//  }
 
   // Do we need simplified 'avgX_' values for par, and cld??
   // ===> YES: the derived variables should probably be based off the avgX
@@ -474,6 +476,41 @@ void Climate::load_from_file(const std::string& fname, int y, int x) {
 
 }
 
+/** Prepares simplified average climate, intended for EQ stage
+*
+* Assumes that averaged climate is produced shortly after
+*  the tseries values have been loaded from the intended
+*  input file.
+*/
+void Climate::prep_avg_climate(){
+  BOOST_LOG_SEV(glg, debug) << "Climate baseline from config: "
+                            << baseline_start << ":" << baseline_end;
+
+  if(   baseline_start > tseries_end_year
+     || baseline_start < tseries_start_year
+     || baseline_end > tseries_end_year
+     || baseline_end < tseries_start_year){
+
+    BOOST_LOG_SEV(glg, fatal) << "Baseline year value exceeds range"
+             << " of input file.\n" << "Acceptable range: "
+             << tseries_start_year << ":" << tseries_end_year << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  if(baseline_end < baseline_start){
+    BOOST_LOG_SEV(glg, fatal) << "Baseline end year " << baseline_end
+                              << " is less than baseline start year "
+                              << baseline_start;
+    exit(EXIT_FAILURE);
+  }
+
+  avgX_tair = avg_over(tair, baseline_start, baseline_end);
+  avgX_prec = avg_over(prec, baseline_start, baseline_end);
+  avgX_nirr = avg_over(nirr, baseline_start, baseline_end);
+  avgX_vapo = avg_over(vapo, baseline_start, baseline_end);
+}
+
+
 /** This loads data from a projected climate data file, overwriting any old climate data*/
 void Climate::load_proj_climate(const std::string& fname, int y, int x){
   BOOST_LOG_SEV(glg, note) << "Climate, loading projected data";
@@ -518,9 +555,7 @@ std::vector<float> Climate::avg_over(const std::vector<float> & var, const int s
   }
 
   BOOST_LOG_SEV(glg, debug) << "result = [" << temutil::vec2csv(result) << "]";
-
   return result;
-
 }
 
 

src/Cohort.cpp

@@ -63,13 +63,17 @@ Cohort::Cohort(int y, int x, ModelData* modeldatapointer):
 
   // might need to set the cd* and the ed* ???
 
-  BOOST_LOG_SEV(glg, debug) << "Setup the NEW STYLE CLIMATE OBJECT ...";
-  // FIX: Historic? Projected?? how to handle both??
-  // Maybe:
-  //this->hist_climate = Climate(modeldatapointer->hist_climate, y, x);
-  //this->proj_climate = Climate(modeldatapointer->proj_climate, y, x);
+  BOOST_LOG_SEV(glg, debug) << "Construct Climate object ...";
+
+  //On construction we assume historic climate, which will be
+  // overwritten by projected climate later when necessary.
   this->climate = Climate(modeldatapointer->hist_climate_file, modeldatapointer->co2_file, y, x);
-  
+
+  this->climate.baseline_start = modeldatapointer->baseline_start;
+  this->climate.baseline_end = modeldatapointer->baseline_end;
+  //Prepare averaged input set for EQ stage
+  this->climate.prep_avg_climate();
+
   // Build a mineral info object
   MineralInfo mineral_info = MineralInfo(modeldatapointer->soil_texture_file, y, x);
 
@@ -853,6 +857,7 @@ void Cohort::updateMonthly_DIMveg(const int & currmind, const bool & dynamic_lai
   for (int ip=0; ip<NUM_PFT; ip++) {
     if (cd.m_veg.vegcov[ip]>0.) {
       cd.m_veg.vegage[ip] = cd.yrsdist;
+
       if (cd.m_veg.vegage[ip]<=0) {
         cd.m_vegd.foliagemx[ip] = 0.;
       }
@@ -1518,4 +1523,3 @@ void Cohort::set_restartdata_from_state() {
     }
   }
 }
-

src/ModelData.cpp

@@ -112,6 +112,8 @@ ModelData::ModelData(Json::Value controldata):force_cmt(-1) {
   caldata_tree_loc  = controldata["calibration-IO"]["caldata_tree_loc"].asString();
 
   dynamic_LAI       = controldata["model_settings"]["dynamic_lai"].asInt(); // checked in Cohort::updateMonthly_DIMVeg
+  baseline_start = controldata["model_settings"]["baseline_start"].asInt();
+  baseline_end   = controldata["model_settings"]["baseline_end"].asInt();
 
   // Unused (11/23/2015)
   //changeclimate = controldata["model_settings"]["dynamic_climate"].asInt();

src/TEMUtilityFunctions.cpp

@@ -659,6 +659,50 @@ namespace temutil {
     return start_year;
   }
 
+  /** Return the calendar end year of a timeseries netcdf file.
+  * Assumes that input files have complete years
+  * Assumes that time units are "days since..."
+  */
+  int get_timeseries_end_year(const std::string& fname){
+
+    int start_year = get_timeseries_start_year(fname);
+
+    int ncid;
+    temutil::nc( nc_open(fname.c_str(), NC_NOWRITE, &ncid) );
+
+    //Information about the time *dimension*
+    int timeD;
+    size_t timeD_len;
+    temutil::nc( nc_inq_dimid(ncid, "time", &timeD) );
+    temutil::nc( nc_inq_dimlen(ncid, timeD, &timeD_len) );
+
+    //Information about the time *variable*
+    int timeV;
+    temutil::nc( nc_inq_varid(ncid, "time", &timeV) );
+
+    size_t start[3];
+    start[0] = timeD_len-1; //Last entry only
+    start[1] = 0;
+    start[2] = 0;
+
+    size_t count[3];
+    count[0] = 1;
+    count[1] = 1;
+    count[2] = 1;
+
+    int time_value;
+    temutil::nc( nc_get_vara_int(ncid, timeV, start, count, &time_value) );
+
+    //Round up, because the time value will be for the beginning
+    // of the last time step and so will not divide evenly by 365.
+    //Casting one of the values to a double to force use of the
+    // proper operator/
+    int year_count = ceil(double(time_value) / DINY);
+
+    int end_year = start_year + year_count;
+    return end_year;
+  }
+
   /** rough draft - look up lon/lat in nc file from y,x coordinates. 
       Assumes that the file has some coordinate dimensions...
   */