CanopyHydrology_kern1_multiple.cc 3.91 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include <array>
#include <sstream>
#include <iterator>
#include <exception>
#include <string>
#include <stdlib.h>
#include <cstring>
#include <vector>
#include <iostream>
10
#include <iomanip>
11
12
13
14
15
16
17
#include <numeric>

#include "utils.hh"
#include "readers.hh"
#include "CanopyHydrology.hh"


18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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


35
36
int main(int argc, char ** argv)
{
37
38
39
40
41
  using ELM::Utils::n_months;
  using ELM::Utils::n_pfts;
  using ELM::Utils::n_grid_cells;
  using ELM::Utils::n_max_times;
  
42
43
44
45
46
47
48
49
50
51
52
53
  // 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
54
55
  ELM::Utils::MatrixState elai;
  ELM::Utils::MatrixState esai;
56
57
  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);
58
59

  // forcing state
60
61
62
  ELM::Utils::MatrixForc forc_rain;
  ELM::Utils::MatrixForc forc_snow;
  ELM::Utils::MatrixForc forc_air_temp;
63
  const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 0, n_grid_cells, forc_rain, forc_snow, forc_air_temp);
64

65
  ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
66
67
  
  // output state by the grid cell
68
69
70
71
72
73
74
  auto qflx_prec_intr = ELM::Utils::MatrixState();
  auto qflx_irrig = ELM::Utils::MatrixState();
  auto qflx_prec_grnd = ELM::Utils::MatrixState();
  auto qflx_snwcp_liq = ELM::Utils::MatrixState();
  auto qflx_snwcp_ice = ELM::Utils::MatrixState();
  auto qflx_snow_grnd_patch = ELM::Utils::MatrixState();
  auto qflx_rain_grnd = ELM::Utils::MatrixState();
75
76

  // output state by the pft
77
  auto h2o_can = ELM::Utils::MatrixState(); h2o_can = 0.;
78
79
80

  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());
81
82
  std::cout << std::setprecision(16)
            << 0 << "\t" << std::accumulate(h2o_can.begin(), h2o_can.end(), 0.)
83
84
85
86
87
88
89
90
91
92
93
94
95
            << "\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
96
        ELM::CanopyHydrology_Interception(dtime,
97
98
99
100
                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,
101
102
103
                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));
104
105
106
107
108
        //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());
109
110
    std::cout << std::setprecision(16)
              << t+1 << "\t" << std::accumulate(h2o_can.begin(), h2o_can.end(), 0.)
111
112
113
114
115
116
              << "\t" << *min_max.first
              << "\t" << *min_max.second << std::endl;

  }
  return 0;
}