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

adds c variant of fortran full-module test

parent fe712b85
......@@ -28,7 +28,78 @@ void CanopyHydrology_Interception(double dtime,
&qflx_prec_intr, &qflx_irrig, &qflx_prec_grnd,
&qflx_snwcp_liq, &qflx_snwcp_ice, & qflx_snow_grnd_patch, &qflx_rain_grnd);
}
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);
}
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
......
......@@ -25,7 +25,57 @@ extern "C" {
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);
}
......
subroutine CanopyHydrology_FracH2OSfc( dtime, min_h2osfc, &
itype, micro_sigma, h2osno, &
ltype, micro_sigma, h2osno, &
h2osfc, h2osoi_liq, frac_sno, frac_sno_eff, &
qflx_h2osfc2topsoi, frac_h2osfc, &
no_update)
use shr_kind_mod, only : &
r8 => shr_kind_r8, &
i4 => shr_kind_i4
use landunit_varcon , only : istcrop, istice, istwet, istsoil, istice_mec
i4 => shr_kind_i4, &
bool => shr_kind_bool
use landunit_varcon, only : istcrop, istice, istwet, istsoil, istice_mec
implicit none
!
! interface
real(r8), intent(in) :: dtime, min_h2osfc, micro_sigma, h2osno
integer(i4), intent(in) :: itype
integer(i4), intent(in), optional :: no_update
integer(i4), intent(in) :: ltype
real(r8), intent(inout) :: h2osfc, h2osoi_liq, frac_sno, frac_sno_eff
real(r8), intent(out) :: qflx_h2osfc2topsoi, frac_h2osfc
logical(bool), intent(in), optional :: no_update
!local
real(r8), parameter :: shr_const_pi=4.0d0*atan(1.0d0)
logical(bool) :: no_update_l
integer :: l
real(r8):: d,fd,dfdd ! temporary variable for frac_h2oscs iteration
real(r8):: sigma ! microtopography pdf sigma in mm
if (.not. present(no_update) ) then
no_update_l = .false.
else
no_update_l = no_update
end if
qflx_h2osfc2topsoi = 0._r8
! h2osfc only calculated for soil vegetated land units
if ( itype == istsoil .or. itype == istcrop) then
if ( ltype == istsoil .or. ltype == istcrop) then
! Use newton-raphson method to iteratively determine frac_h20sfc
! based on amount of surface water storage (h2osfc) and
......@@ -55,7 +66,7 @@ subroutine CanopyHydrology_FracH2OSfc( dtime, min_h2osfc, &
h2osfc=0._r8
endif
if (.not. present(no_update)) then
if (.not.no_update_l) then
! adjust fh2o, fsno when sum is greater than zero
if (frac_sno > (1._r8 - frac_h2osfc) .and. h2osno > 0) then
......
subroutine CanopyHydrology_SnowWater( dtime, &
qflx_floodg, &
ctype, ltype, urbpoi, do_capsnow, oldfflag, &
ltype, ctype, urbpoi, do_capsnow, oldfflag, &
forc_t, t_grnd, qflx_snow_grnd_col, qflx_snow_melt, n_melt, frac_h2osfc, & ! forcing
snow_depth, h2osno, int_snow, 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, &
frac_sno_eff, frac_sno)
frac_sno, frac_sno_eff)
use shr_kind_mod, only : &
r8 => shr_kind_r8, &
......@@ -20,28 +20,31 @@ subroutine CanopyHydrology_SnowWater( dtime, &
implicit none
!
! parameters
real(r8), parameter :: rpi=4.0d0*atan(1.0d0)
real(r8), parameter :: tfrz=273.15
real(r8) :: zlnd = 0.01_r8
real(r8), parameter :: zlnd = 0.01_r8
!
! interface variables
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(i4), intent(in) :: ltype, ctype, oldfflag
logical(bool), intent(in) :: urbpoi, do_capsnow
real(r8), intent(in) :: forc_t, t_grnd, qflx_snow_grnd_col, qflx_snow_melt, n_melt, frac_h2osfc
integer, intent(inout) :: snl
real(r8), intent(inout) :: snow_depth , h2osno, int_snow
real(r8), intent(inout), dimension(-nlevsno+1:0) :: swe_old
real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
real(r8), intent(inout), dimension(-nlevsno+1:0) :: dz
real(r8), intent(inout), dimension(-nlevsno+1:0) :: t_soisno, frac_iceold
integer, intent(inout) :: snl
real(r8), intent(inout), dimension(-nlevsno+1:0) :: dz, z, zi
integer, intent(out) :: newnode
real(r8), intent(out) :: qflx_floodc
real(r8), intent(out) :: qflx_snow_h2osfc
real(r8), intent(out), dimension(-nlevsno+1:0) :: swe_old
real(r8), intent(out), dimension(-nlevsno+1:0) :: z, zi
real(r8), intent(out), dimension(-nlevsno+1:0) :: t_soisno, frac_iceold
real(r8), intent(out) :: frac_sno_eff, frac_sno
! local variables
......
#include <array>
#include <sstream>
#include <iterator>
#include <exception>
#include <string>
#include <stdlib.h>
#include <cstring>
#include <vector>
#include <iostream>
#include <iomanip>
#include <numeric>
#include "utils.hh"
#include "readers.hh"
#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;
static const int n_levels_snow = 5;
using MatrixStatePFT = MatrixStatic<n_grid_cells, n_pfts>;
using MatrixStateSoilColumn = MatrixStatic<n_grid_cells, n_levels_snow>;
using MatrixForc = MatrixStatic<n_max_times,n_grid_cells>;
using VectorColumn = VectorStatic<n_grid_cells>;
using VectorColumnInt = VectorStatic<n_grid_cells,int>;
} // namespace
} // namespace
int main(int argc, char ** argv)
{
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;
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;
const int oldfflag = 0;
// fixed magic parameters for fracH2Osfc
const double micro_sigma = 0.1;
const double min_h2osfc = 1.0e-8;
// NOTE: confusing parameters that need to be fixed --etc
const double qflx_snow_melt = 0.;
const double n_melt = 0.7;
// phenology input
ELM::Utils::MatrixStatePFT elai;
ELM::Utils::MatrixStatePFT 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 input
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::MatrixForc forc_irrig; forc_irrig = 0.;
double qflx_floodg = 0.0;
// mesh input (though can also change as snow layers evolve)
//
// NOTE: in a real case, these would be populated, but we don't actually
// need them to be for these kernels. --etc
auto z = ELM::Utils::MatrixStateSoilColumn(); z = 0.;
auto zi = ELM::Utils::MatrixStateSoilColumn(); zi = 0.;
auto dz = ELM::Utils::MatrixStateSoilColumn(); dz = 0.;
// state variables that require ICs and evolve (in/out)
auto h2ocan = ELM::Utils::MatrixStatePFT(); h2ocan = 0.;
auto swe_old = ELM::Utils::MatrixStateSoilColumn(); swe_old = 0.;
auto h2osoi_liq = ELM::Utils::MatrixStateSoilColumn(); h2osoi_liq = 0.;
auto h2osoi_ice = ELM::Utils::MatrixStateSoilColumn(); h2osoi_ice = 0.;
auto t_soisno = ELM::Utils::MatrixStateSoilColumn(); t_soisno = 0.;
auto frac_iceold = ELM::Utils::MatrixStateSoilColumn(); frac_iceold = 0.;
auto t_grnd = ELM::Utils::VectorColumn(); t_grnd = 0.;
auto h2osno = ELM::Utils::VectorColumn(); h2osno = 0.;
auto snow_depth = ELM::Utils::VectorColumn(); snow_depth= 0.;
auto snl = ELM::Utils::VectorColumnInt(); snl = 0; // note this tracks the snow_depth
auto h2osfc = ELM::Utils::VectorColumn(); h2osfc = 0.1;
auto frac_h2osfc = ELM::Utils::VectorColumn(); frac_h2osfc = 0.;
// output fluxes by pft
auto qflx_prec_intr = ELM::Utils::MatrixStatePFT();
auto qflx_irrig = ELM::Utils::MatrixStatePFT();
auto qflx_prec_grnd = ELM::Utils::MatrixStatePFT();
auto qflx_snwcp_liq = ELM::Utils::MatrixStatePFT();
auto qflx_snwcp_ice = ELM::Utils::MatrixStatePFT();
auto qflx_snow_grnd_patch = ELM::Utils::MatrixStatePFT();
auto qflx_rain_grnd = ELM::Utils::MatrixStatePFT();
// NOTE: 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(); int_snow = 0.;
// output fluxes, state by the column
auto qflx_snow_grnd_col = ELM::Utils::VectorColumn();
auto qflx_snow_h2osfc = ELM::Utils::VectorColumn();
auto qflx_h2osfc2topsoi = ELM::Utils::VectorColumn();
auto qflx_floodc = ELM::Utils::VectorColumn();
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" << 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;
// 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) {
// PFT level operations
for (size_t p = 0; p != n_pfts; ++p) {
//
// Calculate interception
//
// 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::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,
h2ocan(g,p), n_irrig_steps_left,
qflx_prec_intr(g,p), qflx_irrig(g,p), qflx_prec_grnd(g,p),
qflx_snwcp_liq(g,p), qflx_snwcp_ice(g,p),
qflx_snow_grnd_patch(g,p), qflx_rain_grnd(g,p));
//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), h2ocan(g,p), qflx_prec_intr[g]);
//
// Calculate fraction of ?LAI? that is wet vs dry.
//
// NOTE: this currently punts on what to do with the fwet/fdry variables.
// Surely they should be something, as such this is dead code.
// --etc
double fwet = 0., fdry = 0.;
ELM::CanopyHydrology_FracWet(frac_veg_nosno, h2ocan(g,p), elai(g,p), esai(g,p), dewmx, fwet, fdry);
} // end PFT loop
// Column level operations
// NOTE: this is effectively an accumulation kernel/task! --etc
qflx_snow_grnd_col[g] = std::accumulate(qflx_snow_grnd_patch[g].begin(),
qflx_snow_grnd_patch[g].end(), 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),
qflx_snow_grnd_col[g], qflx_snow_melt, n_melt, frac_h2osfc[g],
snow_depth[g], h2osno[g], int_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,
qflx_floodc[g], qflx_snow_h2osfc[g], frac_sno_eff[g], frac_sno[g]);
// Calculate Fraction of Water to the Surface?
//
// FIXME: Fortran black magic... h2osoi_liq is a vector, but the
// interface specifies a single double. For now passing the 0th
// entry. --etc
ELM::CanopyHydrology_FracH2OSfc(dtime, min_h2osfc, ltype, micro_sigma,
h2osno[g], h2osfc[g], h2osoi_liq[g][0], frac_sno[g], frac_sno_eff[g],
qflx_h2osfc2topsoi[g], frac_h2osfc[g]);
} // end grid cell loop
// 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;
} // end timestep loop
return 0;
}
......@@ -8,10 +8,12 @@ CPPFLAGS += -I$(CSRCDIR) -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
......@@ -30,6 +32,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 readers.hh utils.hh
$(LINKER) $(FCFLAGS) -L$(OBJECT)fortran -lelm -L$(NETCDF_ROOT)/lib -lnetcdf -L$(STD_LIB_ROOT)/lib -lc++ $< -o $@
......
......@@ -7,12 +7,34 @@ namespace ELM {
namespace Utils {
template<size_t N, typename T=double>
class VectorStatic {
public:
VectorStatic() {}
VectorStatic(T t) { *this = t; }
T& operator()(size_t i) { return d_[i]; }
const T& operator()(size_t i) const { return d_[i]; }
T& operator[](size_t i) { return d_[i]; }
const T& operator[](size_t i) const { return d_[i]; }
void operator=(T t) {
for (size_t i=0; i!=N; ++i) {
d_[i] = t;
}
}
private:
std::array<T,N> d_;
};
template<size_t ROW, size_t COL, typename T=double>
class MatrixStatic {
public:
MatrixStatic() {}
MatrixStatic(T t) { *this = t; }
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]; }
......
......@@ -135,15 +135,16 @@ program test_CanopyHydrology_module
call CanopyHydrology_SnowWater( dtime, &
qflx_floodg, &
ctype, ltype, urbpoi, do_capsnow, oldfflag, &
ltype, ctype, urbpoi, do_capsnow, oldfflag, &
forc_t, t_grnd, qflx_snow_grnd_col, qflx_snow_melt, n_melt, frac_h2osfc, & ! forcing
snow_depth, h2osno, int_snow, 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, &
frac_sno_eff, frac_sno)
frac_sno, frac_sno_eff)
! FIXME: Fortran black magic... h2osoi_liq is a vector, but
! the interface specifies a single double. --etc
call CanopyHydrology_FracH2OSfc( dtime, min_h2osfc, &
ltype, micro_sigma, h2osno, &
h2osfc, h2osoi_liq, frac_sno, frac_sno_eff, &
......
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