Compilation errors and allowing meters grid

Config_file.cfg

@@ -7,20 +7,24 @@ run_type           equilibrium              # Options: test, equilibrium, transi
 fsm_on             0                        # 1 to enable Fill-Spill-Merge for routing surface water is enabled; 0 otherwise.
 evap_mode          0                        # 1 to use a grid of potential evaporation for lakes; 0 to remove all surface water.
 infiltration_on    0                        # true is 1, false is 0. Only recommend true for high-resolution input data. 
+runoff_ratio_on    0                        # set runoff_ratio_on to 1 to supply a runoff ratio array, or 0 to assume all P-ET infiltrates in the cell
 
 # Space and time steps and extentswithin the model
+grid_type          1                        # 0 for latlon grid, 1 for xy (meters) grid
+res_meters         1                        # cell resolution in meters
 cells_per_degree   120                      # Number of cells in one degree. This example is for 5-arcminute cells. 
 southern_edge      -52                      # Southern-most latitude
 deltat             31536000                 # seconds per timestep
 total_cycles       1000                     # how many times FSM should run before completion. 
+cycles_to_save     250                      # save output every X cycles
 
 # Model parameters
 fdepth_a           200                      # e-folding depth for transmissivity calculations
 fdepth_b           150                      # e-folding depth for transmissivity calculations
 fdepth_fmin        2                        # e-folding depth coefficients 
 
 # Options relating to model convergence
-picard_iterations  5                        # How many picard iterations to perform within the groundwater matrix calculations
+# picard_iterations  5                        # How many picard iterations to perform within the groundwater matrix calculations
 maxiter            500                      # how many times GW should run before FSM runs
 
 # Starting and ending times, used to select input files

src/dephier.hpp

@@ -312,7 +312,7 @@ DepressionHierarchy<elev_t> GetDepressionHierarchy(
   // begin looking for depressions. We add all of these ocean cells to the
   // priority queue now.
   uint64_t ocean_cells = 0;
-#pragma omp parallel for default(none) shared(dem,dx,dy,label,neighbours) collapse(2) reduction(+:ocean_cells) reduction(merge:ocean_seeds)
+#pragma omp parallel for shared(dem,dx,dy,label,neighbours) collapse(2) reduction(+:ocean_cells) reduction(merge:ocean_seeds)
   for (int y = 0; y < dem.height(); y++)
     for (int x = 0; x < dem.width(); x++) {
       // Ensure the input only has OCEAN and NO_DEP labels.
@@ -366,7 +366,7 @@ DepressionHierarchy<elev_t> GetDepressionHierarchy(
   // finds and this shouldn't slow things down too much!
   int pit_cell_count = 0;
   progress.start(dem.size());
-#pragma omp parallel for default(none) shared(dem,dx,dy,label,neighbours,progress) private(ocean_seeds) collapse(2) reduction(+:pit_cell_count) reduction(merge:land_seeds)
+#pragma omp parallel for shared(dem,dx,dy,label,neighbours,progress) private(ocean_seeds) collapse(2) reduction(+:pit_cell_count) reduction(merge:land_seeds)
   for (int y = 0; y < dem.height(); y++)     // Look at all the cells
     for (int x = 0; x < dem.width(); x++) {  // Yes, all of them
       ++progress;

src/irf.cpp

@@ -263,7 +263,11 @@ void cell_size_area(Parameters& params, ArrayPack& arp) {
   const auto meters_per_degree = earth_radius * deg_to_rad;
 
   // N-S Meters per cell; distance between lines of latitude is a constant
-  params.cellsize_n_s_metres = meters_per_degree / params.cells_per_degree;
+  if (params.grid_type == 0) {
+    params.cellsize_n_s_metres = meters_per_degree / params.cells_per_degree;
+  } else {
+    params.cellsize_n_s_metres = params.res_meters;
+  }
 
   // initialise some arrays
 
@@ -272,34 +276,41 @@ void cell_size_area(Parameters& params, ArrayPack& arp) {
   // cell area (metres squared)
   arp.cell_area.resize(params.ncells_y);
 
-  // used to calculate cell latitude in radians.
-  // southern edge of the domain in degrees, plus the number of cells up from this
-  // location/the number of cells per degree, converted to radians.
-  const auto cell_position_latitude = [&](const auto cell_idx) {
-    return (cell_idx / params.cells_per_degree + params.southern_edge) * deg_to_rad;
-  };
-
-  for (int j = 0; j < params.ncells_y; j++) {
-    // latitude, in radians, at the southern edge of a cell:
-    const double latitude_radians_S = cell_position_latitude(j);
-    // latitude, in radians, at the northern edge of a cell (add a cell; equal to the southern edge of the next cell):
-    const double latitude_radians_N = cell_position_latitude(j + 1);
-
-    // distance between lines of longitude varies with latitude.
-
-    // distance at the northern edge of the cell for the given latitude:
-    const double cellsize_e_w_metres_N = params.cellsize_n_s_metres * std::cos(latitude_radians_N);
-    // distance at the southern edge of the cell for the given latitude:
-    const double cellsize_e_w_metres_S = params.cellsize_n_s_metres * std::cos(latitude_radians_S);
-
-    arp.cellsize_e_w_metres[j] = (cellsize_e_w_metres_N + cellsize_e_w_metres_S) / 2.;
-
-    // cell area computed as a trapezoid, using unchanging north-south distance,
-    // and east-west average distance.
-    arp.cell_area[j] = params.cellsize_n_s_metres * arp.cellsize_e_w_metres[j];
-
-    if (arp.cell_area[j] < 0) {
-      throw std::runtime_error("Cell with a negative area was found!");
+  if (params.grid_type == 1) {
+    for (int32_t j = 0; j < params.ncells_y; j++) {
+      arp.cellsize_e_w_metres[j] = params.res_meters;
+      arp.cell_area[j] = params.res_meters * params.res_meters;
+    }
+  } else {
+    // used to calculate cell latitude in radians.
+    // southern edge of the domain in degrees, plus the number of cells up from this
+    // location/the number of cells per degree, converted to radians.
+    const auto cell_position_latitude = [&](const auto cell_idx) {
+      return (cell_idx / params.cells_per_degree + params.southern_edge) * deg_to_rad;
+    };
+
+    for (int j = 0; j < params.ncells_y; j++) {
+      // latitude, in radians, at the southern edge of a cell:
+      const double latitude_radians_S = cell_position_latitude(j);
+      // latitude, in radians, at the northern edge of a cell (add a cell; equal to the southern edge of the next cell):
+      const double latitude_radians_N = cell_position_latitude(j + 1);
+
+      // distance between lines of longitude varies with latitude.
+
+      // distance at the northern edge of the cell for the given latitude:
+      const double cellsize_e_w_metres_N = params.cellsize_n_s_metres * std::cos(latitude_radians_N);
+      // distance at the southern edge of the cell for the given latitude:
+      const double cellsize_e_w_metres_S = params.cellsize_n_s_metres * std::cos(latitude_radians_S);
+
+      arp.cellsize_e_w_metres[j] = (cellsize_e_w_metres_N + cellsize_e_w_metres_S) / 2.;
+
+      // cell area computed as a trapezoid, using unchanging north-south distance,
+      // and east-west average distance.
+      arp.cell_area[j] = params.cellsize_n_s_metres * arp.cellsize_e_w_metres[j];
+
+      if (arp.cell_area[j] < 0) {
+        throw std::runtime_error("Cell with a negative area was found!");
+      }
     }
   }
 }

src/parameters.cpp

@@ -22,6 +22,10 @@ Parameters::Parameters(const std::string& config_file) {
       continue;
     }
 
+    if (line.find("#")==0) {
+      continue;
+    }
+
     std::stringstream ss(line);
     std::string key;
     ss >> key;
@@ -44,6 +48,8 @@ Parameters::Parameters(const std::string& config_file) {
       ss >> fdepth_fmin;
     } else if (key == "fsm_on") {
       ss >> fsm_on;
+    } else if (key == "grid_type") {
+      ss >> grid_type;
     } else if (key == "infiltration_on") {
       ss >> infiltration_on;
     } else if (key == "maxiter") {
@@ -52,6 +58,8 @@ Parameters::Parameters(const std::string& config_file) {
       ss >> outfile_prefix;
     } else if (key == "region") {
       ss >> region;
+    } else if (key == "res_meters") {
+      ss >> res_meters;
     } else if (key == "run_type") {
       ss >> run_type;
     } else if (key == "runoff_ratio_on") {
@@ -98,12 +106,21 @@ void Parameters::check() const {
     }
   };
 
-  check_positive("cells_per_degree", cells_per_degree);
+  check_binary(
+	       grid_type, "set grid_type to 0 for a latlon grid, or 1 for an xy (meters) grid.");
+
+  if (grid_type == 0) {
+    check_positive("cells_per_degree", cells_per_degree);
+    if (std::isnan(southern_edge) || southern_edge < -90 || southern_edge > 90) {
+      throw std::runtime_error("please enter a value between -90 and 90 degrees for the southern_edge!");
+    }
+  } else {
+    check_positive("res_meters", res_meters);
+  }
+
   check_positive("cycles_to_save", cycles_to_save);
   check_positive("deltat", deltat);
-  if (std::isnan(southern_edge) || southern_edge < -90 || southern_edge > 90) {
-    throw std::runtime_error("please enter a value between -90 and 90 degrees for the southern_edge!");
-  }
+
   check_binary(
       evap_mode,
       "set evap_mode to 0 to remove all surface water, or 1 to use a grid of potential evaporation for lakes.");
@@ -148,7 +165,9 @@ std::string Parameters::get_path(const std::string& layer_name) const {
 }
 
 void Parameters::print() const {
-  std::cout << "c cells_per_degree       = " << cells_per_degree << std::endl;
+  if (grid_type == 0) {
+    std::cout << "c cells_per_degree       = " << cells_per_degree << std::endl;
+  } 
   std::cout << "c cycles_to_save         = " << cycles_to_save << std::endl;
   std::cout << "c deltat                 = " << deltat << std::endl;
   std::cout << "c evap_mode              = " << evap_mode << std::endl;
@@ -160,9 +179,14 @@ void Parameters::print() const {
   std::cout << "c maxiter                = " << maxiter << std::endl;
   std::cout << "c outfile_prefix         = " << outfile_prefix << std::endl;
   std::cout << "c region                 = " << region << std::endl;
+  if (grid_type == 1) {
+    std::cout << "c res_meters             = " << res_meters << std::endl;
+  }
   std::cout << "c run_type               = " << run_type << std::endl;
   std::cout << "c runoff_ratio_on        = " << runoff_ratio_on << std::endl;
-  std::cout << "c southern_edge          = " << southern_edge << std::endl;
+  if (grid_type == 0) {
+    std::cout << "c southern_edge          = " << southern_edge << std::endl;
+  }
   std::cout << "c supplied_wt            = " << supplied_wt << std::endl;
   std::cout << "c surfdatadir            = " << surfdatadir << std::endl;
   std::cout << "c textfilename           = " << textfilename << std::endl;

src/parameters.hpp

@@ -26,6 +26,7 @@ struct Parameters {
   std::string time_end       = UNINIT_STR;
 
   double cells_per_degree = -1;
+  double res_meters = -1;
 
   double UNDEF = -1.0e7;
 
@@ -34,6 +35,7 @@ struct Parameters {
   int32_t evap_mode       = -1;
   int32_t fsm_on          = -1;
   int32_t runoff_ratio_on = -1;
+  int32_t grid_type = -1;
 
   double deltat          = std::numeric_limits<double>::signaling_NaN();
   double fdepth_a        = -1.;

src/run_dephier.cpp

@@ -43,43 +43,54 @@ int main(int argc, char** argv) {
   const auto meters_per_degree = earth_radius * deg_to_rad;
 
   // distance between lines of latitude is a constant.
-  params.cellsize_n_s_metres = meters_per_degree / params.cells_per_degree;
+  if (params.grid_type == 0) {
+    params.cellsize_n_s_metres = meters_per_degree / params.cells_per_degree;
+  } else {
+    params.cellsize_n_s_metres = params.res_meters;
+  }
 
   // initialise some arrays
   // size of a cell in the east-west direction at the centre of the cell (metres)
   arp.cellsize_e_w_metres.resize(params.ncells_y);
   // cell area (metres squared)
   arp.cell_area.resize(params.ncells_y);
 
-  // used to calculate cell latitude in radians.
-  // southern edge of the domain in degrees, plus the number of cells up from this
-  // location/the number of cells per degree, converted to radians.
-  const auto cell_position_latitude = [&](const auto cell_idx) {
-    return (cell_idx / params.cells_per_degree + params.southern_edge) * deg_to_rad;
-  };
-  for (int32_t j = 0; j < params.ncells_y; j++) {
-    // southern edge of the domain in degrees, plus the number of cells up
-    // from this location/the number of cells per degree, converted to radians.
-    // latitude at the southern edge of a cell (subtract half a cell):
-    const double latitude_radians_S = cell_position_latitude(j);
-    // latitude at the northern edge of a cell (add half a cell):
-    const double latitude_radians_N = cell_position_latitude(j + 1);
-
-    // distance at the northern edge of the cell for the given latitude:
-    const double cellsize_e_w_metres_N = params.cellsize_n_s_metres * std::cos(latitude_radians_N);
-    // distance at the southern edge of the cell for the given latitude:
-    const double cellsize_e_w_metres_S = params.cellsize_n_s_metres * std::cos(latitude_radians_S);
-
-    // distance between lines of longitude varies with latitude.
-    // This is the distance at the centre of a cell for a given latitude:
-    arp.cellsize_e_w_metres[j] = (cellsize_e_w_metres_N + cellsize_e_w_metres_S) / 2.;
-
-    // cell area computed as a trapezoid, using unchanging north-south distance,
-    // and east-west distances at the northern and southern edges of the cell:
-    arp.cell_area[j] = params.cellsize_n_s_metres * arp.cellsize_e_w_metres[j];
-
-    if (arp.cell_area[j] < 0) {
-      throw std::runtime_error("Cell with a negative area was found!");
+  if (params.grid_type == 1) {
+    for (int32_t j = 0; j < params.ncells_y; j++) {
+      arp.cellsize_e_w_metres[j] = params.res_meters;
+      arp.cell_area[j] = params.res_meters * params.res_meters;
+    }
+  } else {
+    // used to calculate cell latitude in radians.
+    // southern edge of the domain in degrees, plus the number of cells up from this
+    // location/the number of cells per degree, converted to radians.
+    const auto cell_position_latitude = [&](const auto cell_idx) {
+      return (cell_idx / params.cells_per_degree + params.southern_edge) * deg_to_rad;
+    };
+    for (int32_t j = 0; j < params.ncells_y; j++) {
+      // southern edge of the domain in degrees, plus the number of cells up
+      // from this location/the number of cells per degree, converted to radians.
+      // latitude at the southern edge of a cell (subtract half a cell):
+      const double latitude_radians_S = cell_position_latitude(j);
+      // latitude at the northern edge of a cell (add half a cell):
+      const double latitude_radians_N = cell_position_latitude(j + 1);
+
+      // distance at the northern edge of the cell for the given latitude:
+      const double cellsize_e_w_metres_N = params.cellsize_n_s_metres * std::cos(latitude_radians_N);
+      // distance at the southern edge of the cell for the given latitude:
+      const double cellsize_e_w_metres_S = params.cellsize_n_s_metres * std::cos(latitude_radians_S);
+
+      // distance between lines of longitude varies with latitude.
+      // This is the distance at the centre of a cell for a given latitude:
+      arp.cellsize_e_w_metres[j] = (cellsize_e_w_metres_N + cellsize_e_w_metres_S) / 2.;
+
+      // cell area computed as a trapezoid, using unchanging north-south distance,
+      // and east-west distances at the northern and southern edges of the cell:
+      arp.cell_area[j] = params.cellsize_n_s_metres * arp.cellsize_e_w_metres[j];
+
+      if (arp.cell_area[j] < 0) {
+        throw std::runtime_error("Cell with a negative area was found!");
+      }
     }
   }