Vargen.cc 18.1 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

18
19
20
/*
   This module contains the following operators:

21
22
      Vargen     const           Create a constant field
      Vargen     random          Field with random values
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
      Vargen     stdatm          Field values for pressure and temperature for the standard atmosphere
24
25
*/

26
#include <cstdlib>
Oliver Heidmann's avatar
Oliver Heidmann committed
27
#include <cassert>
Ralf Müller's avatar
Ralf Müller committed
28
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
29

30
#include "cdo_options.h"
31
#include "process_int.h"
32
#include "cdo_cdi_wrapper.h"
33
#include "cdo_zaxis.h"
34
#include "param_conversion.h"
35
#include <mpim_grid.h>
36
#include "griddes.h"
37
#include "constants.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
#include "stdnametable.h"
39
#include "timebase.h"
40
#include "param_conversion.h"
41

42
43
static constexpr double etopo_scale = 3.0;
static constexpr double etopo_offset = 11000.0;
44
static constexpr unsigned short etopo[] = {
45
#include "etopo.dat"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47

48
49
static constexpr double temp_scale = 500.0;
static constexpr double temp_offset = -220.0;
50
static constexpr unsigned short temp[] = {
51
#include "temp.dat"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53

54
55
static constexpr double mask_scale = 1.0;
static constexpr double mask_offset = 0.0;
56
static constexpr unsigned short mask[] = {
57
#include "mask.dat"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59

Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
// Some Constants for creating temperatur and pressure for the standard atmosphere
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
62
63
64
65
constexpr double T_ZERO = 213.0;
constexpr double T_DELTA = 75.0;
constexpr double SCALEHEIGHT = 10000.0;  // [m]
constexpr double P_ZERO = 1013.25;       // surface pressure [hPa]
constexpr double CC_R = 287.05;          // specific gas constant for air
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
67
constexpr double TMP4PRESSURE = (C_EARTH_GRAV * SCALEHEIGHT) / (CC_R * T_ZERO);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
static double
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
std_atm_temperatur(const double height)
70
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
  // Compute the temperatur for the given height (in meters) according to the solution of the hydrostatic atmosphere
72
  return (T_ZERO + T_DELTA * std::exp((-1) * (height / SCALEHEIGHT)));
73
74
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
static double
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
std_atm_pressure(const double height)
77
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
  // Compute the pressure for the given height (in meters) according to the solution of the hydrostatic atmosphere
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
  return (P_ZERO
          * std::exp((-1) * TMP4PRESSURE * std::log((std::exp(height / SCALEHEIGHT) * T_ZERO + T_DELTA) / (T_ZERO + T_DELTA))));
81
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
static void
84
conv_generic_grid(const int gridID, const size_t gridsize, Varray<double> &xvals2D, Varray<double> &yvals2D)
85
{
86
87
  const auto xsize = gridInqXsize(gridID);
  const auto ysize = gridInqYsize(gridID);
88

89
  assert(gridsize == xsize * ysize);
90

Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
  Varray<double> xcoord(xsize), ycoord(ysize);
92
93
  gridInqXvals(gridID, &xcoord[0]);
  gridInqYvals(gridID, &ycoord[0]);
94

95
96
  const auto xrange = varrayRange(xsize, xcoord);
  const auto yrange = varrayRange(ysize, ycoord);
97
98
99

  for (size_t j = 0; j < ysize; ++j)
    for (size_t i = 0; i < xsize; ++i)
100
      {
101
102
        xvals2D[j * xsize + i] = xcoord[i] * M_PI / xrange;
        yvals2D[j * xsize + i] = ycoord[j] * M_PI / yrange;
103
104
105
      }
}

106
static void
107
remap_nn_reg2d_reg2d(size_t nx, size_t ny, const double *restrict data, const int gridID, double *restrict array)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
  if (gridInqType(gridID) != GRID_LONLAT) cdoAbort("Internal error, wrong grid type!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110

111
112
  auto nxvals = gridInqXsize(gridID);
  auto nyvals = gridInqYsize(gridID);
113
  Varray<double> xvals(nxvals), yvals(nyvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114

115
116
  gridInqXvals(gridID, &xvals[0]);
  gridInqYvals(gridID, &yvals[0]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117

118
119
120
  // Convert lat/lon units if required
  cdo_grid_to_degree(gridID, CDI_XAXIS, nxvals, &xvals[0], "grid center lon");
  cdo_grid_to_degree(gridID, CDI_YAXIS, nyvals, &yvals[0], "grid center lat");
121

122
#ifdef _OPENMP
123
#pragma omp parallel for default(none) shared(nx, ny, nxvals, nyvals, xvals, yvals, data, array)
124
#endif
125
  for (size_t j = 0; j < nyvals; j++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
    {
127
      const auto yval = yvals[j];
128
      for (size_t i = 0; i < nxvals; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
        {
130
          auto xval = xvals[i];
131
132
133
134
135
136
137
          if (xval >= 180) xval -= 360;
          if (xval < -180) xval += 360;
          size_t ii = (xval + 180) * 2;
          size_t jj = (yval + 90) * 2;
          if (ii >= nx) ii = nx - 1;
          if (jj >= ny) jj = ny - 1;
          array[j * nxvals + i] = data[jj * nx + ii];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
140
141
        }
    }
}

142
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
remap_nn_reg2d_nonreg2d(const size_t nx, const size_t ny, const double *restrict data, const int gridID, double *restrict array)
144
{
145
  auto gridID2 = gridID;
146
  const auto gridsize = gridInqSize(gridID2);
147
  Varray<double> xvals(gridsize), yvals(gridsize);
148

Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
  if (gridInqType(gridID2) == GRID_GME || gridInqType(gridID2) == GRID_GAUSSIAN_REDUCED) gridID2 = gridToUnstructured(gridID2, 0);
150

Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
  if (gridInqType(gridID2) != GRID_UNSTRUCTURED && gridInqType(gridID2) != GRID_CURVILINEAR)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
    gridID2 = gridToCurvilinear(gridID2, 0);
153

154
155
  gridInqXvals(gridID2, &xvals[0]);
  gridInqYvals(gridID2, &yvals[0]);
156

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
  // Convert lat/lon units if required
158
159
  cdo_grid_to_degree(gridID2, CDI_XAXIS, gridsize, &xvals[0], "grid center lon");
  cdo_grid_to_degree(gridID2, CDI_YAXIS, gridsize, &yvals[0], "grid center lat");
160

161
  for (size_t i = 0; i < gridsize; i++)
162
    {
163
164
      auto xval = xvals[i];
      auto yval = yvals[i];
165
166
167
168
169
170
171
      if (xval >= 180) xval -= 360;
      if (xval < -180) xval += 360;
      size_t ii = (xval + 180) * 2;
      size_t jj = (yval + 90) * 2;
      if (ii >= nx) ii = nx - 1;
      if (jj >= ny) jj = ny - 1;
      array[i] = data[jj * nx + ii];
172
173
    }

174
  if (gridID != gridID2) gridDestroy(gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
176
}

177
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
remap_nn_reg2d(const size_t nx, const size_t ny, const double *restrict data, const int gridID, double *restrict array)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
{
180
  if (gridInqType(gridID) == GRID_LONLAT)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
183
    remap_nn_reg2d_reg2d(nx, ny, data, gridID, array);
  else
    remap_nn_reg2d_nonreg2d(nx, ny, data, gridID, array);
184
185
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
static int
randomInit(const int operatorID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
  unsigned int seed = Options::Random_Seed;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
192
  operatorInputArg(cdoOperatorEnter(operatorID));
  if (operatorArgc() < 1) cdoAbort("Too few arguments!");
  if (operatorArgc() > 2) cdoAbort("Too many arguments!");
193
  const auto gridID = cdoDefineGrid(cdoOperatorArgv(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
  if (operatorArgc() == 2)
    {
196
      const auto idum = parameter2int(cdoOperatorArgv(1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
198
      if (idum >= 0 && idum < 0x7FFFFFFF) seed = idum;
    }
199
  std::srand(seed);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
201
202
  return gridID;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
static void
204
randomCompute(const size_t gridsize, Varray<double> &array)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
{
206
  for (size_t i = 0; i < gridsize; ++i) array[i] = ((double) std::rand()) / ((double) RAND_MAX);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
static void
210
sincosCompute(const size_t gridsize, Varray<double> &array, const Varray<double> &xvals, const Varray<double> &yvals)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
{
212
  for (size_t i = 0; i < gridsize; ++i) array[i] = std::cos(1.0 * xvals[i]) * std::sin(2.0 * yvals[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
coshillCompute(const size_t gridsize, Varray<double> &array, const Varray<double> &xvals, const Varray<double> &yvals)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
{
218
  for (size_t i = 0; i < gridsize; ++i) array[i] = 2 - std::cos(std::acos(std::cos(xvals[i]) * std::cos(yvals[i])) / 1.2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
testfieldCompute(const size_t gridsize, Varray<double> &array, const Varray<double> &xvals, const Varray<double> &yvals)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224
225
226
227
{
  double xyz[3];
  for (size_t i = 0; i < gridsize; ++i)
    {
      gcLLtoXYZ(xvals[i], yvals[i], xyz);
228
229
230
      const auto x = xyz[0];
      const auto y = xyz[1];
      const auto z = xyz[2];
231
      array[i] = 1.0 + std::pow(x, 8.0) + std::exp(2.0 * y * y * y) + std::exp(2.0 * x * x) + 10.0 * x * y * z;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
233
234
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
static void
236
unpackData(const size_t datasize, Varray<double> &data, const double scale, const double offset, const unsigned short *zdata)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
238
239
240
{
  for (size_t i = 0; i < datasize; ++i) data[i] = zdata[i] / scale - offset;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
static int
definePointGrid()
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
  const auto gridID = gridCreate(GRID_LONLAT, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
246
  gridDefXsize(gridID, 1);
  gridDefYsize(gridID, 1);
247
  const double value = 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
252
253
  gridDefXvals(gridID, &value);
  gridDefYvals(gridID, &value);

  return gridID;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
255
static int
defineZaxis(const bool lstdatm, const int nlevels, double *levels)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
256
257
258
259
260
261
262
{
  int zaxisID = -1;

  if (lstdatm)
    {
      zaxisID = zaxisCreate(ZAXIS_HEIGHT, nlevels);
      zaxisDefLevels(zaxisID, levels);
263
      cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_NAME, "level");
264
265
      cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_LONGNAME, "Level");
      cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_UNITS, "m");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
268
    }
  else
    {
269
      zaxisID = zaxisFromName("surface");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
271
272
273
274
    }

  return zaxisID;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
276
static void
definePressureAttributes(const int vlistID, const int varID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277
{
278
  cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, "P");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
  vlistDefVarParam(vlistID, varID, cdiEncodeParam(1, 255, 255));
280
281
282
  cdiDefKeyString(vlistID, varID, CDI_KEY_STDNAME, "air_pressure");
  cdiDefKeyString(vlistID, varID, CDI_KEY_LONGNAME, "pressure");
  cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, "hPa");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
284
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
285
286
static void
defineTemperatureAttributes(const int vlistID, const int varID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
287
{
288
  cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, "T");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
  vlistDefVarParam(vlistID, varID, cdiEncodeParam(130, 128, 255));
290
291
292
  cdiDefKeyString(vlistID, varID, CDI_KEY_STDNAME, var_stdname(air_temperature));
  cdiDefKeyString(vlistID, varID, CDI_KEY_LONGNAME, "temperature");
  cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, "K");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293
}
294

295
296
void *
Vargen(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
299
300
  constexpr size_t nlat = 360;
  constexpr size_t nlon = 720;
  double lon[nlon], lat[nlat];
301
  int ntimesteps, nlevels = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
  int gridID = -1, gridIDdata = -1;
303
304
  double rstart = 0.0, rstop = 0.0, rinc = 0.0;
  double rconst = 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
305
  std::vector<double> levels;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306

307
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308

Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
  // clang-format off
310
311
312
  const auto RANDOM    = cdoOperatorAdd("random",    0, 0, "grid description file or name, <seed>");
  const auto SINCOS    = cdoOperatorAdd("sincos",    0, 0, "grid description file or name");
  const auto COSHILL   = cdoOperatorAdd("coshill",   0, 0, "grid description file or name");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
  const auto TESTFIELD = cdoOperatorAdd("testfield", 0, 0, "grid description file or name");
314
315
316
317
318
319
  const auto CONST     = cdoOperatorAdd("const",     0, 0, "constant value, grid description file or name");
  const auto SEQ       = cdoOperatorAdd("seq",       0, 0, "start, end, <increment>");
  const auto TOPO      = cdoOperatorAdd("topo",      0, 0, nullptr);
  const auto TEMP      = cdoOperatorAdd("temp",      0, 0, nullptr);
  const auto MASK      = cdoOperatorAdd("mask",      0, 0, nullptr);
  const auto STDATM    = cdoOperatorAdd("stdatm",    0, 0, "height levels [m]");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321

322
  const auto operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323

324
  if (operatorID == RANDOM)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
326
      gridID = randomInit(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
  else if (operatorID == SINCOS || operatorID == COSHILL || operatorID == TESTFIELD)
329
330
331
    {
      operatorInputArg(cdoOperatorEnter(operatorID));
      operatorCheckArgc(1);
332
      gridID = cdoDefineGrid(cdoOperatorArgv(0));
333
    }
334
  else if (operatorID == CONST)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
      operatorInputArg(cdoOperatorEnter(operatorID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
      operatorCheckArgc(2);
338
339
      rconst = parameter2double(cdoOperatorArgv(0));
      gridID = cdoDefineGrid(cdoOperatorArgv(1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
    }
341
  else if (operatorID == TOPO || operatorID == TEMP || operatorID == MASK)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
    {
343
      gridIDdata = gridCreate(GRID_LONLAT, nlon * nlat);
344
345
346
      gridDefXsize(gridIDdata, nlon);
      gridDefYsize(gridIDdata, nlat);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
348
      for (size_t i = 0; i < nlon; i++) lon[i] = -179.75 + i * 0.5;
      for (size_t i = 0; i < nlat; i++) lat[i] = -89.75 + i * 0.5;
349
350
351
352
353
354

      gridDefXvals(gridIDdata, lon);
      gridDefYvals(gridIDdata, lat);

      gridID = gridIDdata;

355
      if (operatorArgc() == 1) gridID = cdoDefineGrid(cdoOperatorArgv(0));
356
      if (operatorArgc() > 1) cdoAbort("Too many arguments!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
357
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
358
  else if (operatorID == SEQ)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
360
    {
      operatorInputArg(cdoOperatorEnter(operatorID));
361
362
      if (operatorArgc() < 2) cdoAbort("Too few arguments!");
      if (operatorArgc() > 3) cdoAbort("Too many arguments!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363

364
365
366
      rstart = parameter2double(cdoOperatorArgv(0));
      rstop = parameter2double(cdoOperatorArgv(1));
      rinc = (operatorArgc() == 3) ? parameter2double(cdoOperatorArgv(2)) : 1;
367
      if (DBL_IS_EQUAL(rinc, 0.0)) cdoAbort("Increment is zero!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368

Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
      gridID = definePointGrid();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
    }
371
  else if (operatorID == STDATM)
372
    {
373
      operatorInputArg(cdoOperatorEnter(operatorID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
375
      levels = cdoArgvToFlt(cdoGetOperArgv());
      nlevels = levels.size();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
376

377
      if (Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378
        for (int i = 0; i < nlevels; ++i) printf("levels %d: %g\n", i, levels[i]);
379

Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
      gridID = definePointGrid();
381
382
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
  const auto zaxisID = defineZaxis(operatorID == STDATM, nlevels, levels.data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384

385
  const auto vlistID = vlistCreate();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386

387
  const auto timetype = (operatorID == SEQ) ? TIME_VARYING : TIME_CONSTANT;
388

389
  auto varID = vlistDefVar(vlistID, gridID, zaxisID, timetype);
390
  /*
391
392
393
     For the standard atmosphere two output variables are generated: pressure
     and temperatur. The first (varID) is pressure, second (varID2) is
     temperatur. Add an additional variable for the standard atmosphere.
394
   */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
  const int varID2 = (operatorID == STDATM) ? vlistDefVar(vlistID, gridID, zaxisID, TIME_CONSTANT) : -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396

Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
  if (operatorID == MASK) vlistDefVarDatatype(vlistID, varID, CDI_DATATYPE_INT8);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398

399
  if (operatorID == STDATM)
400
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
401
402
      definePressureAttributes(vlistID, varID);
      defineTemperatureAttributes(vlistID, varID2);
403
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404
405
  else
    {
406
407
408
      cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, cdoOperatorName(operatorID));
      if (operatorID == TOPO) cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, "m");
      if (operatorID == TEMP) cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, "K");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410

Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
  const auto taxisID = cdoTaxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
413
  vlistDefTaxis(vlistID, taxisID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
415
  if (operatorID == RANDOM || operatorID == SINCOS || operatorID == COSHILL || operatorID == TESTFIELD || operatorID == CONST
      || operatorID == TOPO || operatorID == TEMP || operatorID == MASK || operatorID == STDATM)
416
417
    vlistDefNtsteps(vlistID, 1);

418
  const auto streamID = cdoOpenWrite(0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419

420
  cdoDefVlist(streamID, vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421

422
  const auto gridsize = gridInqSize(gridID);
423
  Varray<double> array(gridsize);
424

Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
  if (operatorID == SEQ)
426
    ntimesteps = 1.001 + ((rstop - rstart) / rinc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
  else
428
429
430
431
    {
      vlistDefNtsteps(vlistID, 0);
      ntimesteps = 1;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
432

433
  const auto julday = date_to_julday(CALENDAR_PROLEPTIC, 10101);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434

435
  const auto nvars = vlistNvars(vlistID);
436

437
  for (int tsID = 0; tsID < ntimesteps; tsID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
      const double rval = rstart + rinc * tsID;
440
      const auto vdate = julday_to_date(CALENDAR_PROLEPTIC, julday + tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
      const int vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
443
      taxisDefVdate(taxisID, vdate);
      taxisDefVtime(taxisID, vtime);
444
      cdoDefTimestep(streamID, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
445

446
      for (varID = 0; varID < nvars; varID++)
447
448
        {
          nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449
          for (int levelID = 0; levelID < nlevels; levelID++)
450
            {
451
              cdoDefRecord(streamID, varID, levelID);
452

453
              if (operatorID == RANDOM)
454
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
                  randomCompute(gridsize, array);
456
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
              else if (operatorID == SINCOS || operatorID == COSHILL || operatorID == TESTFIELD)
458
                {
459
                  Varray<double> xvals(gridsize), yvals(gridsize);
460

461
                  if (grid_is_distance_generic(gridID))
462
463
464
465
466
                    {
                      conv_generic_grid(gridID, gridsize, xvals, yvals);
                    }
                  else
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
                      if (gridInqType(gridID) == GRID_GME) gridID = gridToUnstructured(gridID, 0);
468

Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
                      if (gridInqType(gridID) != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_CURVILINEAR)
470
471
                        gridID = gridToCurvilinear(gridID, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
473
                      gridInqXvals(gridID, xvals.data());
                      gridInqYvals(gridID, yvals.data());
474

Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
                      // Convert lat/lon units if required
476
477
                      cdo_grid_to_radian(gridID, CDI_XAXIS, gridsize, xvals.data(), "grid center lon");
                      cdo_grid_to_radian(gridID, CDI_YAXIS, gridsize, yvals.data(), "grid center lat");
478
                    }
479
480
481

                  if (operatorID == SINCOS)
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
                      sincosCompute(gridsize, array, xvals, yvals);
483
484
485
                    }
                  else if (operatorID == COSHILL)
                    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
                      coshillCompute(gridsize, array, xvals, yvals);
487
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488
489
490
491
                  else if (operatorID == TESTFIELD)
                    {
                      testfieldCompute(gridsize, array, xvals, yvals);
                    }
492
493
                }
              else if (operatorID == CONST)
494
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
                  for (size_t i = 0; i < gridsize; i++) array[i] = rconst;
496
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
              else if (operatorID == TOPO || operatorID == TEMP || operatorID == MASK)
498
                {
499
                  const auto datasize = gridInqSize(gridIDdata);
500
                  Varray<double> data(datasize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515

                  // clang-format off
                  if      (operatorID == TOPO) unpackData(datasize, data, etopo_scale, etopo_offset, etopo);
                  else if (operatorID == TEMP) unpackData(datasize, data, temp_scale, temp_offset, temp);
                  else if (operatorID == MASK) unpackData(datasize, data, mask_scale, mask_offset, mask);
                  // clang-format on

                  if (gridID != gridIDdata && gridIDdata != -1)
                    {
                      remap_nn_reg2d(nlon, nlat, data.data(), gridID, array.data());
                    }
                  else
                    {
                      for (size_t i = 0; i < gridsize; ++i) array[i] = data[i];
                    }
516
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
              else if (operatorID == SEQ)
518
519
520
                {
                  array[0] = rval;
                }
521
              else if (operatorID == STDATM)
522
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
                  array[0] = (varID == varID2) ? std_atm_temperatur(levels[levelID]) : std_atm_pressure(levels[levelID]);
524
525
                }

526
              cdoWriteRecord(streamID, array.data(), 0);
527
528
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
530
    }

531
  cdoStreamClose(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532

533
534
  vlistDestroy(vlistID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
535
536
  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
}