Commit f3c9a2c0 authored by Ethan Coon's avatar Ethan Coon
Browse files

working multiple-grid-cell and multiple-pft example for Kern1

parent 3c5e40f8
subroutine CanopyHydrologyKern1( dtime, &
forc_rain, forc_snow, irrig_rate, &
forc_rain, forc_snow, forc_irrig, &
ltype, ctype, urbpoi, do_capsnow, &
elai, esai, dewmx, frac_veg_nosno, &
h2ocan, n_irrig_steps_left, &
......@@ -31,7 +31,7 @@ subroutine CanopyHydrologyKern1( dtime, &
real(r8), intent(out) :: qflx_irrig
real(r8), intent(out) :: qflx_prec_grnd
real(r8), intent(out) :: qflx_snwcp_liq, qflx_snwcp_ice, qflx_snow_grnd_patch, qflx_rain_grnd
real(r8), intent(in) :: irrig_rate
real(r8), intent(in) :: forc_irrig
!local variables
real(r8) :: fpi, xrun, h2ocanmx
......@@ -44,7 +44,7 @@ subroutine CanopyHydrologyKern1( dtime, &
! Canopy interception and precipitation onto ground surface
! Add precipitation to leaf water
if (ltype==istsoil .or. ltype==istwet .or. urbpoi .or. &
ltype==istcrop) then
......@@ -58,6 +58,7 @@ subroutine CanopyHydrologyKern1( dtime, &
if (ctype /= icol_sunwall .and. ctype /= icol_shadewall) then
if (frac_veg_nosno == 1 .and. (forc_rain + forc_snow) > 0._r8) then
! determine fraction of input precipitation that is snow and rain
fracsnow = forc_snow/(forc_snow + forc_rain)
fracrain = forc_rain/(forc_snow + forc_rain)
......@@ -79,6 +80,8 @@ 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
! Water storage of intercepted precipitation and dew
h2ocan = max(0._r8, h2ocan + dtime*qflx_prec_intr)
......@@ -129,7 +132,7 @@ subroutine CanopyHydrologyKern1( dtime, &
! Determine whether we're irrigating here; set qflx_irrig appropriately
if (n_irrig_steps_left > 0) then
qflx_irrig = irrig_rate
qflx_irrig = forc_irrig
n_irrig_steps_left = n_irrig_steps_left - 1
else
qflx_irrig = 0._r8
......
......@@ -60,7 +60,9 @@ program CanopyHydrology_kern1_multiple
! end setup stage
! -----------------------------
h2ocan_pft = 0.0d0
print*, "Time", "Total Canopy Water", "Min Water", "MaxWater"
print*, 0, sum(h2ocan_pft), minval(h2ocan_pft), maxval(h2ocan_pft)
do itime=1,28*48 ! February is shortest month
do g=1,ngrcs ! grid cell loop
......@@ -83,21 +85,26 @@ program CanopyHydrology_kern1_multiple
qflx_prec_intr, qflx_irrig, qflx_prec_grnd, &
qflx_snwcp_liq, qflx_snwcp_ice, qflx_snow_grnd_patch, qflx_rain_grnd)
!print*, g, p, forc_rain, forc_snow, elai, esai, h2ocan, qflx_prec_intr
h2ocan_pft(p,g) = h2ocan
qflx_snow_grnd_col = qflx_snow_grnd_col + qflx_snow_grnd_patch
end do ! PFT loop
end do ! grid cell loop
end do ! grid cell loop
print*, itime, sum(h2ocan_pft), minval(h2ocan_pft), maxval(h2ocan_pft)
end do ! time loop
print*, "Final canopy water:"
print*, "Grid Cell", "PFT", "H2O"
do g=1,ngrcs
do p=1,npfts
print*, g, p, h2ocan_pft(p,g)
end do
end do
! print*, ""
! print*, "Final canopy water:"
! print*, "Grid Cell", "PFT", "H2O"
! do g=1,ngrcs
! do p=1,npfts
! print*, g, p, h2ocan_pft(p,g)
! end do
! end do
stop
contains
......
.PHONY: links library CanopyHydrology_kern1_single
SRCDIR = ../src/fortran/
OBJECT = ../src/
default: all
include $(OBJECT)config/Makefile.config
FCFLAGS += -I$(SRCDIR) -I$(SRCDIR)$(ELM_UTILS_DIR) -I$(NETCDF_ROOT)/include
all: links library CanopyHydrology_kern1_single
TESTS = test_CanopyHydrology_kern1_single \
test_CanopyHydrology_kern1_multiple
EXEC_TESTS = CanopyHydrology_kern1_single \
CanopyHydrology_kern1_multiple
.PHONY: links library test
default: all
all: links library $(TESTS)
test: $(EXEC_TESTS)
echo "done"
CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
./test_CanopyHydrology_kern1_single
......@@ -16,19 +27,16 @@ CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
CanopyHydrology_kern1_multiple: test_CanopyHydrology_kern1_multiple
./test_CanopyHydrology_kern1_multiple
test_%: %.o
$(LINKER) $(FCFLAGS) -L$(OBJECT)fortran -lelm -L$(NETCDF_ROOT)/lib -lnetcdff $< -o $@
test: test_CanopyHydrology_kern1_single test_CanopyHydrology_kern1_multiple
echo "done"
clean:
@$(ELM_CLEAN)
$(RM) test_*
allclean:
@$(ELM_CLEAN)
$(RM) test_*
......
#include <array>
#include <sstream>
#include <iterator>
#include <exception>
#include <string>
#include <stdlib.h>
#include <cstring>
#include <vector>
#include <iostream>
#include <numeric>
#include "utils.hh"
#include "readers.hh"
using namespace std;
#include "CanopyHydrology.hh"
#define handle_error( status, what ) \
do { \
if ( status ) \
{ \
std::cout \
<< __FILE__ \
<< ':' \
<< __LINE__ \
<< ':' \
<< what \
<< " failed with rc = " \
<< status \
<< ':' \
<< nc_strerror( status ) \
<< '\n' ; \
abort() ; \
} \
} while ( 0 )
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
// fixed magic parameters for now
const int ctype = 1;
const int ltype = 1;
const bool urbpoi = false;
const bool do_capsnow = false;
const int frac_veg_nosno = 1;
int n_irrig_steps_left = 0;
const double dewmx = 0.1;
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::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);
const int n_times = ELM::Utils::read_forcing("links/forcing", n_max_times, n_grid_cells, forc_rain, forc_snow, forc_air_temp);
ELM::Utils::Matrix<> forc_irrig(n_max_times, n_grid_cells); 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.);
// output state by the pft
auto h2o_can = ELM::Utils::Matrix<>(n_grid_cells, n_pfts); 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());
std::cout << 0 << "\t" << std::accumulate(h2o_can.begin(), h2o_can.end(), 0.)
<< "\t" << *min_max.first
<< "\t" << *min_max.second << std::endl;
// main loop
// -- the timestep loop cannot/should not be parallelized
for (size_t t = 0; t != n_times; ++t) {
// grid cell and/or pft loop can be parallelized
for (size_t g = 0; g != n_grid_cells; ++g) {
for (size_t p = 0; p != n_pfts; ++p) {
// 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,
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,
h2o_can(g,p), n_irrig_steps_left,
qflx_prec_intr[g], qflx_irrig[g], qflx_prec_grnd[g],
qflx_snwcp_liq[g], qflx_snwcp_ice[g],
qflx_snow_grnd_patch[g], qflx_rain_grnd[g]);
//printf("%i %i %16.8g %16.8g %16.8g %16.8g %16.8g %16.8g\n", g, p, forc_rain(t,g), forc_snow(t,g), elai(g,p), esai(g,p), h2o_can(g,p), qflx_prec_intr[g]);
}
}
auto min_max = std::minmax_element(h2o_can.begin(), h2o_can.end());
std::cout << t << "\t" << std::accumulate(h2o_can.begin(), h2o_can.end(), 0.)
<< "\t" << *min_max.first
<< "\t" << *min_max.second << std::endl;
}
return 0;
}
......@@ -13,132 +13,132 @@ using namespace std;
#include "CanopyHydrology.hh"
#define handle_error( status, what ) \
do { \
if ( status ) \
{ \
std::cout \
<< __FILE__ \
<< ':' \
<< __LINE__ \
<< ':' \
<< what \
<< " failed with rc = " \
<< status \
<< ':' \
<< nc_strerror( status ) \
<< '\n' ; \
abort() ; \
} \
} while ( 0 )
#define handle_error( status, what ) \
do { \
if ( status ) \
{ \
std::cout \
<< __FILE__ \
<< ':' \
<< __LINE__ \
<< ':' \
<< what \
<< " failed with rc = " \
<< status \
<< ':' \
<< nc_strerror( status ) \
<< '\n' ; \
abort() ; \
} \
} while ( 0 )
int main(int argc, char ** argv)
{
int ctype = 1;
const int nmonths = 12;
const int npfts = 17;
bool urbpoi = false,do_capsnow=false;
const double dtime = 1800.0;
double elai = 0, esai = 0, dewmx = 0, forc_snow = 0, forc_rain = 0, qflx_prec_intr = 0, qflx_irrig = 0, qflx_prec_grnd = 0, qflx_snwcp_liq = 0, qflx_snwcp_ice = 0, qflx_snow_grnd_patch = 0, qflx_rain_grnd = 0,irrig_rate = 0 ,h2ocan = 0;
int frac_veg_nosno = 1;
int ipft,itime ;
size_t ntimes=0;
const int ltype = 1;
double monthly_lai[nmonths][npfts] = {};
double monthly_sai[nmonths][npfts] = {};
int n_irrig_steps_left = 0;
double * total_precip = NULL;
int ctype = 1;
const int nmonths = 12;
const int npfts = 17;
bool urbpoi = false,do_capsnow=false;
const double dtime = 1800.0;
double elai = 0, esai = 0, dewmx = 0, forc_snow = 0, forc_rain = 0, qflx_prec_intr = 0, qflx_irrig = 0, qflx_prec_grnd = 0, qflx_snwcp_liq = 0, qflx_snwcp_ice = 0, qflx_snow_grnd_patch = 0, qflx_rain_grnd = 0,irrig_rate = 0 ,h2ocan = 0;
int frac_veg_nosno = 1;
int ipft,itime ;
size_t ntimes=0;
const int ltype = 1;
double monthly_lai[nmonths][npfts] = {};
double monthly_sai[nmonths][npfts] = {};
int n_irrig_steps_left = 0;
double * total_precip = NULL;
int ncid= 0, varid= 0, dimid= 0, status= 0;
int ncid= 0, varid= 0, dimid= 0, status= 0;
/*! use netcdf*/
/*! use netcdf*/
/*! This will be the netCDF ID for the file and data variable.*/
/*! This will be the netCDF ID for the file and data variable.*/
status = nc_open("links/surfacedataWBW.nc", NC_NOWRITE, &ncid);
handle_error( status, "nc_create" ) ;
status = nc_open("links/surfacedataWBW.nc", NC_NOWRITE, &ncid);
handle_error( status, "nc_create" ) ;
status = nc_inq_varid(ncid, "MONTHLY_LAI", &varid);
handle_error( status, "nc_inq_varid" ) ;
status = nc_inq_varid(ncid, "MONTHLY_LAI", &varid);
handle_error( status, "nc_inq_varid" ) ;
array<size_t,4> start{ nmonths,npfts, 1, 1 };
array<size_t,4> count{ 0,0,0,0 };
array<size_t,4> start{ nmonths,npfts, 1, 1 };
array<size_t,4> count{ 0,0,0,0 };
status = nc_get_vara_double(ncid, varid, count.data(), start.data(), *monthly_lai);
handle_error( status, "nc_get_vara_double monthly_lai" ) ;
status = nc_get_vara_double(ncid, varid, count.data(), start.data(), *monthly_lai);
handle_error( status, "nc_get_vara_double monthly_lai" ) ;
status = nc_inq_varid(ncid, "MONTHLY_SAI", &varid);
handle_error( status, "nc_inq_varid" ) ;
/*!do ipft=1,npfts */
/*! print *, monthly_lai(ipft,12), monthly_lai(ipft,3), monthly_lai(ipft,6) */
/*!end do */
ipft = 1;
std::cout << "[ " << monthly_lai[12][ipft] << " , " << monthly_lai[3][ipft] << " , " << monthly_lai[6][ipft] << " ]" << std::endl;
status = nc_get_vara_double(ncid, varid, count.data(), start.data(), *monthly_sai);
handle_error( status, "nc_get_vara_double monthly_sai" ) ;
status = nc_inq_varid(ncid, "MONTHLY_SAI", &varid);
handle_error( status, "nc_inq_varid" ) ;
/*!do ipft=1,npfts */
/*! print *, monthly_lai(ipft,12), monthly_lai(ipft,3), monthly_lai(ipft,6) */
/*!end do */
ipft = 1;
std::cout << "[ " << monthly_lai[12][ipft] << " , " << monthly_lai[3][ipft] << " , " << monthly_lai[6][ipft] << " ]" << std::endl;
status = nc_get_vara_double(ncid, varid, count.data(), start.data(), *monthly_sai);
handle_error( status, "nc_get_vara_double monthly_sai" ) ;
status = nc_close(ncid);
handle_error( status, "nc_close1" ) ;
status = nc_close(ncid);
handle_error( status, "nc_close1" ) ;
status = nc_open("links/forcing7.nc", NC_NOWRITE, &ncid);
handle_error( status, "nc_create" ) ;
status = nc_open("links/forcing7.nc", NC_NOWRITE, &ncid);
handle_error( status, "nc_create" ) ;
status = nc_inq_dimid(ncid, "time", &dimid);
handle_error( status, "nc_inq_dimid" ) ;
status = nc_inq_dimid(ncid, "time", &dimid);
handle_error( status, "nc_inq_dimid" ) ;
status = nc_inq_dimlen(ncid, dimid, &ntimes );
handle_error( status, "nc_inq_dimlen" ) ;
array<size_t,3> start3{ntimes, 1, 1};
array<size_t,3> count3{ 0,0,0 };
total_precip = (double*)malloc(ntimes*1*1*sizeof(double));
status = nc_inq_varid(ncid,"PRECTmms", &varid);
handle_error( status, "nc_inq_varid" ) ;
status = nc_get_vara_double(ncid, varid, count3.data(), start3.data(), total_precip);
handle_error( status, "nc_get_vara_double total_precip" ) ;
status = nc_close(ncid);
handle_error( status, "nc_close2" ) ;
//cout.precision(10);
h2ocan = 0.0;
for(itime = 0; itime < ntimes; itime += 1){
forc_rain = total_precip[itime];
forc_snow = 0.0;
dewmx = 0.1;
elai = monthly_lai[5][7];
esai = monthly_sai[5][7];
frac_veg_nosno = 1;
irrig_rate = 0.0;
n_irrig_steps_left = 0;
urbpoi = false;
do_capsnow = false;
ELM::CanopyHydrologyKern1(dtime, forc_snow, 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);
status = nc_inq_dimlen(ncid, dimid, &ntimes );
handle_error( status, "nc_inq_dimlen" ) ;
array<size_t,3> start3{ntimes, 1, 1};
array<size_t,3> count3{ 0,0,0 };
total_precip = (double*)malloc(ntimes*1*1*sizeof(double));
status = nc_inq_varid(ncid,"PRECTmms", &varid);
handle_error( status, "nc_inq_varid" ) ;
status = nc_get_vara_double(ncid, varid, count3.data(), start3.data(), total_precip);
handle_error( status, "nc_get_vara_double total_precip" ) ;
status = nc_close(ncid);
handle_error( status, "nc_close2" ) ;
//cout.precision(10);
h2ocan = 0.0;
for(itime = 0; itime < ntimes; itime += 1){
forc_rain = total_precip[itime];
forc_snow = 0.0;
dewmx = 0.1;
elai = monthly_lai[5][7];
esai = monthly_sai[5][7];
frac_veg_nosno = 1;
irrig_rate = 0.0;
n_irrig_steps_left = 0;
urbpoi = false;
do_capsnow = false;
ELM::CanopyHydrologyKern1(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);
//printf("[ %d,%E,%4.4f,%E,%E]\n",itime, forc_rain, h2ocan, qflx_prec_grnd, qflx_prec_intr);
std::cout << "[ " << itime << " , " << forc_rain<< " , " << h2ocan<< " , " << qflx_prec_grnd<< " , " << qflx_prec_intr << " ]" << std::endl;
//printf("[ %d,%E,%4.4f,%E,%E]\n",itime, forc_rain, h2ocan, qflx_prec_grnd, qflx_prec_intr);
std::cout << "[ " << itime << " , " << forc_rain<< " , " << h2ocan<< " , " << qflx_prec_grnd<< " , " << qflx_prec_intr << " ]" << std::endl;
}
}
return 0;
return 0;
}
.PHONY: links library CanopyHydrology_kern1_single
SRCDIR = ../src/fortran/
CSRCDIR = ../src/cpp/
OBJECT = ../src/
default: all
include $(OBJECT)config/Makefile.config
CPPFLAGS += -I$(CSRCDIR) -I$(NETCDF_ROOT)/include
all: links library CanopyHydrology_kern1_single
CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
TESTS = test_CanopyHydrology_kern1_single \
test_CanopyHydrology_kern1_multiple
EXEC_TESTS = CanopyHydrology_kern1_single \
CanopyHydrology_kern1_multiple
.PHONY: links library test
default: all
all: links library $(TESTS)
test: $(EXEC_TESTS)
echo "done"
CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single readers.hh utils.hh
CanopyHydrology_kern1_multiple: test_CanopyHydrology_kern1_multiple readers.hh utils.hh
./test_CanopyHydrology_kern1_multiple
test_%: %.o
$(LINKER) $(FCFLAGS) -L$(OBJECT)fortran -lelm -L$(NETCDF_ROOT)/lib -lnetcdf -L$(STD_LIB_ROOT)/lib -lc++ $< -o $@
test: test_CanopyHydrology_kern1_single
./test_CanopyHydrology_kern1_single
clean:
@$(ELM_CLEAN)
......
//! A set of utilities for testing ELM kernels in C++
#ifndef ELM_KERNEL_TEST_NETCDF_HH_
#define ELM_KERNEL_TEST_NETCDF_HH_
#include "netcdf.h"
#define NC_HANDLE_ERROR( status, what ) \
do { \
if ( status ) \
{ \
std::cout \
<< __FILE__ \
<< ':' \
<< __LINE__ \
<< ':' \
<< what \
<< " failed with rc = " \
<< status \
<< ':' \
<< nc_strerror( status ) \
<< '\n' ; \
abort() ; \
} \
} while ( 0 )
namespace ELM {
namespace Utils {
//
// Reads n_grid_cells worth of phenology data from NetCDF files, each with
// n_pft PFTs, into LAI and SAI matrices with a potential offset (in grid
// cells) given by offset.
// -----------------------------------------------------------------------------
template<typename Matrix_t>
void read_phenology(const std::string& fname,
size_t n_grid_cells, size_t n_pfts,
size_t offset,
Matrix_t& lai, Matrix_t& sai) {
int ncid = -1;
auto status = nc_open(fname.c_str(), NC_NOWRITE, &ncid);
NC_HANDLE_ERROR(status, std::string("nc_open")+" \""+fname+"\"");
auto start = std::array<size_t,4>{0,0,0,0};
auto count = std::array<size_t,4>{n_grid_cells, n_pfts, 1, 1};