Commit 13e0f43b authored by Ethan Coon's avatar Ethan Coon
Browse files

adds compile-time variant of Matrix for use in Legion

parent c835cf5f
......@@ -5,22 +5,24 @@
namespace ELM {
void CanopyHydrologyKern1(double dtime,
const double& forc_rain, const double& forc_snow, const double& irrig_rate,
const int& ltype, const int& ctype,
const bool& urbpoi, const bool& do_capsnow,
const double& elai, const double& esai,
const double& dewmx, const int& frac_veg_nosno,
double& h2ocan,
int& n_irrig_steps_left,
double& qflx_prec_intr,
double& qflx_irrig,
double& qflx_prec_grnd,
double& qflx_snwcp_liq,
double& qflx_snwcp_ice,
double& qflx_snow_grnd_patch,
double& qflx_rain_grnd) {
canopyhydrologykern1_(&dtime, &forc_rain, &forc_snow, &irrig_rate,
void CanopyHydrology_Interception(double dtime,
const double& forc_rain,
const double& forc_snow,
const double& irrig_rate,
const int& ltype, const int& ctype,
const bool& urbpoi, const bool& do_capsnow,
const double& elai, const double& esai,
const double& dewmx, const int& frac_veg_nosno,
double& h2ocan,
int& n_irrig_steps_left,
double& qflx_prec_intr,
double& qflx_irrig,
double& qflx_prec_grnd,
double& qflx_snwcp_liq,
double& qflx_snwcp_ice,
double& qflx_snow_grnd_patch,
double& qflx_rain_grnd) {
canopyhydrology_interception_(&dtime, &forc_rain, &forc_snow, &irrig_rate,
&ltype, &ctype, &urbpoi, &do_capsnow, &elai, &esai,
&dewmx, &frac_veg_nosno, &h2ocan, &n_irrig_steps_left,
&qflx_prec_intr, &qflx_irrig, &qflx_prec_grnd,
......
......@@ -3,7 +3,7 @@
extern "C" {
void canopyhydrologykern1_(const double* dtime,
void canopyhydrology_interception_(const double* dtime,
const double* forc_rain,
const double* forc_snow,
const double* irrig_rate,
......
subroutine CanopyHydrologyKern1( dtime, &
subroutine CanopyHydrology_Interception( dtime, &
forc_rain, forc_snow, forc_irrig, &
ltype, ctype, urbpoi, do_capsnow, &
elai, esai, dewmx, frac_veg_nosno, &
......@@ -41,10 +41,8 @@ subroutine CanopyHydrologyKern1( dtime, &
real(r8) :: fracsnow
real(r8) :: fracrain
! Canopy interception and precipitation onto ground surface
! Add precipitation to leaf water
if (ltype==istsoil .or. ltype==istwet .or. urbpoi .or. &
ltype==istcrop) then
......@@ -80,7 +78,7 @@ subroutine CanopyHydrologyKern1( dtime, &
! Intercepted precipitation [mm/s]
qflx_prec_intr = (forc_snow + forc_rain) * fpi
!print*, forc_rain, forc_snow, elai, esai, qflx_prec_intr, fpi
! print*, forc_rain, forc_snow, elai, esai, qflx_prec_intr, fpi
! Water storage of intercepted precipitation and dew
......@@ -162,5 +160,5 @@ subroutine CanopyHydrologyKern1( dtime, &
return
end subroutine CanopyHydrologyKern1
end subroutine CanopyHydrology_Interception
subroutine CanopyHydrologyKern2( dtime, oldfflag, newnode, &
ctype, qflx_floodc, qflx_snow_h2osfc, snow_depth, snl, &
swe_old, h2osoi_liq, h2osoi_ice, dz, z, zi, &
t_soisno, frac_iceold, &
do_capsnow, frac_h2osfc, qflx_snow_grnd_col, frac_sno, int_snow, forc_t, &
h2osno,qflx_snow_melt, n_melt, frac_sno_eff, t_grnd, &
qflx_floodg, &
ltype, urbpoi) bind(C)
subroutine CanopyHydrology_SnowWater( dtime, &
qflx_floodg, &
ctype, ltype, urbpoi, do_capsnow, oldfflag, &
forc_t, t_grnd, qflx_snow_grnd_col, qflx_snow_melt, n_melt, frac_h2osfc, & ! forcing
snow_depth, h2osno, swe_old, &
h2osoi_liq, h2osoi_ice, t_soisno, frac_iceold, & ! state
snl, dz, z, zi, newnode, & ! snow mesh for initialization
qflx_floodc, qflx_snow_h2osfc, &
int_snow, frac_sno_eff, frac_sno)
use shr_kind_mod, only : &
r8 => shr_kind_r8, &
i4 => shr_kind_i4, &
......@@ -26,28 +27,29 @@ subroutine CanopyHydrologyKern2( dtime, oldfflag, newnode, &
real(r8), parameter :: tfrz=273.15
real(r8) :: zlnd = 0.01_r8
real(r8), intent(in) :: dtime
integer, intent(in) :: oldfflag, ctype , ltype
logical, intent(in) :: do_capsnow, urbpoi
real(r8), intent(in) :: qflx_floodg
real(r8), intent(in) :: frac_h2osfc, qflx_snow_grnd_col, forc_t , qflx_snow_melt, n_melt
real(r8), intent(in) :: t_grnd
integer, intent(inout) :: snl
real(r8), intent(inout) :: snow_depth , h2osno
real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
real(r8), intent(inout), dimension(-nlevsno+1:0) :: dz
integer, intent(out) :: newnode
real(r8), intent(out) :: qflx_floodc
real(r8), intent(in) :: qflx_floodg
real(r8), intent(out) :: qflx_snow_h2osfc
real(r8), intent(out), dimension(-nlevsno+1:0) :: swe_old
real(r8), intent(inout) :: snow_depth , h2osno
real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
real(r8), intent(inout), dimension(-nlevsno+1:0) :: dz
real(r8), intent(out), dimension(-nlevsno+1:0) :: z, zi
real(r8), intent(out), dimension(-nlevsno+1:0) :: t_soisno, frac_iceold
logical, intent(in) :: do_capsnow, urbpoi
real(r8), intent(in) :: frac_h2osfc, qflx_snow_grnd_col, forc_t , qflx_snow_melt, n_melt
real(r8), intent(out) :: int_snow, frac_sno_eff, frac_sno
real(r8), intent(in) :: t_grnd
! local variables
real(r8) :: temp_intsnow, temp_snow_depth, z_avg, fmelt, dz_snowf, snowmelt
real(r8) :: newsnow, bifall, fsnow_new, accum_factor, fsno_new, smr
real(r8) :: newsnow, bifall, accum_factor, fsno_new, smr
integer :: j
! apply gridcell flood water flux to non-lake columns
......@@ -236,5 +238,5 @@ subroutine CanopyHydrologyKern2( dtime, oldfflag, newnode, &
dz(snl+1) = dz(snl+1)+dz_snowf*dtime
end if
end subroutine CanopyHydrologyKern2
end subroutine CanopyHydrology_SnowWater
......@@ -17,8 +17,10 @@ OBJS = \
column_varcon.o \
clm_varpar.o \
clm_varctl.o \
CanopyHydrologyKern1.o \
CanopyHydrologyKern2.o
CanopyHydrology_Interception.o \
CanopyHydrology_SnowWater.o \
CanopyHydrology_FracWet.o \
CanopyHydrology_FracH2OSfc.o
all: $(OBJS)
$(AR) cr libelm.a $(OBJS)
......
......@@ -15,15 +15,30 @@
#include "CanopyHydrology.hh"
namespace ELM {
namespace Utils {
static const int n_months = 12;
static const int n_pfts = 17;
static const int n_max_times = 31 * 24 * 2; // max days per month times hours per
// day * half hour timestep
static const int n_grid_cells = 24;
using MatrixState = MatrixStatic<n_grid_cells, n_pfts>;
using MatrixForc = MatrixStatic<n_max_times,n_grid_cells>;
} // namespace
} // namespace
int main(int argc, char ** argv)
{
// dimensions
const int n_months = 12;
const int n_pfts = 17;
const int n_grid_cells = 24;
const int n_max_times = 31 * 24 * 2; // max days per month times hours per
// day * half hour timestep
using ELM::Utils::n_months;
using ELM::Utils::n_pfts;
using ELM::Utils::n_grid_cells;
using ELM::Utils::n_max_times;
// fixed magic parameters for now
const int ctype = 1;
const int ltype = 1;
......@@ -36,30 +51,30 @@ int main(int argc, char ** argv)
const double dtime = 1800.0;
// phenology state
ELM::Utils::Matrix<> elai(n_grid_cells, n_pfts);
ELM::Utils::Matrix<> esai(n_grid_cells, n_pfts);
ELM::Utils::MatrixState elai;
ELM::Utils::MatrixState esai;
ELM::Utils::read_phenology("../links/surfacedataWBW.nc", n_months, n_pfts, 0, elai, esai);
ELM::Utils::read_phenology("../links/surfacedataBRW.nc", n_months, n_pfts, n_months, elai, esai);
// forcing state
ELM::Utils::Matrix<> forc_rain(n_max_times, n_grid_cells);
ELM::Utils::Matrix<> forc_snow(n_max_times, n_grid_cells);
ELM::Utils::Matrix<> forc_air_temp(n_max_times, n_grid_cells);
ELM::Utils::MatrixForc forc_rain;
ELM::Utils::MatrixForc forc_snow;
ELM::Utils::MatrixForc forc_air_temp;
const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 0, n_grid_cells, forc_rain, forc_snow, forc_air_temp);
ELM::Utils::Matrix<> forc_irrig(n_max_times, n_grid_cells); forc_irrig = 0.;
ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
// output state by the grid cell
auto qflx_prec_intr = std::vector<double>(n_grid_cells, 0.);
auto qflx_irrig = std::vector<double>(n_grid_cells, 0.);
auto qflx_prec_grnd = std::vector<double>(n_grid_cells, 0.);
auto qflx_snwcp_liq = std::vector<double>(n_grid_cells, 0.);
auto qflx_snwcp_ice = std::vector<double>(n_grid_cells, 0.);
auto qflx_snow_grnd_patch = std::vector<double>(n_grid_cells, 0.);
auto qflx_rain_grnd = std::vector<double>(n_grid_cells, 0.);
auto qflx_prec_intr = std::array<double,n_grid_cells>();
auto qflx_irrig = std::array<double,n_grid_cells>();
auto qflx_prec_grnd = std::array<double,n_grid_cells>();
auto qflx_snwcp_liq = std::array<double,n_grid_cells>();
auto qflx_snwcp_ice = std::array<double,n_grid_cells>();
auto qflx_snow_grnd_patch = std::array<double,n_grid_cells>();
auto qflx_rain_grnd = std::array<double,n_grid_cells>();
// output state by the pft
auto h2o_can = ELM::Utils::Matrix<>(n_grid_cells, n_pfts); h2o_can = 0.;
auto h2o_can = ELM::Utils::MatrixState(); h2o_can = 0.;
std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water" << std::endl;
auto min_max = std::minmax_element(h2o_can.begin(), h2o_can.end());
......@@ -78,7 +93,7 @@ int main(int argc, char ** argv)
// NOTE: this currently punts on what to do with the qflx variables!
// Surely they should be either accumulated or stored on PFTs as well.
// --etc
ELM::CanopyHydrologyKern1(dtime,
ELM::CanopyHydrology_Interception(dtime,
forc_rain(t,g), forc_snow(t,g), forc_irrig(t,g),
ltype, ctype, urbpoi, do_capsnow,
elai(g,p), esai(g,p), dewmx, frac_veg_nosno,
......
......@@ -14,6 +14,19 @@
#include "readers.hh"
#include "CanopyHydrology.hh"
namespace ELM {
namespace Utils {
static const int n_months = 12;
static const int n_pfts = 17;
using MatrixState = MatrixStatic<n_months, n_pfts>;
static const int n_max_times = 31 * 24 * 2; // max days per month times hours per
// day * half hour timestep
using MatrixForc = MatrixStatic<n_max_times,1>;
} // namespace
} // namespace
int main(int argc, char ** argv)
{
......@@ -35,14 +48,14 @@ int main(int argc, char ** argv)
const double dtime = 1800.0;
// phenology state
ELM::Utils::Matrix<> elai(n_months, n_pfts);
ELM::Utils::Matrix<> esai(n_months, n_pfts);
ELM::Utils::MatrixState elai;
ELM::Utils::MatrixState esai;
ELM::Utils::read_phenology("../links/surfacedataWBW.nc", n_months, n_pfts, 0, elai, esai);
// forcing state
ELM::Utils::Matrix<> forc_rain(n_max_times, 1);
ELM::Utils::Matrix<> forc_snow(n_max_times, 1);
ELM::Utils::Matrix<> forc_air_temp(n_max_times, 1);
ELM::Utils::MatrixForc forc_rain;
ELM::Utils::MatrixForc forc_snow;
ELM::Utils::MatrixForc forc_air_temp;
const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 6, 1, forc_rain, forc_snow, forc_air_temp);
double h2ocan = 0.0;
......@@ -57,10 +70,10 @@ int main(int argc, char ** argv)
std::cout << "Timestep, forc_rain, h2ocan, qflx_prec_grnd, qflx_prec_intr" << std::endl;
for(size_t itime = 0; itime < n_times; itime += 1) {
// note this call puts all precip as rain for testing
double total_precip = forc_rain(itime, 0) + forc_snow(itime, 0);
ELM::CanopyHydrologyKern1(dtime, total_precip, 0., 0.,
double total_precip = forc_rain[itime][0] + forc_snow[itime][0];
ELM::CanopyHydrology_Interception(dtime, total_precip, 0., 0.,
ltype, ctype, urbpoi, do_capsnow,
elai(5,7), esai(5,7), dewmx, frac_veg_nosno,
elai[5][7], esai[5][7], dewmx, frac_veg_nosno,
h2ocan, n_irrig_steps_left,
qflx_prec_intr, qflx_irrig, qflx_prec_grnd,
qflx_snwcp_liq, qflx_snwcp_ice,
......
......@@ -21,7 +21,7 @@ default: all
all: links library $(TESTS)
test: $(EXEC_TESTS)
@echo "done"
python ../compare_to_gold.py $(TESTS)
CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
......
......@@ -134,13 +134,13 @@ int read_forcing(const std::string& fname,
// allocate the precip to rain or snow
for (int lcv_t=0; lcv_t!=ntimes_l; ++lcv_t) {
if (data_temp[lcv_t] < 273.15) {
snow(lcv_t, lcv_gc) = data_precip[lcv_t];
rain(lcv_t, lcv_gc) = 0.;
temp(lcv_t, lcv_gc) = data_temp[lcv_t];
snow[lcv_t][lcv_gc] = data_precip[lcv_t];
rain[lcv_t][lcv_gc] = 0.;
temp[lcv_t][lcv_gc] = data_temp[lcv_t];
} else {
snow(lcv_t, lcv_gc) = 0.;
rain(lcv_t, lcv_gc) = data_precip[lcv_t];
temp(lcv_t, lcv_gc) = data_temp[lcv_t];
snow[lcv_t][lcv_gc] = 0.;
rain[lcv_t][lcv_gc] = data_precip[lcv_t];
temp[lcv_t][lcv_gc] = data_temp[lcv_t];
}
}
}
......
......@@ -7,8 +7,36 @@ namespace ELM {
namespace Utils {
enum struct Ordering { C, FORTRAN };
template<size_t ROW, size_t COL, typename T=double>
class MatrixStatic {
public:
MatrixStatic() {}
T& operator()(size_t i, size_t j) { return d_[i][j]; }
const T& operator()(size_t i, size_t j) const { return d_[i][j]; }
std::array<T,COL>& operator[](size_t i) { return d_[i]; }
const std::array<T,COL>& operator[](size_t i) const { return d_[i]; }
void operator=(T t) {
for (size_t i=0; i!=ROW; ++i) {
for (size_t j=0; j!=COL; ++j) {
d_[i][j] = t;
}
}
}
double const * begin() const { return &d_[0][0]; }
double const * end() const { return &d_[ROW-1][COL-1] +1; }
private:
std::array<std::array<T,COL>,ROW> d_;
};
enum struct Ordering { C, FORTRAN };
template <typename T=double, Ordering O=Ordering::C>
class Matrix {
......
......@@ -77,7 +77,7 @@ program CanopyHydrology_kern1_multiple
esai = surfdata(g)%sai(p)
h2ocan = h2ocan_pft(p,g)
call CanopyHydrologyKern1( dtime, &
call CanopyHydrology_Interception( dtime, &
forc_rain, forc_snow, irrig_rate, &
ltype, ctype, urbpoi, do_capsnow, &
elai, esai, dewmx, frac_veg_nosno, &
......
......@@ -78,7 +78,7 @@ program CanopyHydrology_kern1_single
n_irrig_steps_left = 0
urbpoi=.false.
do_capsnow = .false.
call CanopyHydrologyKern1( dtime, &
call CanopyHydrology_Interception( dtime, &
forc_rain, forc_snow, irrig_rate, &
ltype, ctype, urbpoi, do_capsnow, &
elai, esai, dewmx, frac_veg_nosno, &
......
......@@ -6,10 +6,12 @@ FCFLAGS += -I$(SRCDIR) -I$(SRCDIR)$(ELM_UTILS_DIR) -I$(NETCDF_ROOT)/include
TESTS = test_CanopyHydrology_kern1_single \
test_CanopyHydrology_kern1_multiple
test_CanopyHydrology_kern1_multiple \
test_CanopyHydrology_module
EXEC_TESTS = CanopyHydrology_kern1_single \
CanopyHydrology_kern1_multiple
CanopyHydrology_kern1_multiple \
CanopyHydrology_module
.PHONY: links library test
......@@ -27,6 +29,9 @@ CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
CanopyHydrology_kern1_multiple: test_CanopyHydrology_kern1_multiple
./test_CanopyHydrology_kern1_multiple &> test_CanopyHydrology_kern1_multiple.stdout
CanopyHydrology_module: test_CanopyHydrology_module
./test_CanopyHydrology_module &> test_CanopyHydrology_module.stdout
test_%: %.o
$(LINKER) $(FCFLAGS) -L$(OBJECT)fortran -lelm -L$(NETCDF_ROOT)/lib -lnetcdff $< -o $@
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment