diff --git a/Config_file.cfg b/Config_file.cfg index 87fa189d..295d9b19 100644 --- a/Config_file.cfg +++ b/Config_file.cfg @@ -7,12 +7,16 @@ 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 @@ -20,7 +24,7 @@ fdepth_b 150 # e-folding depth for transmissivity 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 diff --git a/src/dephier.hpp b/src/dephier.hpp index 0bdc0f34..63856553 100644 --- a/src/dephier.hpp +++ b/src/dephier.hpp @@ -312,7 +312,7 @@ DepressionHierarchy 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 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; diff --git a/src/irf.cpp b/src/irf.cpp index ddadb55d..911bf572 100644 --- a/src/irf.cpp +++ b/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!"); + } } } } diff --git a/src/parameters.cpp b/src/parameters.cpp index 3eee4693..f668e7f9 100644 --- a/src/parameters.cpp +++ b/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; diff --git a/src/parameters.hpp b/src/parameters.hpp index a1635424..249365f1 100644 --- a/src/parameters.hpp +++ b/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::signaling_NaN(); double fdepth_a = -1.; diff --git a/src/run_dephier.cpp b/src/run_dephier.cpp index 9071519b..38eac1c9 100644 --- a/src/run_dephier.cpp +++ b/src/run_dephier.cpp @@ -43,7 +43,11 @@ 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) @@ -51,35 +55,42 @@ int main(int argc, char** argv) { // 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!"); + } } }