afterburnerlib.cc 94.1 KB
Newer Older
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>
6
7
8
9
10
11
12
13
14
15
16
  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.
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
17

18
19
#include <algorithm>

Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
21
#include "cdi.h"

22
#include "cimdOmp.h"
23
24
#define streamDefRecord cdoDefRecord
#define streamWriteRecord cdoWriteRecord
25

Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
27
28
#include "afterburner.h"
#include "constants.h"
#include "compare.h"
29
#include "vertical_interp.h"
30

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
int afterDebug = 0;

33
34
static char *
FieldName(int code, const char *text)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
37
{
  static char name[256];

38
  snprintf(name, sizeof(name), "[%3d].%s", code, text);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39

40
  return name;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
43
}

/* Free array space */
44
45
static void *
FreeMemory(void *ptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
{
47
  Free(ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48

Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
}

52
53
static void
FreeSpectral(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
{
55
  for (int code = MaxCodes - 1; code >= 0; --code)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
    if (vars[code].spectral) vars[code].spectral = (double *) FreeMemory(vars[code].spectral);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
58
}

59
60
static void
FreeFourier(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
{
62
  for (int code = 0; code < MaxCodes; code++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
    if (vars[code].fourier) vars[code].fourier = (double *) FreeMemory(vars[code].fourier);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
}

66
67
static void
FreeHybrid(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
{
69
  for (int code = 0; code < MaxCodes; code++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
    if (vars[code].hybrid) vars[code].hybrid = (double *) FreeMemory(vars[code].hybrid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
72
}

73
74
static void
FreeGrid(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
{
76
  for (int code = 0; code < MaxCodes; code++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
    if (vars[code].grid) vars[code].grid = (double *) FreeMemory(vars[code].grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
}

80
81
static void
FreeSamp(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
{
83
84
  for (int code = 0; code < MaxCodes; code++)
    if (vars[code].samp) vars[code].samp = (int *) FreeMemory(vars[code].samp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
86
87
}

/* alloc_dp - Allocate space for double array */
88
89
static double *
alloc_dp(int words, const char *array_name)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
  double *result = nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

93
  if (words > 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
    {
95
      result = (double *) Malloc(words * sizeof(double));
96
      if (result == nullptr) cdoSysError(array_name, "No Memory!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
    }

99
  return result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
}

102
103
static void
IniQuaSum(double *dest, const double *restrict src, int len)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
  for (int i = 0; i < len; i++) dest[i] = src[i] * src[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
107
}

108
109
static void
AddQuaSum(double *dest, const double *restrict src, int len)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
  for (int i = 0; i < len; i++) dest[i] += src[i] * src[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
113
}

114
115
static void
VarQuaSum(double *Variance, const double *restrict Sum, int len, int n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
{
117
  const double rn1 = 1.0 / (n - 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
  for (int i = 0; i < len; i++) Variance[i] = (Variance[i] - Sum[i] * Sum[i] * n) * rn1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120

121
  for (int i = 0; i < len; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
    Variance[i] = (Variance[i] > 0.0) ? std::sqrt(Variance[i]) : 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
124
}

125
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
AddVector(double *dest, const double *src, size_t len, size_t *nmiss, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
  if (*nmiss)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
    {
130
131
      arrayAddArrayMV(len, dest, src, missval);
      *nmiss = arrayNumMV(len, dest, missval);
132
      if (*nmiss == 0) *nmiss = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
134
135
    }
  else
    {
136
      arrayAddArray(len, dest, src);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
139
    }
}

140
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
Add2Vectors(double *dest, const double *restrict srcA, const double *restrict srcB, size_t len)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
  for (size_t i = 0; i < len; i++) dest[i] = srcA[i] + srcB[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
145
}

146
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
Sub2Vectors(double *dest, const double *restrict srcA, const double *restrict srcB, size_t len)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
  for (size_t i = 0; i < len; i++) dest[i] = srcA[i] - srcB[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150
151
}

152
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
MultVectorScalar(double *dest, const double *restrict src, double factor, size_t len, size_t nmiss, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
155
  if (nmiss)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
    {
157
      for (size_t i = 0; i < len; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
        dest[i] = (IS_EQUAL(src[i], missval)) ? missval : src[i] * factor;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160
161
    }
  else
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
      for (size_t i = 0; i < len; i++) dest[i] = src[i] * factor;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
164
165
    }
}

166
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
DivVectorIvector(double *dest, const double *restrict src, const int *samp, size_t len, size_t *nmiss, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
169
170
{
  *nmiss = 0;

171
  for (size_t i = 0; i < len; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
    {
173
174
175
176
177
      if (IS_EQUAL(src[i], missval) || samp[i] == 0)
        {
          dest[i] = missval;
          *nmiss = *nmiss + 1;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
      else
179
        dest[i] = src[i] / samp[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181
182
    }
}

183
184
void
after_read_vct(const char *vctfile, double **vct, int *nvct)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
187
188
189
190
{
  int n;
  char line[1024];
  double va, vb;

  FILE *fp = fopen(vctfile, "r");
191
  if (fp == nullptr) cdoSysError("Open failed on %s", vctfile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192

Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
  while (fgets(line, 1023, fp)) nvct++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195

  *nvct *= 2;
196
  *vct = (double *) Malloc(*nvct * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
198

  rewind(fp);
199
  for (int i = 0; i < *nvct / 2; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
201
202
    {
      fgets(line, 1023, fp);
      sscanf(line, "%d %lg %lg", &n, &va, &vb);
203
204
      *vct[i] = va;
      *vct[i + *nvct / 2] = vb;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
  fprintf(stdout, "  Reading VCT for %d hybrid levels from file %s\n", *nvct / 2 - 1, vctfile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
209
210

  fclose(fp);
}

211
void
212
after_gp2sp(const AfterControl &globs, struct Variable *vars, int ccode)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
{
  struct Variable *var = &vars[ccode];
215

Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
  if (var->spectral == nullptr)
217
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
      if (var->hybrid == nullptr)
219
220
221
222
223
        {
          fprintf(stderr, "%d.hybrid not found\n", ccode);
          exit(99);
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
      if (var->fourier == nullptr)
225
        {
226
          long fieldSize = globs.DimFC * var->hlev;
227
          var->fourier = alloc_dp(fieldSize, "gp2sp.fourier");
228
          after_GP2FC(var->hybrid, var->fourier, globs.Latitudes, globs.Longitudes, var->hlev, globs.Fouriers);
229
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230

231
232
      var->spectral = alloc_dp(globs.Dim3SP, "gp2sp.spectral");
      fc2sp(var->fourier, var->spectral, globs.pold, var->hlev, globs.Latitudes, globs.Fouriers, globs.Truncation);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
234
235
    }
}

236
237
void
after_GP2FC(double *gp, double *fc, long nlat, long nlon, long nlev, long nfc)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
{
  static long ifax[10];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
240
  static double *trig = nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241

242
  if (ifax[9] != nlon)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
    {
244
      if (trig) Free(trig);
245
      trig = (double *) Malloc(nlon * sizeof(double));
246
247
      int status = fft_set(trig, ifax, nlon);
      if (status < 0) exit(1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
252
    }

  gp2fc(trig, ifax, gp, fc, nlat, nlon, nlev, nfc);
}

253
254
void
after_FC2GP(double *fc, double *gp, long nlat, long nlon, long nlev, long nfc)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
256
{
  static long ifax[10];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
257
  static double *trig = nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258

259
  if (ifax[9] != nlon)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
    {
261
      if (trig) Free(trig);
262
      trig = (double *) Malloc(nlon * sizeof(double));
263
264
      int status = fft_set(trig, ifax, nlon);
      if (status < 0) exit(1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
266
267
268
269
    }

  fc2gp(trig, ifax, fc, gp, nlat, nlon, nlev, nfc);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
// HUMTEST
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271

272
static void
273
sh2rh(int AnalysisData, double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level, double *fullpresshybrid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274
{
275
276
277
278
279
280
281
  int lp, i;
  int lpi, lfp;
  double es, qsatr;

  /* ***************************************************** */
  /* Define constants for calculation in presence of water */
  /* ***************************************************** */
282
283
  constexpr double RGAMW = (C_RCW - C_RCPV) / C_RV;
  constexpr double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
  const double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT);
285
286
287
288
289
290
291

  /* ***************************************************** */
  /* Define constants for calculation in presence of  ice  */
  /* ***************************************************** */
  /*
  const double RGAMS = (C_RCS - C_RCPV) / C_RV;
  const double RBETS = C_RLSTT / C_RV + RGAMS * C_RTT;
292
  const double RALPS = std::log(C_RESTT) + RBETS / C_RTT + RGAMS * std::log(C_RTT);
293
  */
294
  const double *fullp = AnalysisData ? level : fullpresshybrid;
295
296
297
298
299
300

  /***************************************************/
  /* Diagnostics of saturation water vapour pressure */
  /* over ice makes no sense, therefore ...          */
  /* Hint of Michael Ponater                08.10.97 */
  /***************************************************/
301
302
  constexpr double RGAM = RGAMW;
  constexpr double RBET = RBETW;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
303
  const double RALP = RALPW;
304
305
306
307
308
309
310
311
312
313
314
315
316
  for (lp = 0; lp < lev; lp++)
    {
      for (i = 0; i < dimgpout; i++)
        {
          lpi = lp * dimgpout + i;
          lfp = (1 - AnalysisData) * lpi + AnalysisData * lp;
          /*
          if (t[lpi] < C_RTT) {
            RGAM = RGAMS; RBET = RBETS; RALP = RALPS;
          } else {
            RGAM = RGAMW; RBET = RBETW; RALP = RALPW;
          }
          */
317
          es = (std::exp(RALP - RBET / t[lpi] - RGAM * std::log(t[lpi]))) / fullp[lfp];
318
319
320
321
322
          // qsat = es / (1. + C_RETV * (1. - es));
          qsatr = (1. + C_RETV * (1. - es)) / es;
          rhum[lpi] = sphum[lpi] * 100. * qsatr;
        }
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
324
}

325
static void
326
rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
{
328
329
330
331
332
333
334
  int lp, i;
  int lpi;
  double es, qsat;

  /* ***************************************************** */
  /* Define constants for calculation in presence of water */
  /* ***************************************************** */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
336
  constexpr double RGAMW = (C_RCW - C_RCPV) / C_RV;
  constexpr double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
  const double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT);
338
339
340
341
342
343

  /* ***************************************************** */
  /* Define constants for calculation in presence of  ice  */
  /* ***************************************************** */
  // const double RGAMS = (C_RCS - C_RCPV) / C_RV;
  // const double RBETS = C_RLSTT / C_RV + RGAMS * C_RTT;
344
  // const double RALPS = std::log(C_RESTT) + RBETS / C_RTT + RGAMS * std::log(C_RTT);
345
346
347
348
349
350
351

  /***************************************************/
  /* Diagnostics of saturation water vapour pressure */
  /* over ice makes no sense, therefore ...          */
  /* Hint of Michael Ponater                08.10.97 */
  /***************************************************/

Uwe Schulzweida's avatar
Uwe Schulzweida committed
352
353
  constexpr double RGAM = RGAMW;
  constexpr double RBET = RBETW;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
354
  const double RALP = RALPW;
355
356
357
358
359
360
361
362
363
364
  for (lp = 0; lp < lev; lp++)
    {
      for (i = 0; i < dimgpout; i++)
        {
          lpi = lp * dimgpout + i;
          /*       if (t[lpi] < C_RTT) { */
          /*          RGAM = RGAMS; RBET = RBETS; RALP = RALPS; */
          /*       }  else { */
          /*          RGAM = RGAMW; RBET = RBETW; RALP = RALPW; */
          /*       } */
365
          es = (std::exp(RALP - RBET / t[lpi] - RGAM * std::log(t[lpi]))) / level[lp];
366
367
368
369
          qsat = es / (1. + C_RETV * (1. - es));
          sphum[lpi] = rhum[lpi] * qsat / 100.;
        }
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
371
}

372
void
373
after_FCrh2FCsh(const AfterControl &globs, struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
{
375
  const long fieldSize = globs.DimGP * globs.NumLevelRequest;
376

Uwe Schulzweida's avatar
Uwe Schulzweida committed
377
378
379
  if (vars[RHUMIDITY].grid == nullptr) vars[RHUMIDITY].grid = alloc_dp(fieldSize, "vars[RHUMIDITY].grid");
  if (vars[TEMPERATURE].grid == nullptr) vars[TEMPERATURE].grid = alloc_dp(fieldSize, "vars[TEMPERATURE].grid");
  if (vars[HUMIDITY].grid == nullptr) vars[HUMIDITY].grid = alloc_dp(fieldSize, "vars[HUMIDITY].grid");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380

381
382
383
384
  after_FC2GP(vars[RHUMIDITY].fourier, vars[RHUMIDITY].grid, globs.Latitudes, globs.Longitudes, vars[RHUMIDITY].plev,
              globs.Fouriers);
  after_FC2GP(vars[TEMPERATURE].fourier, vars[TEMPERATURE].grid, globs.Latitudes, globs.Longitudes, vars[TEMPERATURE].plev,
              globs.Fouriers);
385

386
387
  rh2sh(vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest, globs.DimGP,
        globs.LevelRequest);
388

389
390
  after_GP2FC(vars[HUMIDITY].grid, vars[HUMIDITY].fourier, globs.Latitudes, globs.Longitudes, vars[HUMIDITY].plev,
              globs.Fouriers);
391
392
393
394

  vars[HUMIDITY].grid = (double *) FreeMemory(vars[HUMIDITY].grid);
  vars[RHUMIDITY].grid = (double *) FreeMemory(vars[RHUMIDITY].grid);
  vars[TEMPERATURE].grid = (double *) FreeMemory(vars[TEMPERATURE].grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
}

397
void
398
after_SPuv2SPdv(const AfterControl &globs, struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
399
{
400
401
402
403
  double *Div, *DivOut, *Vor, *VorOut;

  Div = DivOut = vars[DIVERGENCE].spectral;
  Vor = VorOut = vars[VORTICITY].spectral;
404
  const long fieldSize = globs.DimFC * globs.NumLevelRequest;
405

Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
407
  if (vars[U_WIND].fourier == nullptr) vars[U_WIND].fourier = alloc_dp(fieldSize, "vars[U_WIND].fourier");
  if (vars[V_WIND].fourier == nullptr) vars[V_WIND].fourier = alloc_dp(fieldSize, "vars[V_WIND].fourier");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408

409
410
411
412
413
414
  sp2fc(vars[U_WIND].spectral, vars[U_WIND].fourier, globs.poli, globs.NumLevelRequest, globs.Latitudes, globs.Fouriers,
        globs.Truncation);
  sp2fc(vars[V_WIND].spectral, vars[V_WIND].fourier, globs.poli, globs.NumLevelRequest, globs.Latitudes, globs.Fouriers,
        globs.Truncation);
  uv2dv(vars[U_WIND].fourier, vars[V_WIND].fourier, Div, Vor, globs.pol2, globs.pol3, globs.NumLevelRequest, globs.Latitudes,
        globs.Truncation);
415
416
417
418

  vars[U_WIND].fourier = (double *) FreeMemory(vars[U_WIND].fourier);
  vars[V_WIND].fourier = (double *) FreeMemory(vars[V_WIND].fourier);

419
  for (int i = 0; i < globs.NumLevelRequest; ++i)
420
    {
421
422
423
424
425
426
      sp2sp(Div, globs.Truncation, DivOut, globs.Truncation);
      sp2sp(Vor, globs.Truncation, VorOut, globs.Truncation);
      Div += globs.DimSP;
      Vor += globs.DimSP;
      DivOut += globs.DimSP;
      VorOut += globs.DimSP;
427
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
429
}

430
void
431
after_FCsh2FCrh(const AfterControl &globs, struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
432
{
433
  const long fieldSize = globs.DimGP * globs.NumLevelRequest;
434

Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
436
437
  if (vars[RHUMIDITY].grid == nullptr) vars[RHUMIDITY].grid = alloc_dp(fieldSize, "vars[RHUMIDITY].grid");
  if (vars[TEMPERATURE].grid == nullptr) vars[TEMPERATURE].grid = alloc_dp(fieldSize, "vars[TEMPERATURE].grid");
  if (vars[HUMIDITY].grid == nullptr) vars[HUMIDITY].grid = alloc_dp(fieldSize, "vars[HUMIDITY].grid");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438

439
440
441
442
  after_FC2GP(vars[HUMIDITY].fourier, vars[HUMIDITY].grid, globs.Latitudes, globs.Longitudes, vars[HUMIDITY].plev,
              globs.Fouriers);
  after_FC2GP(vars[TEMPERATURE].fourier, vars[TEMPERATURE].grid, globs.Latitudes, globs.Longitudes, vars[TEMPERATURE].plev,
              globs.Fouriers);
443

444
445
  sh2rh(globs.AnalysisData, vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest,
        globs.DimGP, globs.LevelRequest, vars[FULL_PRESS].hybrid);
446

447
448
  after_GP2FC(vars[RHUMIDITY].grid, vars[RHUMIDITY].fourier, globs.Latitudes, globs.Longitudes, vars[RHUMIDITY].plev,
              globs.Fouriers);
449
450
451
452

  vars[HUMIDITY].grid = (double *) FreeMemory(vars[HUMIDITY].grid);
  vars[RHUMIDITY].grid = (double *) FreeMemory(vars[RHUMIDITY].grid);
  vars[TEMPERATURE].grid = (double *) FreeMemory(vars[TEMPERATURE].grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
454
455
}
/* ENDE HUMTEST */

456
457
static void
CheckAnalyses(struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
458
{
459
  for (int code = 0; code < 272; code++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
461
    if (vars[code].needed && code != DIVERGENCE && code != VORTICITY && code != STREAM && code != U_WIND && code != HUMIDITY
        && code != VELOPOT && code != V_WIND && code != RHUMIDITY && code != GEOPOTHEIGHT && code != PS
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
        && vars[code].spectral == nullptr && vars[code].grid == nullptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463
      {
464
        afterAbort("Code  %3d not found", code);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
465
466
467
468
      }
}

/*  Process Pressure Level data */
469
void
470
after_processPL(AfterControl &globs, struct Variable *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
472
473
474
{
  int code, l;
  long fieldSize;
  int lindex, nlevel;
475
  long offset;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476

477
478
  globs.MeanCount++;
  globs.TermCount++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479

480
  if (globs.MeanCount == 1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
    {
482
      if (globs.Debug) fprintf(stderr, "CheckAnalyses: %d %d\n", globs.TermCount, globs.MeanCount);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
      CheckAnalyses(vars);
484
      globs.StartDate = globs.OldDate;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
    }
486
  if (globs.TermCount > 120) globs.Debug = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
487
488
489
490
491

  /* ============================== */
  /* Computations in spectral space */
  /* ============================== */

492
  if (vars[TEMPERATURE].needed)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
494
    {
      vars[TEMPERATURE].hlev = 2;
495
      vars[TEMPERATURE].plev = globs.NumLevelRequest;
496
      vars[TEMPERATURE].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
498
    }

499
  if (vars[GEOPOTHEIGHT].comp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
501
    {
      vars[GEOPOTHEIGHT].hlev = 2;
502
      vars[GEOPOTHEIGHT].plev = globs.NumLevelRequest;
503
      vars[GEOPOTHEIGHT].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
505
    }

506
  if (vars[GEOPOTHEIGHT].comp && vars[GEOPOTENTIAL].detected)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
508
      if (vars[GEOPOTHEIGHT].spectral == nullptr)
509
510
        vars[GEOPOTHEIGHT].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "GEOPOTHEIGHT.spectral");
      MultVectorScalar(vars[GEOPOTHEIGHT].spectral, vars[GEOPOTENTIAL].spectral, C_RG, globs.DimSP * globs.NumLevelRequest, 0, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
513
      vars[GEOPOTENTIAL].needed = vars[GEOPOTENTIAL].selected;
    }

514
  if (globs.Type == 50 && vars[HUMIDITY].needed && vars[HUMIDITY].spectral == nullptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
515
    {
516
      vars[HUMIDITY].plev = globs.NumLevelRequest;
517
      vars[HUMIDITY].sfit = true;
518
      vars[HUMIDITY].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[HUMIDITY].spectral");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
      /*
520
        SPrh2SPsh();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
521
522
523
524
525
      */
      vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
      vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
526
  if (vars[U_WIND].spectral && vars[V_WIND].spectral && (vars[DIVERGENCE].comp || vars[VORTICITY].comp))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
    {
528
      vars[DIVERGENCE].hlev = vars[VORTICITY].hlev = 2;
529
      vars[DIVERGENCE].plev = vars[VORTICITY].plev = globs.NumLevelRequest;
530
      vars[DIVERGENCE].sfit = vars[VORTICITY].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
531
      if (vars[DIVERGENCE].spectral == nullptr)
532
        vars[DIVERGENCE].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[DIVERGENCE].spectral");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
      if (vars[VORTICITY].spectral == nullptr)
534
        vars[VORTICITY].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VORTICITY].spectral");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
535
536
537
      after_SPuv2SPdv(globs, vars);
    }

538
  if (vars[U_WIND].comp || vars[V_WIND].comp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
539
    {
540
      vars[U_WIND].hlev = vars[V_WIND].hlev = 2;
541
      vars[U_WIND].plev = vars[V_WIND].plev = globs.NumLevelRequest;
542
      vars[U_WIND].sfit = vars[V_WIND].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
543
      if (vars[U_WIND].spectral == nullptr)
544
        vars[U_WIND].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[U_WIND].spectral");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
      if (vars[V_WIND].spectral == nullptr)
546
547
548
        vars[V_WIND].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[V_WIND].spectral");
      dv2uv(vars[DIVERGENCE].spectral, vars[VORTICITY].spectral, vars[U_WIND].spectral, vars[V_WIND].spectral, globs.dv2uv_f1,
            globs.dv2uv_f2, globs.Truncation, globs.DimSP, globs.NumLevelRequest);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
    }

551
  if (vars[VELOPOT].comp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
    {
553
      vars[VELOPOT].hlev = 2;
554
      vars[VELOPOT].plev = globs.NumLevelRequest;
555
      vars[VELOPOT].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
      if (vars[VELOPOT].spectral == nullptr)
557
558
        vars[VELOPOT].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VELOPOT].spectral");
      dv2ps(vars[DIVERGENCE].spectral, vars[VELOPOT].spectral, globs.NumLevelRequest, globs.Truncation);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559
560
    }

561
  if (vars[STREAM].comp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
562
    {
563
      vars[STREAM].hlev = 2;
564
      vars[STREAM].plev = globs.NumLevelRequest;
565
      vars[STREAM].sfit = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
      if (vars[STREAM].spectral == nullptr)
567
568
        vars[STREAM].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[STREAM].spectral");
      dv2ps(vars[VORTICITY].spectral, vars[STREAM].spectral, globs.NumLevelRequest, globs.Truncation);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
569
570
571
572
573
574
    }

  /* --------------------------- */
  /*  Output of spectral fields  */
  /* --------------------------- */

575
  if (globs.Type == 50)
576
577
578
579
    {
      for (code = 0; code < MaxCodes; code++)
        if (vars[code].selected)
          {
580
            if (!vars[code].spectral) afterAbort("Code %d not available on spectral space!", code);
581
582
583
584

            nlevel = zaxisInqSize(vars[code].ozaxisID);
            for (lindex = 0; lindex < nlevel; lindex++)
              {
585
586
587
                offset = lindex * globs.DimSP;
                streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
                streamWriteRecord(globs.ostreamID, vars[code].spectral + offset, 0);
588
589
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
590
591
592
593
594
595
596
597
598
599

      FreeSpectral(vars);
      return;
    }

  /* =============================== */
  /* Transformation to fourier space */
  /* Computations in fourier space   */
  /* =============================== */

600
  if (globs.Type >= 60)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
601
602
    {
      for (code = 0; code < MaxCodes; code++)
603
604
        if (vars[code].needed && vars[code].spectral)
          {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
            if (vars[code].fourier == nullptr)
606
              {
607
                fieldSize = vars[code].plev * globs.DimFC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
608
                vars[code].fourier = alloc_dp(fieldSize, FieldName(code, "fourier"));
609
              }
610
611
            sp2fc(vars[code].spectral, vars[code].fourier, globs.poli, vars[code].plev, globs.Latitudes, globs.Fouriers,
                  globs.Truncation);
612
613
          }
      if (vars[U_WIND].needed && vars[U_WIND].fourier)
614
        scaluv(vars[U_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
615
      if (vars[V_WIND].needed && vars[V_WIND].fourier)
616
        scaluv(vars[V_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
617
618

      /* HUMTEST */
619
      if (globs.Type < 70 && vars[HUMIDITY].needed && vars[HUMIDITY].fourier == nullptr)
620
        {
621
          vars[HUMIDITY].plev = globs.NumLevelRequest;
622
          vars[HUMIDITY].sfit = true;
623
          vars[HUMIDITY].fourier = alloc_dp(globs.DimFC * globs.NumLevelRequest, "vars[HUMIDITY].fourier");
624
625
626
627
628
629
630

          after_FCrh2FCsh(globs, vars);

          vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
          vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
        }

631
      if (globs.Type < 70 && vars[RHUMIDITY].needed && vars[RHUMIDITY].fourier == nullptr)
632
        {
633
          vars[RHUMIDITY].plev = globs.NumLevelRequest;
634
          vars[RHUMIDITY].sfit = true;
635
          vars[RHUMIDITY].fourier = alloc_dp(globs.DimFC * globs.NumLevelRequest, "vars[RHUMIDITY].fourier");
636
637
638
639
640
641

          after_FCsh2FCrh(globs, vars);

          vars[HUMIDITY].needed = vars[HUMIDITY].selected;
          vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
643
644
645
646
647
648
649
650
      /* ENDE HUMTEST */
    }

  FreeSpectral(vars);

  /* -------------------------- */
  /*  Output of fourier fields  */
  /* -------------------------- */

651
  if (globs.Type == 60)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
652
    {
653
654
655
      for (code = 0; code < MaxCodes; code++)
        if (vars[code].selected)
          {
656
            if (!vars[code].fourier) afterAbort("Code %d not available on fourier space!", code);
657
658
659
660

            nlevel = zaxisInqSize(vars[code].ozaxisID);
            for (lindex = 0; lindex < nlevel; lindex++)
              {
661
662
663
                offset = lindex * globs.DimFC;
                streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
                streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
664
665
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
666
667
668
669
670
671
672
673
674

      FreeFourier(vars);
      return;
    }

  /* ----------------------- */
  /*  Output of zonal means  */
  /* ----------------------- */

675
  if (globs.Type == 61)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
    {
677
678
679
      for (code = 0; code < MaxCodes; code++)
        if (vars[code].selected)
          {
680
            if (!vars[code].fourier) afterAbort("Code %d not available on zonal mean!", code);
681
682
683
684

            nlevel = zaxisInqSize(vars[code].ozaxisID);
            for (lindex = 0; lindex < nlevel; lindex++)
              {
685
686
687
                offset = lindex * globs.DimFC;
                streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
                streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
688
689
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
690
691
692

      FreeFourier(vars);
      return;
693
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
695
696
697
698

  /* ============================ */
  /* Transformation to gridpoints */
  /* ============================ */

699
  if (vars[PS].comp && vars[LNPS].grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
700
    {
701
702
      if (vars[PS].grid == nullptr) vars[PS].grid = alloc_dp(globs.DimGP, "Ps");
      for (l = 0; l < globs.DimGP; l++) vars[PS].grid[l] = std::exp(vars[LNPS].grid[l]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703
704
    }

705
  if (globs.Type >= 70)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706
707
    {
      for (code = 0; code < MaxCodes; code++)
708
709
        if (vars[code].needed && vars[code].fourier)
          {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
710
            if (vars[code].grid == nullptr)
711
              {
712
                fieldSize = vars[code].plev * globs.DimGP;
713
714
                vars[code].grid = alloc_dp(fieldSize, FieldName(code, "grid"));
              }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
715

716
            after_FC2GP(vars[code].fourier, vars[code].grid, globs.Latitudes, globs.Longitudes, vars[code].plev, globs.Fouriers);
717
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
718
719
    }

720
  FreeFourier(vars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
722
723
724
725
726

  /* HUMTEST */
  /* -------------------------------- */
  /*  Computation in gridpoint space  */
  /* -------------------------------- */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
  if (vars[RHUMIDITY].needed && vars[RHUMIDITY].grid == nullptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
728
    {
729
      vars[RHUMIDITY].plev = globs.NumLevelRequest;
730
      vars[RHUMIDITY].sfit = true;
731
732
733
      vars[RHUMIDITY].grid = alloc_dp(globs.DimGP * globs.NumLevelRequest, "vars[RHUMIDITY].grid");
      sh2rh(globs.AnalysisData, vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest,
            globs.DimGP, globs.LevelRequest, vars[FULL_PRESS].hybrid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
734
735
736
737
      vars[HUMIDITY].needed = vars[HUMIDITY].selected;
      vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
738
  if (vars[HUMIDITY].needed && vars[HUMIDITY].grid == nullptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
    {
740
      vars[HUMIDITY].plev = globs.NumLevelRequest;
741
      vars[HUMIDITY].sfit = true;
742
743
744
      vars[HUMIDITY].grid = alloc_dp(globs.DimGP * globs.NumLevelRequest, "vars[HUMIDITY].grid");
      rh2sh(vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest, globs.DimGP,
            globs.LevelRequest);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
746
747
748
749
750
751
752
753
      vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
      vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
    }
  /* HUMTEST ENDE */

  /* -------------------------- */
  /*  Computation of Means      */
  /* -------------------------- */

754
  if (globs.Mean)
755
756
757
    for (code = 0; code < MaxCodes; code++)
      if (vars[code].needed && vars[code].grid)
        {
758
          fieldSize = globs.DimGP * vars[code].plev;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759
          if (vars[code].mean == nullptr) vars[code].mean = alloc_dp(fieldSize, FieldName(code, "mean"));
760

761
          if (globs.MeanCount == 1)
762
763
            arrayCopy(fieldSize, vars[code].grid, vars[code].mean);
          else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
            AddVector(vars[code].mean, vars[code].grid, fieldSize, &vars[code].nmiss, vars[code].missval);
765

766
          if (globs.EndOfInterval)
767
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
768
              if (vars[code].samp == nullptr)
769
                MultVectorScalar(vars[code].mean, vars[code].mean, 1.0 / globs.MeanCount, fieldSize, vars[code].nmiss,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
770
                                 vars[code].missval);
771
              else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
772
                DivVectorIvector(vars[code].mean, vars[code].mean, vars[code].samp, fieldSize, &vars[code].nmiss,
773
774
775
                                 vars[code].missval);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
776
777
778
779
780

  /* -------------------------- */
  /*  Computation of Variances  */
  /* -------------------------- */

781
  if (globs.Mean > 1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
782
783
    for (code = 0; code < MaxCodes; code++)
      if (vars[code].needed && vars[code].mean)
784
        {
785
          fieldSize = globs.DimGP * vars[code].plev;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
786
          if (vars[code].variance == nullptr) vars[code].variance = alloc_dp(fieldSize, FieldName(code, "var"));
787
          if (globs.MeanCount == 1)
788
789
790
791
            IniQuaSum(vars[code].variance, vars[code].grid, fieldSize);
          else
            AddQuaSum(vars[code].variance, vars[code].grid, fieldSize);

792
          if (globs.EndOfInterval) VarQuaSum(vars[code].variance, vars[code].mean, fieldSize, globs.MeanCount);
793
794
        }

795
  if (globs.Mean && !globs.EndOfInterval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
796
797
798
799
800
801
802
803
804
    {
      FreeGrid(vars);
      return;
    }

  /* ---------------------------------------------- */
  /*  Output of pressure level means and variances  */
  /* ---------------------------------------------- */

805
  if (globs.Type == 70 && globs.Mean && globs.EndOfInterval)
806
807
808
809
810
811
812
    {
      for (code = 0; code < MaxCodes; code++)
        if (vars[code].selected)
          {
            nlevel = zaxisInqSize(vars[code].ozaxisID);
            for (lindex = 0; lindex < nlevel; lindex++)
              {
813
814
                offset = lindex * globs.DimGP;
                if (globs.Mean != 2)
815
                  {
816
817
                    streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
                    streamWriteRecord(globs.ostreamID, vars[code].mean + offset, vars[code].nmiss);
818
                  }
819
                if (globs.Mean >= 2)
820
                  {
821
822
                    streamDefRecord(globs.ostreamID2, vars[code].ovarID2, lindex);
                    streamWriteRecord(globs.ostreamID2, vars[code].variance + offset, vars[code].nmiss);
823
824
825
                  }
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
826
827
828
829
830
831
832
833
834
835

      FreeSamp(vars);
      FreeGrid(vars);
      return;
    }

  /* -------------------------------- */
  /*  Output of pressure level grids  */
  /* -------------------------------- */

836
  if (globs.Type == 70)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
837
    {
838
839
840
841
842
843
      for (code = 0; code < MaxCodes; code++)
        if (vars[code].selected)
          {
            nlevel = zaxisInqSize(vars[code].ozaxisID);
            for (lindex = 0; lindex < nlevel; lindex++)
              {
844
845
846
                offset = lindex * globs.DimGP;
                streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
                streamWriteRecord(globs.ostreamID, vars[code].grid + offset, vars[code].nmiss);
847
848
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
849
850
851
852
853
854

      FreeGrid(vars);
      return;
    }
}

855
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
856
theta(double *pthetaf, double *pthetah, double *ph, double *ps, double *tf, double *ts, int levels, int dimgp, int dim3gp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
857
{
858
859
860
861
862
  double *thetah = pthetah;
  double *thetaf = pthetaf;

  const double kappa = PlanetRD / C_RCPD;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
863
  for (int h = 0; h < dimgp; h++) thetah[h] = 0.0;
864
865
866
  thetah += dimgp;
  for (int l = 0; l < levels - 1; l++)
    {
867
      for (int h = 0; h < dimgp; h++) thetah[h] = 0.5 * (tf[h] + tf[h + dimgp]) * std::pow((ps[h] / ph[h]), kappa);
868
869
870
871
872
873
874
875

      ph += dimgp;
      tf += dimgp;
      thetah += dimgp;
    }

  arrayCopy(dimgp, ts, thetah);
  thetah = pthetah;