Wind.cc 13.6 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
21
22
/*
   This module contains the following operators:

      Wind       uv2dv           U and V wind to divergence and vorticity
      Wind       dv2uv           Divergence and vorticity to U and V wind
23
      Wind       dv2ps           Divergence and vorticity to velocity potential and stream function
24
25
*/

Ralf Mueller's avatar
Ralf Mueller committed
26
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
27

28
#include "cdo_vlist.h"
29
#include "process_int.h"
30
#include "param_conversion.h"
31
#include <mpim_grid.h>
32
#include "griddes.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
#include "specspace.h"
34
#include "util_string.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35

36
37
38
39
40
41
42
43

static void
defineAttributesDV(int vlistID2, int gridID2, int varID1, int varID2)
{
  vlistChangeVarGrid(vlistID2, varID1, gridID2);
  vlistChangeVarGrid(vlistID2, varID2, gridID2);
  vlistDefVarParam(vlistID2, varID1, cdiEncodeParam(155, 128, 255));
  vlistDefVarParam(vlistID2, varID2, cdiEncodeParam(138, 128, 255));
44
45
46
47
48
49
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_NAME, "sd");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_NAME, "svo");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_LONGNAME, "divergence");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_LONGNAME, "vorticity");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_UNITS, "1/s");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_UNITS, "1/s");
50
51
52
53
54
55
56
57
58
}

static void
defineAttributesUV(int vlistID2, int gridID2, int varID1, int varID2)
{
  vlistChangeVarGrid(vlistID2, varID1, gridID2);
  vlistChangeVarGrid(vlistID2, varID2, gridID2);
  vlistDefVarParam(vlistID2, varID1, cdiEncodeParam(131, 128, 255));
  vlistDefVarParam(vlistID2, varID2, cdiEncodeParam(132, 128, 255));
59
60
61
62
63
64
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_NAME, "u");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_NAME, "v");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_LONGNAME, "u-velocity");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_LONGNAME, "v-velocity");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_UNITS, "m/s");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_UNITS, "m/s");
65
66
67
68
69
70
71
}

static void
defineAttributesPS(int vlistID2, int varID1, int varID2)
{
  vlistDefVarParam(vlistID2, varID1, cdiEncodeParam(149, 128, 255));
  vlistDefVarParam(vlistID2, varID2, cdiEncodeParam(148, 128, 255));
72
73
74
75
76
77
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_NAME, "velopot");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_NAME, "stream");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_LONGNAME, "velocity potential");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_LONGNAME, "streamfunction");
  cdiDefKeyString(vlistID2, varID1, CDI_KEY_UNITS, "m^2/s");
  cdiDefKeyString(vlistID2, varID2, CDI_KEY_UNITS, "m^2/s");
78
79
}

80
81
void *
Wind(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
{
83
  int nrecs;
84
  size_t nlev = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
  int gridID1 = -1, gridID2 = -1;
86
  size_t nmiss;
87
  long ntr = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
  int varID1 = -1, varID2 = -1;
89
90
  SP_Transformation spTrans;
  DV_Transformation dvTrans;
91
  char varname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

93
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94

95
  const auto lcopy = unchangedRecord();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96

Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
  // clang-format off
98
99
100
101
102
  const auto UV2DV  = cdoOperatorAdd("uv2dv",  0, 0, nullptr);
  const auto UV2DVL = cdoOperatorAdd("uv2dvl", 0, 0, nullptr);
  const auto DV2UV  = cdoOperatorAdd("dv2uv",  0, 0, nullptr);
  const auto DV2UVL = cdoOperatorAdd("dv2uvl", 0, 0, nullptr);
  const auto DV2PS  = cdoOperatorAdd("dv2ps",  0, 0, nullptr);
103

104
  const auto operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105

106
107
108
  const bool luv2dv = (operatorID == UV2DV || operatorID == UV2DVL);
  const bool ldv2uv = (operatorID == DV2UV || operatorID == DV2UVL);
  const bool linear = (operatorID == UV2DVL || operatorID == DV2UVL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109

110
111
  int (*nlat2ntr)(int) = linear ? nlat_to_ntr_linear : nlat_to_ntr;
  const char *ctype = linear ? "l" : "";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
114
  if ((luv2dv || ldv2uv) && operatorArgc() == 1)
    {
115
      std::string type = parameter2word(cdoOperatorArgv(0));
116
117
118
119
      if      (type == "linear")    { nlat2ntr = nlat_to_ntr_linear; ctype = "l"; }
      else if (type == "cubic")     { nlat2ntr = nlat_to_ntr_cubic; ctype = "c"; }
      else if (type == "quadratic") { nlat2ntr = nlat_to_ntr; }
      else cdoAbort("Unsupported type: %s\n", type.c_str());
120
121
    }
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122

123
  const auto streamID1 = cdoOpenRead(0);
124

125
126
  const auto vlistID1 = cdoStreamInqVlist(streamID1);
  const auto vlistID2 = vlistDuplicate(vlistID1);
127

128
129
  const auto taxisID1 = vlistInqTaxis(vlistID1);
  const auto taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
  vlistDefTaxis(vlistID2, taxisID2);

132
  // find variables
133
  const auto nvars = vlistNvars(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
  for (int varID = 0; varID < nvars; varID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
    {
136
      const auto param = vlistInqVarParam(vlistID2, varID);
137
      int pnum, pcat, pdis;
138
      cdiDecodeParam(param, &pnum, &pcat, &pdis);
139
      int code = pnum;
140
141
      if (operatorID == UV2DV || operatorID == UV2DVL)
        {
142
          // search for u and v wind
143
144
145
          if (pdis != 255 || code <= 0)
            {
              vlistInqVarName(vlistID1, varID, varname);
146
              cstrToLowerCase(varname);
147
148
              if (cstrIsEqual(varname, "u")) code = 131;
              if (cstrIsEqual(varname, "v")) code = 132;
149
150
            }

151
152
          if (code == 131) varID1 = varID;
          if (code == 132) varID2 = varID;
153
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
      else if (operatorID == DV2UV || operatorID == DV2UVL || operatorID == DV2PS)
155
        {
156
          // search for divergence and vorticity
157
158
159
          if (pdis != 255)  // GRIB2
            {
              vlistInqVarName(vlistID1, varID, varname);
160
              cstrToLowerCase(varname);
161
162
              if (cstrIsEqual(varname, "d")) code = 155;
              if (cstrIsEqual(varname, "vo")) code = 138;
163
164
165
166
            }
          else if (code <= 0)
            {
              vlistInqVarName(vlistID1, varID, varname);
167
              cstrToLowerCase(varname);
168
169
              if (cstrIsEqual(varname, "sd")) code = 155;
              if (cstrIsEqual(varname, "svo")) code = 138;
170
171
            }

172
173
          if (code == 155) varID1 = varID;
          if (code == 138) varID2 = varID;
174
        }
175
      else
176
        cdoAbort("Unexpected operatorID %d", operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
178
    }

179
180
  auto gridIDsp = vlistGetFirstSpectralGrid(vlistID1);
  auto gridIDgp = vlistGetFirstGaussianGrid(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
183
  // define output grid
  if (luv2dv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
    {
185
186
      if (varID1 == -1) cdoWarning("U-wind not found!");
      if (varID2 == -1) cdoWarning("V-wind not found!");
187
188
189
190
191

      if (varID1 != -1 && varID2 != -1)
        {
          gridID1 = vlistInqVarGrid(vlistID1, varID1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
          if (gridInqType(gridID1) != GRID_GAUSSIAN) cdoAbort("U-wind is not on Gaussian grid!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
          if (gridID1 != vlistInqVarGrid(vlistID1, varID2)) cdoAbort("U and V wind must have the same grid represention!");
194

195
196
197
198
          const auto numLPE = gridInqNP(gridID1);
          const long nlon = gridInqXsize(gridID1);
          const long nlat = gridInqYsize(gridID1);
          const long ntr1 = nlat2ntr(nlat);
199

200
          if (nlat != (numLPE * 2)) cdoAbort("U and V fields on Gaussian grid are not global!");
201

202
203
          if (gridIDsp != -1)
            if (ntr1 != gridInqTrunc(gridIDsp)) gridIDsp = -1;
204

205
206
207
208
209
210
          if (gridIDsp == -1)
            {
              gridIDsp = gridCreate(GRID_SPECTRAL, (ntr1 + 1) * (ntr1 + 2));
              gridDefTrunc(gridIDsp, ntr1);
              gridDefComplexPacking(gridIDsp, 1);
            }
211
212
213
214
215

          if (gridIDsp == -1) cdoAbort("No Gaussian grid data found!");

          gridID2 = gridIDsp;

216
          defineAttributesDV(vlistID2, gridID2, varID1, varID2);
217
218

          ntr = gridInqTrunc(gridID2);
219
          nlev = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID1));
220

221
          spTrans.init(nlon, nlat, ntr, 1, nlev);
222
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
    }
224
  else if (ldv2uv)
225
    {
226
227
      if (varID1 == -1) cdoWarning("Divergence not found!");
      if (varID2 == -1) cdoWarning("Vorticity not found!");
228
229
230
231
232

      if (varID1 != -1 && varID2 != -1)
        {
          gridID1 = vlistInqVarGrid(vlistID1, varID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
          if (gridInqType(gridID1) != GRID_SPECTRAL) cdoAbort("Vorticity is not on spectral grid!");
234
235

          if (gridID1 != vlistInqVarGrid(vlistID1, varID1))
236
            cdoAbort("Divergence and vorticity must have the same grid represention!");
237
238
239

          if (gridIDgp != -1)
            {
240
              const long nlat = gridInqYsize(gridIDgp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
              const long ntr1 = nlat2ntr(nlat);
              if (gridInqTrunc(gridIDsp) != ntr1) gridIDgp = -1;
243
244
245
246
247
            }

          if (gridIDgp == -1)
            {
              char gridname[20];
248
              snprintf(gridname, sizeof(gridname), "t%s%dgrid", ctype, gridInqTrunc(gridIDsp));
249
              gridIDgp = gridFromName(gridname);
250
251
252
253
            }

          gridID2 = gridIDgp;

254
          defineAttributesUV(vlistID2, gridID2, varID1, varID2);
255

256
257
          const long nlon = gridInqXsize(gridID2);
          const long nlat = gridInqYsize(gridID2);
258
          ntr = gridInqTrunc(gridID1);
259
          nlev = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID1));
260

261
          spTrans.init(nlon, nlat, ntr, 0, nlev);
262
          dvTrans.init(ntr);
263
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
    }
265
266
  else if (operatorID == DV2PS)
    {
267
268
      if (varID1 == -1) cdoWarning("Divergence not found!");
      if (varID2 == -1) cdoWarning("Vorticity not found!");
269
270
271
272
273

      if (varID1 != -1 && varID2 != -1)
        {
          gridID1 = vlistInqVarGrid(vlistID1, varID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
274
          if (gridInqType(gridID1) != GRID_SPECTRAL) cdoAbort("Vorticity is not on spectral grid!");
275
276

          if (gridID1 != vlistInqVarGrid(vlistID1, varID1))
277
            cdoAbort("Divergence and vorticity must have the same grid represention!");
278

279
          defineAttributesPS(vlistID2, varID1, varID2);
280
281
282
283

          ntr = gridInqTrunc(gridID1);
          gridID2 = gridID1;
        }
284
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285

286
  const auto streamID2 = cdoOpenWrite(1);
287
  cdoDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288

289
  const auto gridsizemax = vlistGridsizeMax(vlistID1);
290
  Varray<double> array1(gridsizemax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
291

292
  Varray<double> ivar1, ivar2, ovar1, ovar2;
293
  if (varID1 != -1 && varID2 != -1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
    {
295
      nlev = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296

297
      auto gridsize = gridInqSize(gridID1);
298
299
      ivar1.resize(nlev * gridsize);
      ivar2.resize(nlev * gridsize);
300

Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
      gridsize = gridInqSize(gridID2);
302
303
      ovar1.resize(nlev * gridsize);
      ovar2.resize(nlev * gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
305

306
  int tsID = 0;
307
  while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
309
    {
      taxisCopyTimestep(taxisID2, taxisID1);
310
      cdoDefTimestep(streamID2, tsID);
311
312
313

      for (int recID = 0; recID < nrecs; recID++)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314
          int varID, levelID;
315
          cdoInqRecord(streamID1, &varID, &levelID);
316

Uwe Schulzweida's avatar
Uwe Schulzweida committed
317
          if ((varID1 != -1 && varID2 != -1) && (varID == varID1 || varID == varID2))
318
            {
319
              cdoReadRecord(streamID1, array1.data(), &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320
              if (nmiss) cdoAbort("Missing values unsupported for spectral data!");
321

322
323
              const auto gridsize = gridInqSize(gridID1);
              const auto offset = gridsize * levelID;
324
325

              if (varID == varID1)
326
                arrayCopy(gridsize, array1.data(), &ivar1[offset]);
327
              else if (varID == varID2)
328
                arrayCopy(gridsize, array1.data(), &ivar2[offset]);
329
330
331
            }
          else
            {
332
              cdoDefRecord(streamID2, varID, levelID);
333
334
              if (lcopy)
                {
335
                  cdoCopyRecord(streamID2, streamID1);
336
337
338
                }
              else
                {
339
340
                  cdoReadRecord(streamID1, array1.data(), &nmiss);
                  cdoWriteRecord(streamID2, array1.data(), nmiss);
341
342
343
344
345
346
                }
            }
        }

      if (varID1 != -1 && varID2 != -1)
        {
347
          if (luv2dv)
348
            trans_uv2dv(spTrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
349
          else if (ldv2uv)
350
            trans_dv2uv(spTrans, dvTrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
351
352
          else if (operatorID == DV2PS)
            {
353
354
              dv2ps(ivar1.data(), ovar1.data(), nlev, ntr);
              dv2ps(ivar2.data(), ovar2.data(), nlev, ntr);
355
356
            }

357
          const auto gridsize = gridInqSize(gridID2);
358
          if (luv2dv || operatorID == DV2PS)
359
            {
360
              for (size_t levelID = 0; levelID < nlev; levelID++)
361
                {
362
                  const auto offset = gridsize * levelID;
363
364
                  cdoDefRecord(streamID2, varID2, levelID);
                  cdoWriteRecord(streamID2, &ovar2[offset], 0);
365
                }
366
              for (size_t levelID = 0; levelID < nlev; levelID++)
367
                {
368
                  const auto offset = gridsize * levelID;
369
370
                  cdoDefRecord(streamID2, varID1, levelID);
                  cdoWriteRecord(streamID2, &ovar1[offset], 0);
371
372
                }
            }
373
          else if (ldv2uv)
374
            {
375
              for (size_t levelID = 0; levelID < nlev; levelID++)
376
                {
377
                  const auto offset = gridsize * levelID;
378
379
                  cdoDefRecord(streamID2, varID1, levelID);
                  cdoWriteRecord(streamID2, &ovar1[offset], 0);
380
                }
381
              for (size_t levelID = 0; levelID < nlev; levelID++)
382
                {
383
                  const auto offset = gridsize * levelID;
384
385
                  cdoDefRecord(streamID2, varID2, levelID);
                  cdoWriteRecord(streamID2, &ovar2[offset], 0);
386
387
388
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389

Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
391
392
      tsID++;
    }

393
394
  cdoStreamClose(streamID2);
  cdoStreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
397

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
399
}