readers.hh 5.4 KB
Newer Older
1
2
3
4
//! A set of utilities for testing ELM kernels in C++
#ifndef ELM_KERNEL_TEST_NETCDF_HH_
#define ELM_KERNEL_TEST_NETCDF_HH_

5
6
7
8
9
10
#include <string>
#include <array>
#include <vector>
#include <sstream>
#include <iostream>

11
12
13
14
15
16
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
57
58
59
60
61
62
63
64
65
#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};

  // get id and read LAI
  int varid = -1;
  status = nc_inq_varid(ncid, "MONTHLY_LAI", &varid);
  NC_HANDLE_ERROR(status, "nc_inq_varid for LAI");

  std::vector<double> data(n_grid_cells * n_pfts);
  status = nc_get_vara_double(ncid, varid, start.data(), count.data(), data.data());
  NC_HANDLE_ERROR(status, "nc_get_vara_double monthly_lai");

  // copy into array -- note this may thrash cache depending upon order of
  // LAI.  Good for C order, bad for Fortran
  for (int i=offset; i!=offset+n_grid_cells; ++i) {
    for (int j=0; j!=n_pfts; ++j) {
66
      lai[i][j] = data[(i-offset)*n_pfts + j];
67
68
69
70
71
72
73
74
75
76
77
78
79
80
    }
  }
  
  
  // get id and read SAI directly into the array
  status = nc_inq_varid(ncid, "MONTHLY_SAI", &varid);
  NC_HANDLE_ERROR( status, "nc_inq_varid for SAI" ) ;
  
  status = nc_get_vara_double(ncid, varid, start.data(), count.data(), data.data());
  NC_HANDLE_ERROR(status, "nc_get_vara_double monthly_sai");
  // copy into array -- note this may thrash cache depending upon order of
  // LAI.  Good for C order, bad for Fortran
  for (int i=offset; i!=offset+n_grid_cells; ++i) {
    for (int j=0; j!=n_pfts; ++j) {
81
      sai[i][j] = data[(i-offset)*n_pfts + j];
82
83
84
85
86
87
88
89
90
91
92
93
94
    }
  }

  status = nc_close(ncid);
  NC_HANDLE_ERROR( status, "nc_close" ) ;
}


//
// Reads NetCDF met forcing data, returning the min number of times read.
// -----------------------------------------------------------------------------
template<typename Matrix_t>
int read_forcing(const std::string& fname,
95
                 size_t n_times, size_t start_grid_cell, size_t n_grid_cells,
96
97
98
99
100
                 Matrix_t& rain, Matrix_t& snow, Matrix_t& temp) {

  size_t min_ntimes = n_times;
  for (size_t lcv_gc=0; lcv_gc!=n_grid_cells; ++lcv_gc) {
    std::stringstream fname_full;
101
    fname_full << fname << lcv_gc+1+start_grid_cell << ".nc";
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    
    int ncid = -1;
    auto status = nc_open(fname_full.str().c_str(), NC_NOWRITE, &ncid);
    NC_HANDLE_ERROR(status, std::string("nc_open")+" \""+fname_full.str()+"\"");

    int dimid = -1;
    status = nc_inq_dimid(ncid, "time", &dimid);
    NC_HANDLE_ERROR(status, "nc_inq_dimid");

    size_t ntimes_l;
    status = nc_inq_dimlen(ncid, dimid, &ntimes_l);
    NC_HANDLE_ERROR(status, "nc_inq_dimlen");

    // note the number of local times may be less than the max number of times
    assert(ntimes_l <= n_times);
    min_ntimes = std::min(ntimes_l, min_ntimes);
    auto start = std::array<size_t,3>{0,0,0};
    auto count = std::array<size_t,3>{ntimes_l, 1, 1};
    std::vector<double> data_precip(ntimes_l);
    std::vector<double> data_temp(ntimes_l);

    int varid = -1;
    status = nc_inq_varid(ncid,"PRECTmms", &varid);
    NC_HANDLE_ERROR(status, "nc_inq_varid");

    status = nc_get_vara_double(ncid, varid, start.data(), count.data(), data_precip.data());
    NC_HANDLE_ERROR( status, "nc_get_vara_double total_precip" );

    status = nc_inq_varid(ncid,"TBOT", &varid);
    NC_HANDLE_ERROR(status, "nc_inq_varid");

    status = nc_get_vara_double(ncid, varid, start.data(), count.data(), data_temp.data());
    NC_HANDLE_ERROR( status, "nc_get_vara_double temperature" );
    
    status = nc_close(ncid);
    NC_HANDLE_ERROR( status, "nc_close" ) ;

    // allocate the precip to rain or snow
    for (int lcv_t=0; lcv_t!=ntimes_l; ++lcv_t) {
      if (data_temp[lcv_t] < 273.15) {
142
143
144
        snow[lcv_t][lcv_gc] = data_precip[lcv_t];
        rain[lcv_t][lcv_gc] = 0.;
        temp[lcv_t][lcv_gc] = data_temp[lcv_t];
145
      } else {
146
147
148
        snow[lcv_t][lcv_gc] = 0.;
        rain[lcv_t][lcv_gc] = data_precip[lcv_t];
        temp[lcv_t][lcv_gc] = data_temp[lcv_t];
149
150
151
152
153
154
155
156
157
158
159
160
      }
    }
  }
  return min_ntimes;
}


} // namespace Utils
} // namespace ELM


#endif