CanopyHydrology_module.cc 9.07 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include <array>
#include <sstream>
#include <iterator>
#include <exception>
#include <string>
#include <stdlib.h>
#include <cstring>
#include <vector>
#include <iostream>
#include <iomanip>
#include <numeric>
12
#include <algorithm>
13
14
#include "utils.hh"
#include "readers.hh"
15

16
#include "CanopyHydrology.hh"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56


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;

Ethan Coon's avatar
Ethan Coon committed
57
58
  // fixed magic parameters for SnowWater
  const double qflx_snow_melt = 0.;
59
60

  // fixed magic parameters for fracH2Osfc  
Ethan Coon's avatar
Ethan Coon committed
61
  const int oldfflag = 0;
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
  const double micro_sigma = 0.1;
  const double min_h2osfc = 1.0e-8;
  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
Ethan Coon's avatar
Ethan Coon committed
85
86
87
  auto z = ELM::Utils::MatrixStateSoilColumn(0.);
  auto zi = ELM::Utils::MatrixStateSoilColumn(0.);
  auto dz = ELM::Utils::MatrixStateSoilColumn(0.);
88
89

  // state variables that require ICs and evolve (in/out)
Ethan Coon's avatar
Ethan Coon committed
90
91
92
93
94
95
96
97
98
99
100
101
102
  auto h2ocan = ELM::Utils::MatrixStatePFT(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 snow_depth = ELM::Utils::VectorColumn(0.);
  auto snl = ELM::Utils::VectorColumnInt(0.); // note this tracks the snow_depth

  auto h2osfc = ELM::Utils::VectorColumn(0.);
  auto frac_h2osfc = ELM::Utils::VectorColumn(0.);
103
104
105
106
107
108
109
110
111
112
113

  
  // 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();

Ethan Coon's avatar
Ethan Coon committed
114
  // FIXME: I have no clue what this is... it is inout on WaterSnow.  For now I
115
  // am guessing the data structure. Ask Scott.  --etc
Ethan Coon's avatar
Ethan Coon committed
116
  auto int_snow = ELM::Utils::VectorColumn(0.);
117
118
119
120
121
122
123
124
125
126
  
  // 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();
  
127
128
129
130
  {
    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.);
131

132
133
134
135
136
137
138
139
140
141
142
    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 << 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;
  }
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
  
  // 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]);

        //
170
        // Calculate fraction of LAI that is wet vs dry.
171
        //
Ethan Coon's avatar
Ethan Coon committed
172
        // FIXME: this currently punts on what to do with the fwet/fdry variables.
173
        // Surely they should be something, as such this is dead code.
174
        // By the PFT?
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
        // --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

211
212
213
214
215
    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.);
216

217
218
219
220
221
222
223
    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;
224
225
226
  } // end timestep loop
  return 0;
}