Commit 148c48d8 authored by Pillai, Himanshu's avatar Pillai, Himanshu
Browse files

Module implemetation in pure C++, Legion and Kokkos

parent 661a6afb
#ifndef ELM_CANOPY_HYDROLOGY_FORT_HH_
#define ELM_CANOPY_HYDROLOGY_FORT_HH_
#include "CanopyHydrology_private.h"
namespace ELM {
inline 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,
&qflx_snwcp_liq, &qflx_snwcp_ice, & qflx_snow_grnd_patch, &qflx_rain_grnd);
}
inline void CanopyHydrology_FracWet(const int& frac_veg_nosno,
const double& h2ocan,
const double& elai,
const double& esai,
const double& dewmx,
double& fwet,
double& fdry) {
canopyhydrology_fracwet_(&frac_veg_nosno, &h2ocan, &elai, &esai, &dewmx, &fwet, &fdry);
}
template<typename Array_d>
void CanopyHydrology_SnowWater(const double& dtime,
const double& qflx_floodg,
const int& ltype,
const int& ctype,
const bool& urbpoi,
const bool& do_capsnow,
const int& oldfflag,
const double& forc_t,
const double& t_grnd,
const double& qflx_snow_grnd_col,
const double& qflx_snow_melt,
const double& n_melt,
const double& frac_h2osfc,
double& snow_depth,
double& h2osno,
double& int_snow,
Array_d& swe_old,
Array_d& h2osoi_liq,
Array_d& h2osoi_ice,
Array_d& t_soisno,
Array_d& frac_iceold,
int& snl,
Array_d& dz,
Array_d& z,
Array_d& zi,
int& newnode,
double& qflx_floodc,
double& qflx_snow_h2osfc,
double& frac_sno_eff,
double& frac_sno) {
canopyhydrology_snowwater_(&dtime, &qflx_floodg, &ltype, &ctype, &urbpoi,
&do_capsnow, &oldfflag, &forc_t, &t_grnd,
&qflx_snow_grnd_col, &qflx_snow_melt, &n_melt,
&frac_h2osfc, &snow_depth, &h2osno, &int_snow,
&swe_old[0], &h2osoi_liq[0], &h2osoi_ice[0], &t_soisno[0],
&frac_iceold[0], &snl, &dz[0], &z[0], &zi[0], &newnode,
&qflx_floodc, &qflx_snow_h2osfc, &frac_sno_eff, &frac_sno);
}
inline void CanopyHydrology_FracH2OSfc(const double& dtime,
const double& min_h2osfc,
const int& ltype,
const double& micro_sigma,
const double& h2osno,
double& h2osfc,
double& h2osoi_liq,
double& frac_sno,
double& frac_sno_eff,
double& qflx_h2osfc2topsoi,
double& frac_h2osfc,
bool no_update=false) {
canopyhydrology_frach2osfc_(&dtime, &min_h2osfc, &ltype, &micro_sigma, &h2osno,
&h2osfc, &h2osoi_liq, &frac_sno, &frac_sno_eff, &qflx_h2osfc2topsoi,
&frac_h2osfc, &no_update);
}
} // namespace
#endif
#ifndef ELM_CANOPY_HYDROLOGY_FORT_PRIVATE_HH_
#define ELM_CANOPY_HYDROLOGY_FORT_PRIVATE_HH_
extern "C" {
void canopyhydrology_interception_(const 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,
const double* h2ocan,
const int* n_irrig_steps_left,
const double* qflx_prec_intr,
const double* qflx_irrig,
const double* qflx_prec_grnd,
const double* qflx_snwcp_liq,
const double* qflx_snwcp_ice,
const double* qflx_snow_grnd_patch,
const double* qflx_rain_grnd);
void canopyhydrology_fracwet_(const int* frac_veg_nosno,
const double* h2ocan,
const double* elai,
const double* esai,
const double* dewmx,
const double* fwet,
const double* fdry);
void canopyhydrology_snowwater_(const double* dtime,
const double* qflx_floodg,
const int* ltype,
const int* ctype,
const bool* urbpoi,
const bool* do_capsnow,
const int* oldfflag,
const double* forc_t,
const double* t_grnd,
const double* qflx_snow_grnd_col,
const double* qflx_snow_melt,
const double* n_melt,
const double* frac_h2osfc,
const double* snow_depth,
const double* h2osno,
const double* int_snow,
const double* swe_old,
const double* h2osoi_liq,
const double* h2osoi_ice,
const double* t_soisno,
const double* frac_iceold,
const int* snl,
const double* dz,
const double* z,
const double* zi,
const int* newnode,
const double* qflx_floodc,
const double* qflx_snow_h2osfc,
const double* frac_sno_eff,
const double* frac_sno);
void canopyhydrology_frach2osfc_(const double* dtime,
const double* min_h2osfc,
const int* ltype,
const double* micro_sigma,
const double* h2osno,
const double* h2osfc,
const double* h2osoi_liq,
const double* frac_sno,
const double* frac_sno_eff,
const double* qflx_h2osfc2topsoi,
const double* frac_h2osfc,
const bool* no_update);
}
#endif
#include <algorithm>
#include <stdio.h>
#include <math.h>
#include <cmath>
#include <iostream>
#include <string>
#include "landunit_varcon.h"
......
#include <stdio.h>
#include <math.h>
#include <cmath>
#include "landunit_varcon.h"
#include "column_varcon.h"
#include "clm_varpar.h"
......
......@@ -3,20 +3,20 @@
# ------------------------------------------------------------------
# Available from: git clone https://code.ornl.gov/isp/acme_inputdata
#ELM_DATA_LOCATION ?= /home/7hp/Downloads/acme_inputdata
ELM_DATA_LOCATION ?= /Users/uec/codes/elm/kernels/repos/acme_inputdata
ELM_DATA_LOCATION ?= /home/7hp/Downloads/acme_inputdata
#ELM_DATA_LOCATION ?= /Users/uec/codes/elm/kernels/repos/acme_inputdata
# assumes you have a working NETCDF. This can easily be the version built by
# ATS as a part of an ATS installation.
NETCDF_ROOT ?= ${AMANZI_TPLS_DIR}
NETCDF_ROOT ?= /usr/local
# assumes you have a working gfortran
FC = gfortran
FC_FLAGS += -fimplicit-none -free -g3 -fbounds-check -fbacktrace
FC_LIB_ROOT = /usr/local/Cellar/gcc/8.2.0/lib/gcc/8
FC_LIB_ROOT = /usr/lib/gcc/x86_64-linux-gnu/8 #/usr/local/Cellar/gcc/8.2.0/lib/gcc/8
# assumes you have a working c++ compiler
CC = clang++
CC = g++
CPP_FLAGS += -g -Wall -Wshadow -std=c++11
STD_LIB_ROOT = /usr
......@@ -34,9 +34,9 @@ CUDA_INCDIRS=-I/usr/local/cuda-10.0/include
# linking flags
CC_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT)/lib -lnetcdf -L$(STD_LIB_ROOT)/lib -lstdc++ -L$(FC_LIB_ROOT) -lgfortran
FC_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT)/lib -lnetcdff
CUDA_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT)/lib -lnetcdf -L$(CUDA_LIBDIRS) -lstdc++ -lopenblas -lpthread -lcudart -lcublas
CC_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT) -lnetcdf -L$(STD_LIB_ROOT)/lib -lstdc++ -L$(FC_LIB_ROOT) -lgfortran
FC_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT) -lnetcdff
CUDA_LD_FLAGS += -L$(SRCDIR) -lelm -L$(NETCDF_ROOT) -lnetcdf -L$(CUDA_LIBDIRS) -lstdc++ -lopenblas -lpthread -lcudart -lcublas
#
# rules
......
......@@ -5,9 +5,9 @@ Error codes:
G | Error in the gold file -- bad testname?
D | Error in my file -- test did not complete?
E | Error in my file size -- test did not complete.
* | Unknown error
F | Failing test
. | Passing test
N/A | Unknown error
FAIL | Failing test
PASS | Passing test
----------------------------------------------------------
"""
......@@ -49,9 +49,9 @@ def run_comparison(testname, full_message=False):
return 'E'
if close:
return '.'
return 'PASS'
else:
return 'F'
return 'FAIL'
if __name__ == "__main__":
......@@ -69,5 +69,5 @@ if __name__ == "__main__":
results.append(run_comparison(test, args.verbose))
except Exception as err:
print err
results.append('*')
results.append('N/A')
print ''.join(results)
......@@ -9,7 +9,7 @@
#include <iostream>
#include <iomanip>
#include <numeric>
#include <algorithm>
#include "utils.hh"
#include "readers.hh"
......@@ -90,19 +90,19 @@ int main(int argc, char ** argv)
auto dz = ELM::Utils::MatrixStateSoilColumn(0.);
// state variables that require ICs and evolve (in/out)
auto h2ocan = ELM::Utils::MatrixStatePFT(0.);
auto h2ocan = ELM::Utils::MatrixStatePFT(); h2ocan = 0.;
auto swe_old = ELM::Utils::MatrixStateSoilColumn(0.);
auto h2osoi_liq = ELM::Utils::MatrixStateSoilColumn(0.);
auto h2osoi_ice = ELM::Utils::MatrixStateSoilColumn(0.);
auto t_soisno = ELM::Utils::MatrixStateSoilColumn(0.);
auto frac_iceold = ELM::Utils::MatrixStateSoilColumn(0.);
auto t_grnd = ELM::Utils::VectorColumn(0.);
auto h2osno = ELM::Utils::VectorColumn(0.);
auto h2osno = ELM::Utils::VectorColumn(0.); h2osno = 0.;
auto snow_depth = ELM::Utils::VectorColumn(0.);
auto snl = ELM::Utils::VectorColumnInt(0.); // note this tracks the snow_depth
auto snow_level = ELM::Utils::VectorColumnInt(0.); // note this tracks the snow_depth
auto h2osfc = ELM::Utils::VectorColumn(0.);
auto frac_h2osfc = ELM::Utils::VectorColumn(0.);
auto frac_h2osfc = ELM::Utils::VectorColumn(0.); frac_h2osfc = 0.;
// output fluxes by pft
......@@ -116,7 +116,7 @@ int main(int argc, char ** argv)
// FIXME: I have no clue what this is... it is inout on WaterSnow. For now I
// am guessing the data structure. Ask Scott. --etc
auto int_snow = ELM::Utils::VectorColumn(0.);
auto integrated_snow = ELM::Utils::VectorColumn(0.);
// output fluxes, state by the column
auto qflx_snow_grnd_col = ELM::Utils::VectorColumn();
......@@ -126,14 +126,22 @@ int main(int argc, char ** argv)
auto frac_sno_eff = ELM::Utils::VectorColumn();
auto frac_sno = ELM::Utils::VectorColumn();
{
std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water\t Total Snow\t Min Snow\t Max Snow\t Avg Frac Sfc\t Min Frac Sfc\t Max Frac Sfc" << std::endl;
auto min_max_water = std::minmax_element(h2ocan.begin(), h2ocan.end());
auto sum_water = std::accumulate(h2ocan.begin(), h2ocan.end(), 0.);
auto min_max_snow = std::minmax_element(h2osno.begin(), h2osno.end());
auto sum_snow = std::accumulate(h2osno.begin(), h2osno.end(), 0.);
auto min_max_frac_sfc = std::minmax_element(frac_h2osfc.begin(), frac_h2osfc.end());
auto avg_frac_sfc = std::accumulate(frac_h2osfc.begin(), frac_h2osfc.end(), 0.) / (frac_h2osfc.end() - frac_h2osfc.begin());
// std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water" << std::endl;
// auto min_max = std::minmax_element(h2ocan.begin(), h2ocan.end());
// std::cout << std::setprecision(16)
// << 0 << "\t" << std::accumulate(h2ocan.begin(), h2ocan.end(), 0.)
// << "\t" << *min_max.first
// << "\t" << *min_max.second << std::endl;
std::cout << std::setprecision(16)
<< 0 << "\t" << sum_water << "\t" << *min_max_water.first << "\t" << *min_max_water.second
<< "\t" << sum_snow << "\t" << *min_max_snow.first << "\t" << *min_max_snow.second
<< "\t" << avg_frac_sfc << "\t" << *min_max_frac_sfc.first << "\t" << *min_max_frac_sfc.second << std::endl;
}
// main loop
// -- the timestep loop cannot/should not be parallelized
......@@ -186,9 +194,9 @@ int main(int argc, char ** argv)
ltype, ctype, urbpoi, do_capsnow, oldfflag,
forc_air_temp(t,g), t_grnd(g),
qflx_snow_grnd_col[g], qflx_snow_melt, n_melt, frac_h2osfc[g],
snow_depth[g], h2osno[g], int_snow[g], swe_old[g],
snow_depth[g], h2osno[g], integrated_snow[g], swe_old[g],
h2osoi_liq[g], h2osoi_ice[g], t_soisno[g], frac_iceold[g],
snl[g], dz[g], z[g], zi[g], newnode,
snow_level[g], dz[g], z[g], zi[g], newnode,
qflx_floodc[g], qflx_snow_h2osfc[g], frac_sno_eff[g], frac_sno[g]);
// Calculate Fraction of Water to the Surface?
......@@ -202,13 +210,19 @@ int main(int argc, char ** argv)
} // end grid cell loop
auto min_max_water = std::minmax_element(h2ocan.begin(), h2ocan.end());
auto sum_water = std::accumulate(h2ocan.begin(), h2ocan.end(), 0.);
auto min_max_snow = std::minmax_element(h2osno.begin(), h2osno.end());
auto sum_snow = std::accumulate(h2osno.begin(), h2osno.end(), 0.);
// auto min_max = std::minmax_element(h2ocan.begin(), h2ocan.end());
// std::cout << std::setprecision(16)
// << t+1 << "\t" << std::accumulate(h2ocan.begin(), h2ocan.end(), 0.)
// << "\t" << *min_max.first
// << "\t" << *min_max.second << std::endl;
auto min_max_frac_sfc = std::minmax_element(frac_h2osfc.begin(), frac_h2osfc.end());
auto avg_frac_sfc = std::accumulate(frac_h2osfc.begin(), frac_h2osfc.end(), 0.) / (frac_h2osfc.end() - frac_h2osfc.begin());
std::cout << std::setprecision(16)
<< 0 << "\t" << sum_water << "\t" << *min_max_water.first << "\t" << *min_max_water.second
<< "\t" << sum_snow << "\t" << *min_max_snow.first << "\t" << *min_max_snow.second
<< "\t" << avg_frac_sfc << "\t" << *min_max_frac_sfc.first << "\t" << *min_max_frac_sfc.second << std::endl;
} // end timestep loop
return 0;
}
......@@ -4,7 +4,7 @@ SRCDIR = $(OBJECT)$(KERNEL_FOLDER)
include $(OBJECT)config/Makefile.config
INC_FLAGS = -I$(AMANZI_TPLS_DIR)/include -I../../src/cpp/ -I../test_c_c
INC_FLAGS = -I$(NETCDF_ROOT)/include -I$(SRCDIR)
TESTS = test_CanopyHydrology_kern1_single \
test_CanopyHydrology_kern1_multiple \
......
......@@ -9,7 +9,7 @@
#include <iostream>
#include <iomanip>
#include <numeric>
#include <algorithm>
#include "utils.hh"
#include "readers.hh"
......
......@@ -4,7 +4,8 @@ OBJECT = ../../src/
include $(OBJECT)config/Makefile.config
INC_FLAGS = -I$(AMANZI_TPLS_DIR)/include -I$(CSRCDIR) -I../tests_c
INC_FLAGS = -I$(NETCDF_ROOT)/include -I$(CSRCDIR) -I../tests_c
TESTS = test_CanopyHydrology_kern1_single \
......
......@@ -204,8 +204,8 @@ int main(int argc, char ** argv)
// typedef Kokkos::Cuda MemSpace;
// typedef Kokkos::RangePolicy<ExecSpace> range_policy;
ViewMatrixType elai( "elai", n_months, n_pfts );
ViewMatrixType esai( "esai", n_months, n_pfts );
ViewMatrixType elai( "elai", n_grid_cells, n_pfts );
ViewMatrixType esai( "esai", n_grid_cells, n_pfts );
ViewMatrixType::HostMirror h_elai = Kokkos::create_mirror_view( elai );
ViewMatrixType::HostMirror h_esai = Kokkos::create_mirror_view( esai );
......@@ -287,12 +287,12 @@ int main(int argc, char ** argv)
Kokkos::deep_copy( qflx_rain_grnd, h_qflx_rain_grnd);
Kokkos::deep_copy( h2o_can, h_h2o_can);
double* end = &h_h2o_can(n_grid_cells-1, n_pfts-1) ;
std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water" << std::endl;
auto min_max = std::minmax_element(&h_h2o_can(0,0), &h_h2o_can(n_grid_cells, n_pfts));//h2o_can1.begin(), h2o_can1.end());
auto min_max = std::minmax_element(&h_h2o_can(0,0), end+1);//h2o_can1.begin(), h2o_can1.end());
std::cout << std::setprecision(16)
<< 0 << "\t" << std::accumulate(&h_h2o_can(0,0), &h_h2o_can(n_grid_cells, n_pfts), 0.) //h2o_can1.begin(), h2o_can1.end(), 0.)
<< 0 << "\t" << std::accumulate(&h_h2o_can(0,0), end+1, 0.) //h2o_can1.begin(), h2o_can1.end(), 0.)
<< "\t" << *min_max.first
<< "\t" << *min_max.second << std::endl;
......@@ -335,9 +335,9 @@ int main(int argc, char ** argv)
}
});
auto min_max = std::minmax_element(&h_h2o_can(0,0), &h_h2o_can(n_grid_cells, n_pfts));//h2o_can1.begin(), h2o_can1.end());
auto min_max = std::minmax_element(&h_h2o_can(0,0), end+1);//h2o_can1.begin(), h2o_can1.end());
std::cout << std::setprecision(16)
<< t+1 << "\t" << std::accumulate(&h_h2o_can(0,0), &h_h2o_can(n_grid_cells, n_pfts), 0.)//h2o_can1.begin(), h2o_can1.end(), 0.)
<< t+1 << "\t" << std::accumulate(&h_h2o_can(0,0), end+1, 0.)//h2o_can1.begin(), h2o_can1.end(), 0.)
<< "\t" << *min_max.first
<< "\t" << *min_max.second << std::endl;
......
......@@ -37,7 +37,7 @@ using VectorColumnInt = VectorStatic<n_grid_cells,int>;
} // namespace
namespace ELM {
void CanopyHydrology_Interception(double dtime,
KOKKOS_INLINE_FUNCTION void CanopyHydrology_Interception(double dtime,
const double& forc_rain,
const double& forc_snow,
const double& irrig_rate,
......@@ -181,7 +181,7 @@ void CanopyHydrology_Interception(double dtime,
namespace ELM {
void CanopyHydrology_FracWet(const int& frac_veg_nosno,
KOKKOS_INLINE_FUNCTION void CanopyHydrology_FracWet(const int& frac_veg_nosno,
const double& h2ocan,
const double& elai,
const double& esai,
......@@ -217,7 +217,7 @@ void CanopyHydrology_FracWet(const int& frac_veg_nosno,
namespace ELM {
template<typename Array_d>
void CanopyHydrology_SnowWater(const double& dtime,
KOKKOS_INLINE_FUNCTION void CanopyHydrology_SnowWater(const double& dtime,
const double& qflx_floodg,
const int& ltype,
const int& ctype,
......@@ -453,7 +453,7 @@ void CanopyHydrology_SnowWater(const double& dtime,
namespace ELM {
void CanopyHydrology_FracH2OSfc(const double& dtime,
KOKKOS_INLINE_FUNCTION void CanopyHydrology_FracH2OSfc(const double& dtime,
const double& min_h2osfc,
const int& ltype,
const double& micro_sigma,
......@@ -561,6 +561,7 @@ int main(int argc, char ** argv)
const double micro_sigma = 0.1;
const double min_h2osfc = 1.0e-8;
const double n_melt = 0.7;
double qflx_floodg = 0.0;
Kokkos::initialize( argc, argv );
{
......@@ -571,8 +572,8 @@ int main(int argc, char ** argv)
typedef Kokkos::View<int*> ViewVectorType1;
// ELM::Utils::MatrixState elai;
// ELM::Utils::MatrixState esai;
ViewMatrixType elai( "elai", n_months, n_pfts );
ViewMatrixType esai( "esai", n_months, n_pfts );
ViewMatrixType elai( "elai", n_grid_cells, n_pfts );
ViewMatrixType esai( "esai", n_grid_cells, n_pfts );
ViewMatrixType::HostMirror h_elai = Kokkos::create_mirror_view( elai );
ViewMatrixType::HostMirror h_esai = Kokkos::create_mirror_view( esai );
ELM::Utils::read_phenology("../links/surfacedataWBW.nc", n_months, n_pfts, 0, h_elai, h_esai);
......@@ -586,8 +587,10 @@ int main(int argc, char ** argv)
ViewMatrixType::HostMirror h_forc_snow = Kokkos::create_mirror_view( forc_snow );
ViewMatrixType::HostMirror h_forc_air_temp = Kokkos::create_mirror_view( forc_air_temp );
const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 0, n_grid_cells, h_forc_rain, h_forc_snow, h_forc_air_temp);
ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
double qflx_floodg = 0.0;
// ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
ViewMatrixType forc_irrig( "forc_irrig", n_max_times,n_grid_cells );
ViewMatrixType::HostMirror h_forc_irrig = Kokkos::create_mirror_view( forc_irrig );
// mesh input (though can also change as snow layers evolve)
......@@ -697,16 +700,19 @@ int main(int argc, char ** argv)
// std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water" << std::endl;
// auto min_max = std::minmax_element(h2ocan.begin(), h2ocan.end());
// auto min_max = std::minmax_element(&h_h2ocan(0,0), end1+1);
// std::cout << std::setprecision(16)
// << 0 << "\t" << std::accumulate(h2ocan.begin(), h2ocan.end(), 0.)
// << 0 << "\t" << std::accumulate(&h_h2ocan(0,0), end1+1, 0.)
// << "\t" << *min_max.first
// << "\t" << *min_max.second << std::endl;
Kokkos::deep_copy( elai, h_elai);
Kokkos::deep_copy( esai, h_esai);
Kokkos::deep_copy( forc_rain, h_forc_rain);
Kokkos::deep_copy( forc_snow, h_forc_snow);
Kokkos::deep_copy( forc_irrig, h_forc_irrig);
Kokkos::deep_copy( forc_air_temp, h_forc_air_temp);
Kokkos::deep_copy( z, h_z);
Kokkos::deep_copy( zi, h_zi);
......@@ -738,6 +744,24 @@ int main(int argc, char ** argv)
Kokkos::deep_copy( frac_sno_eff, h_frac_sno_eff);
Kokkos::deep_copy( frac_sno, h_frac_sno);
double* end1 = &h_h2ocan(n_grid_cells-1, n_pfts-1) ;
double* end2 = &h_h2osno(n_grid_cells-1) ;
double* end3 = &h_frac_h2osfc(n_grid_cells-1) ;
std::cout << "Time\t Total Canopy Water\t Min Water\t Max Water\t Total Snow\t Min Snow\t Max Snow\t Avg Frac Sfc\t Min Frac Sfc\t Max Frac Sfc" << std::endl;
auto min_max_water = std::minmax_element(&h_h2ocan(0,0), end1+1);
auto sum_water = std::accumulate(&h_h2ocan(0,0), end1+1, 0.);
auto min_max_snow = std::minmax_element(&h_h2osno(0), end2+1);
auto sum_snow = std::accumulate(&h_h2osno(0), end2+1, 0.);
auto min_max_frac_sfc = std::minmax_element(&h_frac_h2osfc(0), end3+1);
auto avg_frac_sfc = std::accumulate(&h_frac_h2osfc(0), end3+1, 0.) / (end3+1 - &h_frac_h2osfc(0));
std::cout << std::setprecision(16)
<< 0 << "\t" << sum_water << "\t" << *min_max_water.first << "\t" << *min_max_water.second
<< "\t" << sum_snow << "\t" << *min_max_snow.first << "\t" << *min_max_snow.second
<< "\t" << avg_frac_sfc << "\t" << *min_max_frac_sfc.first << "\t" << *min_max_frac_sfc.second << std::endl;
// main loop
// -- the timestep loop cannot/should not be parallelized
......@@ -779,16 +803,16 @@ int main(int argc, char ** argv)
// Column level operations
double* qpatch = &qflx_snow_grnd_patch(n_grid_cells-1, n_pfts-1);
// NOTE: this is effectively an accumulation kernel/task! --etc
qflx_snow_grnd_col(g) = std::accumulate(&qflx_snow_grnd_patch(0,0),
&qflx_snow_grnd_patch(n_grid_cells, n_pfts), 0.);
qflx_snow_grnd_col(g) = std::accumulate(&qflx_snow_grnd_patch(0,0), qpatch+1, 0.);
// Calculate ?water balance? on the snow column, adding throughfall,
// removing melt, etc.
//
// local outputs
int newnode;
ELM::CanopyHydrology_SnowWater(dtime, qflx_floodg,
ltype, ctype, urbpoi, do_capsnow, oldfflag,
forc_air_temp(t,g), t_grnd(g),
......@@ -810,12 +834,26 @@ int main(int argc, char ** argv)
}); // end grid cell loop
// auto min_max = std::minmax_element(h2ocan.begin(), h2o