Afterburner.cc 57.2 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
20
21
22
23
24
25
26
27
28
/* ============================================================= */
/*                                                               */
/* postprocessing program for ECHAM data and ECMWF analysis data */
/*                                                               */
/* Luis     Kornblueh   - MPI    Hamburg                         */
/* Uwe      Schulzweida - MPI    Hamburg                         */
/* Arno     Hellbach    - DKRZ   Hamburg                         */
/* Edilbert Kirk        - MI Uni Hamburg                         */
/* Michael  Ponater     - DLR    Oberpfaffenhofen                */
/*                                                               */
/* ============================================================= */

29
#ifdef HAVE_CONFIG_H
30
31
32
#include "config.h"
#endif

33
34
35
36
#ifdef _OPENMP
#include <omp.h>
#endif

37
38
#include <cdi.h>

39
#include "cdo_default_values.h"
40

41
#define streamOpenWrite cdoOpenWrite
42
43
#define streamDefVlist cdoDefVlist
#define streamDefTimestep cdoDefTimestep
44
45
46
47

#include "afterburner.h"
#include "constants.h"
#include "compare.h"
48
#include "cdo_options.h"
49
#include "mpmo_color.h"
50
#include "cdi_lockedIO.h"
51

Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
int scan_par_obsolete(char *namelist, const char *name, int def);
53
int scan_par(int verbose, char *namelist, const char *name, int def);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
void scan_code(char *namelist, struct Variable *vars, int maxCodes, int *numCodes);
55
int scan_time(int verbose, char *namelist, int *hours, int max_hours);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
void scan_darray(char *namelist, const char *name, double *values, int maxValues, int *numValues);
57

58
char zaxistypename[CDI_MAX_NAME];
59

Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
struct RARG
61
{
62
63
  int lana, nrecs;
  struct Variable *vars;
64
  AfterControl *globs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
};
66

67
cdo::Task *afterReadTask = nullptr;
68

69
static bool lstdout = true;
70

71
72
73
74
75
76
77
static int Source = 0;

static int ofiletype = -1;

static int DataType = -1;

static char *filename;
78
static const char **ifiles;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
static char *ifile = nullptr;
static const char *ofile = nullptr;
81

82
83
static int ofileidx = 0;

84
static int specGridID = -1;
85
86
87
88
static int gaussGridID = -1;
static int iVertID = -1;
static int oVertID = -1;

89
static bool Lhybrid2pressure = false;
90

91
static int TsID;
92
bool afterReadAsync = true;
93

94
95
96
97
#define TIMESTEP_INTERVAL -1
#define MONTHLY_INTERVAL 0
#define DAILY_INTERVAL 1
#define UNLIM_INTERVAL 2
98

99
#define MaxHours 24
100
static int nrqh;
101
static int hours[MaxHours + 1];
102

103
static std::vector<double> LevelFound;
104

Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
106
107
108
109
110
111
static inline bool
filetype_is_netcdf(int filetype)
{
  return filetype == CDI_FILETYPE_NC || filetype == CDI_FILETYPE_NC2 || filetype == CDI_FILETYPE_NC4
    || filetype == CDI_FILETYPE_NC4C || filetype == CDI_FILETYPE_NC5;
}

112
113
static void
lprintf(FILE *fp)
114
115
116
117
118
{
  int num = 67;
  int cval = '-';

  fprintf(fp, " ");
119
  for (int inum = 0; inum < num; inum++) fprintf(fp, "%c", cval);
120
121
122
  fprintf(fp, "\n");
}

123
124
static void
FreeMean(struct Variable *vars)
125
{
126
127
  for (int code = 0; code < MaxCodes; code++)
    if (vars[code].mean)
128
      {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
        Free(vars[code].mean);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
        vars[code].mean = nullptr;
131
132
133
      }
}

134
static void
135
after_PostProcess(AfterControl *globs)
136
{
137
138
139
140
141
  if (globs->EndOfInterval)
    {
      if (lstdout)
        {
          if (globs->OutputInterval == DAILY_INTERVAL)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
            fprintf(stdout, " Processed Day %2d  Month %2d  Year %04d", globs->OldDate.dy, globs->OldDate.mo, globs->OldDate.yr);
143
          else if (globs->OutputInterval == MONTHLY_INTERVAL)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
            fprintf(stdout, " Processed Month %2d  Year %04d", globs->OldDate.mo, globs->OldDate.yr);
145
          else if (globs->OutputInterval == UNLIM_INTERVAL)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
147
            fprintf(stdout, " Processed range from %6.4d-%2.2d-%2.2d to %6.4d-%2.2d-%2.2d", globs->StartDate.yr,
                    globs->StartDate.mo, globs->StartDate.dy, globs->OldDate.yr, globs->OldDate.mo, globs->OldDate.dy);
148
149
150
151
152
153

          if (globs->Mean)
            fprintf(stdout, "  (Mean of %3d Terms)\n", globs->MeanCount);
          else
            fprintf(stdout, "   Terms %3d\n", globs->MeanCount);
        }
154

155
      globs->EndOfInterval = false;
156
157
158
159
160
161
162
      globs->MeanCount = 0;
    }
}

/* ================= */
/* switch input file */
/* ================= */
163
static void
164
after_SwitchFile(AfterControl *globs)
165
{
166
  bool echam4 = false;
167
  int n;
168
  char y3, y2, y1, y0;
169
170
  char m1, m0;
  char d1, d0;
171
172
173

  streamClose(globs->istreamID);

174
  if (globs->Multi > 0)
175
    {
176
      int i = strlen(ifile);
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
      if (i < 10)
        {
          fprintf(stderr, " Not a valid filename: %s \n", ifile);
          exit(1);
        }

      if (ifile[i - 3] == '.')
        {
          echam4 = true;
          y3 = ifile[i - 9];
          y2 = ifile[i - 8];
          y1 = ifile[i - 7];
          y0 = ifile[i - 6];
          m1 = ifile[i - 5];
          m0 = ifile[i - 4];
          d1 = ifile[i - 2];
          d0 = ifile[i - 1];
        }
195
      else
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
        {
          y3 = ifile[i - 6];
          y2 = ifile[i - 5];
          y1 = ifile[i - 4];
          y0 = ifile[i - 3];
          m1 = ifile[i - 2];
          m0 = ifile[i - 1];
          d1 = '0';
          d0 = '1';
        }

      for (n = 0; n < globs->DayIn; n++)
        {
          if (d0 == '9')
            {
              d0 = '0';
              d1++;
            }
          else
            d0++;
          if (d1 == '3' && d0 > '0')
            {
              d1 = '0';
              d0 = '1';
              if (m1 == '0')
                {
                  if (m0 == '9')
                    {
                      m0 = '0';
                      m1 = '1';
                    }
                  else
                    m0++;
                }
              else
                {
                  if (m0 < '2')
                    m0++;
                  else
                    {
                      m1 = '0';
                      m0 = '1';
                      y0++;
                      if (y0 > '9')
                        {
                          y0 = '0';
                          y1++;
                        }
                      if (y1 > '9')
                        {
                          y1 = (char) '0';
                          if (isdigit((int) y2))
                            y2++;
                          else
                            y2 = '1';
                          if (y2 > '9')
                            {
                              y2 = (char) '0';
                              if (isdigit((int) y3))
                                y3++;
                              else
                                y3 = '1';
                            }
                        }
                    }
                }
            }
        }

      if (echam4)
        {
          ifile[i - 9] = y3;
          ifile[i - 8] = y2;
          ifile[i - 7] = y1;
          ifile[i - 6] = y0;
          ifile[i - 5] = m1;
          ifile[i - 4] = m0;
          ifile[i - 2] = d1;
          ifile[i - 1] = d0;
        }
276
      else
277
278
279
280
281
282
283
284
        {
          ifile[i - 6] = y3;
          ifile[i - 5] = y2;
          ifile[i - 4] = y1;
          ifile[i - 3] = y0;
          ifile[i - 2] = m1;
          ifile[i - 1] = m0;
        }
285
286
287
288

      globs->Multi--;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
  if (globs->Nfiles > 0) ifile = (char *) ifiles[--globs->Nfiles];
290
291
292

  fprintf(stderr, " Continuation file: %s\n", ifile);

293
  globs->istreamID = streamOpenReadLocked(ifile);
294
295

  globs->ivlistID = streamInqVlist(globs->istreamID);
296
  globs->taxisID = vlistInqTaxis(globs->ivlistID);
297
298
}

299
static int64_t
300
after_getDate(struct Date datetime)
301
302
303
304
{
  return cdiEncodeDate(datetime.yr, datetime.mo, datetime.dy);
}

305
306
static int
after_getTime(struct Date datetime)
307
308
309
310
{
  return cdiEncodeTime(datetime.hr, datetime.mn, 0);
}

311
static void
312
after_setDateTime(struct Date *datetime, int64_t date, int time)
313
314
315
316
317
318
{
  int sec;
  cdiDecodeDate(date, &datetime->yr, &datetime->mo, &datetime->dy);
  cdiDecodeTime(time, &datetime->hr, &datetime->mn, &sec);
}

319
320
static void
after_printProcessStatus(int tsID)
321
{
322
  static bool counthead = false;
323

324
  if (tsID == -1)
325
    {
326
327
      if (stdout_is_tty)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
          fprintf(stdout, "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
329
330
          fflush(stdout);
        }
331

332
      counthead = false;
333
334
335
    }
  else
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
      if (!counthead)
337
338
        {
          if (stdout_is_tty) fprintf(stdout, " Process timestep :       ");
339

340
341
          counthead = true;
        }
342

343
344
345
346
347
      if (stdout_is_tty)
        {
          fprintf(stdout, "\b\b\b\b\b\b%6d", tsID);
          fflush(stdout);
        }
348
349
350
    }
}

351
static int
352
after_setNextDate(AfterControl *globs)
353
{
354
  int nrecs = 0;
355

356
  bool righttime = false;
357
  while (true)
358
359
    {
      nrecs = streamInqTimestep(globs->istreamID, TsID);
360
361
362
363
364
365
366
367
368
369
370
371
372
      if (nrecs == 0 && (globs->Multi > 0 || globs->Nfiles > 0))
        {
          if (lstdout) after_printProcessStatus(-1);

          after_SwitchFile(globs);

          if (globs->istreamID >= 0)
            {
              TsID = 0;
              nrecs = streamInqTimestep(globs->istreamID, TsID);
            }
        }
      if (nrecs == 0) break;
373

374
      const auto vdate = taxisInqVdate(globs->taxisID);
375
      const auto vtime = taxisInqVtime(globs->taxisID);
376
377
      after_setDateTime(&globs->NextDate, vdate, vtime);

378
379
380
381
382
383
      for (int i = 0; i < nrqh; i++)
        if (hours[i] < 0 || hours[i] == globs->NextDate.hr)
          {
            righttime = true;
            break;
          }
384

385
386
      if (righttime)
        break;
387
      else
388
        TsID += 1;
389
390
    }

391
  return nrecs;
392
393
394
395
}

static int num_recs = 0;

396
397
static void *
after_readTimestep(void *arg)
398
{
399
  int varID, gridID, zaxisID, levelID, timeID;
400
  size_t nmiss;
401
  auto rarg = (RARG *) arg;
402

403
404
405
406
  auto nrecs = rarg->nrecs;
  auto analysisData = rarg->lana;
  auto vars = rarg->vars;
  auto globs = rarg->globs;
407

408
  for (int code = 0; code < MaxCodes; code++) vars[code].nmiss0 = 0;
409

410
  for (int recID = 0; recID < nrecs; recID++)
411
412
413
    {
      streamInqRecord(globs->istreamID, &varID, &levelID);

414
      const auto code = vlistInqVarCode(globs->ivlistID, varID);
415
      if (code <= 0 || code >= MaxCodes) continue;
416
417
418

      /* Skip records containing unneeded codes */

419
      if (!vars[code].needed0) continue;
420
421
422

      vlistInqVar(globs->ivlistID, varID, &gridID, &zaxisID, &timeID);

423
      const auto leveltype = zaxisInqType(zaxisID);
424

425
426
      /* Skip records with unselected levels */

427
      int levelOffset = -1;
428
      /*
429
        if ( vars[code].ozaxisID != vars[code].izaxisID && ! Lhybrid2pressure )
430
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
      if ((vars[code].ozaxisID != vars[code].izaxisID) && (leveltype == ZAXIS_PRESSURE))
432
        {
433
          const auto level = (int) zaxisInqLevel(zaxisID, levelID);
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
          for (int i = 0; i < globs->NumLevelRequest; ++i)
            {
              if (IS_EQUAL(globs->LevelRequest[i], level))
                {
                  levelOffset = i;
                  break;
                }
            }

          if (levelOffset < 0) continue;

          zaxisID = vars[code].ozaxisID;
          levelID = levelOffset;
        }

      if (globs->Debug)
        {
          fprintf(stderr, "T%d", globs->Truncation);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
          fprintf(stderr, "  Code %3d   Level%6d   %6.4d-%2.2d-%2.2d  %2.2d:%2.2d:00\n", code,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
455
                  (int) zaxisInqLevel(zaxisID, levelID), globs->OldDate.yr, globs->OldDate.mo, globs->OldDate.dy, globs->OldDate.hr,
                  globs->OldDate.mn);
456
457
458
        }

      if (analysisData)
459
460
        {
          streamReadRecord(globs->istreamID, globs->Field, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
461
          after_AnalysisAddRecord(globs, vars, code, gridID, zaxisID, levelID, nmiss);
462
        }
463
      else
464
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
465
          double *dataptr = after_get_dataptr(vars, code, gridID, zaxisID, levelID);
466
          streamReadRecord(globs->istreamID, dataptr, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
          after_EchamAddRecord(globs, vars, code, gridID, zaxisID, levelID, nmiss);
468
        }
469

Uwe Schulzweida's avatar
Uwe Schulzweida committed
470
      if (iVertID != -1 && oVertID != -1 && (vars[code].izaxisID == iVertID)) vars[code].ozaxisID = oVertID;
471
472
473
    }

  TsID++;
474
475
476

  // printf("%3d  date = %d  time = %04d\n", TsID, vdate, vtime);

477
478
  num_recs = after_setNextDate(globs);

479
  return (void *) &num_recs;
480
481
}

482
static void
483
after_defineNextTimestep(AfterControl *globs)
484
485
{
  static int otsID = 0;
486
487
  const auto vdate = after_getDate(globs->OldDate);
  const auto vtime = after_getTime(globs->OldDate);
488
489
490
  taxisDefVdate(globs->taxisID2, vdate);
  taxisDefVtime(globs->taxisID2, vtime);

491
  if (globs->Mean != 2)
492
    {
493
494
      if (otsID == 0)
        {
495
          const auto nvars = vlistNvars(globs->ovlistID);
496
          if (nvars == 0) afterAbort("No variable selected!");
497
498
499
500
501
          vlistDefTaxis(globs->ovlistID, globs->taxisID2);
          streamDefVlist(globs->ostreamID, globs->ovlistID);
        }
      taxisDefNumavg(globs->taxisID2, globs->MeanCount + 1);
      streamDefTimestep(globs->ostreamID, otsID);
502
503
504
505
506
    }

  otsID++;
}

507
static void
508
after_setEndOfInterval(AfterControl *globs, int nrecs)
509
{
510
  if (nrecs == 0)
511
    {
512
      globs->EndOfInterval = true;
513
514
515
    }
  else
    {
516
517
518
519
520
      if (globs->OutputInterval == DAILY_INTERVAL)
        globs->EndOfInterval = globs->NewDate.dy != globs->OldDate.dy;
      else if (globs->OutputInterval == MONTHLY_INTERVAL)
        globs->EndOfInterval = globs->NewDate.mo != globs->OldDate.mo;
      else if (globs->OutputInterval == UNLIM_INTERVAL)
521
        globs->EndOfInterval = false;
522
      else
523
        afterAbort("output interval %d not implemented!", globs->OutputInterval);
524
525
526
    }
}

527
528
static void
after_moveTimestep(struct Variable *vars)
529
530
531
{
  int code;

532
  for (code = 0; code < MaxCodes; code++) vars[code].nmiss = vars[code].nmiss0;
533

534
535
  for (code = 0; code < MaxCodes; code++)
    if (vars[code].hybrid0)
536
      {
537
        vars[code].hybrid = vars[code].hybrid0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538
        vars[code].hybrid0 = nullptr;
539
540
      }

541
542
  for (code = 0; code < MaxCodes; code++)
    if (vars[code].spectral0)
543
      {
544
        vars[code].spectral = vars[code].spectral0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
        vars[code].spectral0 = nullptr;
546
547
      }

548
549
  for (code = 0; code < MaxCodes; code++)
    if (vars[code].grid0)
550
      {
551
        vars[code].grid = vars[code].grid0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
        vars[code].grid0 = nullptr;
553
554
555
      }
}

556
557
static void
after_check_content(struct Variable *vars, int timestep)
558
{
559
  for (int code = 0; code < 272; code++)
560
561
    {
      /*  if ( code == GEOPOTENTIAL ) continue; */
562
563
564
565
566
567
568
569
570
571
572
573
574
575
      if (code == SLP) continue;
      if (code == GEOPOTHEIGHT) continue;
      if (code == STREAM) continue;
      if (code == VELOPOT) continue;
      if (code == U_WIND) continue;
      if (code == V_WIND) continue;
      if (code == OMEGA) continue;
      if (code == RHUMIDITY) continue;
      if (code == LOW_CLOUD) continue;
      if (code == MID_CLOUD) continue;
      if (code == HIH_CLOUD) continue;
      if (code == PS) continue;
      if (code == HUMIDITY)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
          if (vars[code].needed && !vars[code].selected && vars[code].spectral == nullptr && vars[code].hybrid == nullptr)
577
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
578
              static bool lwarn = true;
579
              if (lwarn) cdoWarning("No humidity in data file, set to zero !");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
              lwarn = false;
581
              vars[code].needed = false;
582
583
            }
        }
584
      else
585
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
          if (vars[code].needed && !vars[code].comp && vars[code].spectral == nullptr && vars[code].hybrid == nullptr)
587
            {
588
              afterAbort("Code  %3d not found at timestep %d!", code, timestep);
589
590
            }
        }
591
592
593
    }
}

594
static void
595
after_control(AfterControl *globs, struct Variable *vars)
596
597
{
  int i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
  int nrecs = 0;
599
600
601
  int code;
  RARG rarg;

602
  if (afterReadAsync) afterReadTask = new cdo::Task;
603

604
  for (code = 0; code < MaxCodes; code++) vars[code].needed0 = vars[code].needed;
605
606
607

  TsID = 0;

608
  bool righttime = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
609
  while (true)
610
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
612
613
614
615
      nrecs = streamInqTimestep(globs->istreamID, TsID);
      if (nrecs <= 0) break;

      const auto vdate = taxisInqVdate(globs->taxisID);
      const auto vtime = taxisInqVtime(globs->taxisID);
616
617
618
      after_setDateTime(&globs->StartDate, vdate, vtime);
      after_setDateTime(&globs->NewDate, vdate, vtime);

619
620
621
622
623
624
      for (i = 0; i < nrqh; i++)
        if (hours[i] < 0 || hours[i] == globs->NewDate.hr)
          {
            righttime = true;
            break;
          }
625

Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
      if (righttime) break;
627

Uwe Schulzweida's avatar
Uwe Schulzweida committed
628
      TsID++;
629
630
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
631
632
633
634
635
  const auto taxis_is_relative = (taxisInqType(globs->taxisID) == TAXIS_RELATIVE);
  const auto rdate = taxis_is_relative ? taxisInqRdate(globs->taxisID) : after_getDate(globs->StartDate);
  const auto rtime = taxis_is_relative ? taxisInqRtime(globs->taxisID) : after_getTime(globs->StartDate);

  if (filetype_is_netcdf(ofiletype))
636
637
638
639
640
641
642
643
644
645
    {
      taxisDefCalendar(globs->taxisID2, CALENDAR_PROLEPTIC);
      taxisDefType(globs->taxisID2, TAXIS_RELATIVE);
      taxisDefTunit(globs->taxisID2, TUNIT_DAY);
      taxisDefRdate(globs->taxisID2, rdate);
      taxisDefRtime(globs->taxisID2, rtime);
    }

  globs->OldDate = globs->NewDate;

646
  bool tsFirst = true;
647

648
  while (nrecs > 0)
649
650
    {
      rarg.nrecs = nrecs;
651
652
      rarg.lana = globs->AnalysisData;
      rarg.vars = vars;
653
654
      rarg.globs = globs;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
655
656
      int status = -1;

657
      if (tsFirst || !afterReadAsync)
658
        {
659
          if (afterReadAsync)
660
            {
661
              afterReadTask->start(after_readTimestep, &rarg);
662
663
664
            }
          else
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
              status = *(int *) after_readTimestep(&rarg);
666
            }
667

668
          if (tsFirst && globs->Type > 0) after_legini_setup(globs, vars);
669

670
          if (afterReadAsync)
671
            {
672
              status = *(int *) afterReadTask->wait();
673
              if (status < 0) cdoAbort("after_readTimestep error! (status = %d)", status);
674
675
676
            }
          tsFirst = false;
        }
677
      else
678
        {
679
          status = *(int *) afterReadTask->wait();
680
          if (status < 0) cdoAbort("after_readTimestep error! (status = %d)", status);
681
        }
682

Uwe Schulzweida's avatar
Uwe Schulzweida committed
683
      nrecs = status;
684
685
686
687
688
689

      globs->MeanCount0 = globs->MeanCount;
      globs->NewDate = globs->NextDate;

      after_moveTimestep(vars);

690
      if (nrecs && afterReadAsync)
691
        {
692
          afterReadTask->start(after_readTimestep, &rarg);
693
        }
694
695
696

      after_setEndOfInterval(globs, nrecs);

697
698
699
700
701
      if (lstdout) after_printProcessStatus(TsID);

      if (lstdout && globs->EndOfInterval) after_printProcessStatus(-1);

      if (globs->Mean == 0 || globs->EndOfInterval)
702
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703
          if (!globs->AnalysisData) after_check_content(vars, globs->TermCount + 1);
704
705
          after_defineNextTimestep(globs);
        }
706
707
708

      if (globs->AnalysisData)
        after_processPL(globs, vars);
709
      else
710
        after_processML(globs, vars);
711
712

      after_PostProcess(globs);
713
714
715
716
717
718
719
720
721

      if (nrecs)
        {
          if (globs->AnalysisData)
            after_AnalysisDependencies(vars, MaxCodes);
          else
            after_EchamDependencies(vars, MaxCodes, globs->Type, Source);
        }

722
723
724
      globs->OldDate = globs->NewDate;
    }

725
  if (afterReadTask) delete afterReadTask;
726
727
}

728
static void
729
after_setLevel(AfterControl *globs)
730
731
{
  int k, l, found;
732
  bool removeLevel[MaxLevel];
733
  double level;
734
  bool checkLevel = true;
735
  // default pressure level
Uwe Schulzweida's avatar
Uwe Schulzweida committed
736
737
  long plevelDefault[]
      = { 100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 25000, 20000, 15000, 10000, 7000, 5000, 3000, 2000, 1000 };
738
  // default height level
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
  long hlevelDefault[] = { 0, 1000, 2000, 5000, 10000, 15000, 20000, 25000, 30000 };
740

741
742
  const int numplevelDefault = sizeof(plevelDefault) / sizeof(plevelDefault[0]);
  const int numhlevelDefault = sizeof(hlevelDefault) / sizeof(hlevelDefault[0]);
743
744

  if (iVertID != -1)
745
    if (zaxisInqType(iVertID) == ZAXIS_HYBRID && globs->Type > 20) Lhybrid2pressure = true;
746

747
  if (globs->Verbose) lprintf(stdout);
748

749
  if (globs->NumLevelRequest == 0)
750
    {
751
752
753
754
      if (iVertID == -1)
        {
          if (globs->Verbose) fprintf(stdout, " No level detected\n");
        }
755
      else
756
757
758
759
760
        {
          if (Lhybrid2pressure)
            {
              if (globs->unitsel == 0)
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
761
                  if (globs->Verbose) fprintf(stdout, " Default pressure level selected:\n");
762
                  globs->NumLevelRequest = numplevelDefault;
763
                  for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = plevelDefault[l];
764
765
766
767
768
                  oVertID = zaxisCreate(ZAXIS_PRESSURE, globs->NumLevelRequest);
                  zaxisDefLevels(oVertID, globs->LevelRequest);
                }
              else
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
769
                  if (globs->Verbose) fprintf(stdout, " Default height level selected:\n");
770
                  globs->NumLevelRequest = numhlevelDefault;
771
                  for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = hlevelDefault[l];
772
773
774
775
776
777
778
779
780
781
782
783
784
785
                  oVertID = zaxisCreate(ZAXIS_HEIGHT, globs->NumLevelRequest);
                  zaxisDefLevels(oVertID, globs->LevelRequest);
                }
            }
          else
            {
              if (globs->Verbose)
                {
                  if (zaxisInqType(iVertID) == ZAXIS_HYBRID)
                    fprintf(stdout, " All detected hybrid level selected:\n");
                  else
                    fprintf(stdout, " All detected pressure level selected:\n");
                }
              globs->NumLevelRequest = globs->NumLevelFound;
786
              for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = LevelFound[l];
787
788
789
              oVertID = iVertID;
            }
        }
790
      checkLevel = false;
791
792
793
    }
  else
    {
794
795
796
797
798
      if (iVertID == -1)
        {
          if (globs->Verbose) fprintf(stdout, " No level detected\n");
          checkLevel = false;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
799
      else if (globs->NumLevelRequest == 1 && IS_EQUAL(globs->LevelRequest[0], 0))
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
        {
          if (globs->Verbose) fprintf(stdout, " No level selected\n");
          globs->NumLevelRequest = 0;
          checkLevel = false;
        }
      else if (globs->Verbose)
        {
          if (Lhybrid2pressure)
            {
              if (globs->unitsel == 0)
                fprintf(stdout, " Selected pressure level:\n");
              else
                fprintf(stdout, " Selected height level:\n");
            }
          else
            {
              if (zaxisInqType(iVertID) == ZAXIS_HYBRID)
                fprintf(stdout, " Selected hybrid level:\n");
              else
                {
                  if (globs->unitsel == 0)
                    fprintf(stdout, " Selected pressure level:\n");
                  else
                    fprintf(stdout, " Selected height level:\n");
                }
            }
        }
    }

  if (globs->Verbose && iVertID != -1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
830
    for (l = 0; l < globs->NumLevelRequest; l++) fprintf(stdout, "  Level %2d = %13.4f\n", l + 1, globs->LevelRequest[l]);
831
832
833

  if (checkLevel)
    {
834
      for (k = 0; k < globs->NumLevelRequest; k++) removeLevel[k] = false;
835
836
837
838
      for (k = 0; k < globs->NumLevelRequest; k++)
        {
          level = globs->LevelRequest[k];
          for (l = k + 1; l < globs->NumLevelRequest; l++)
839
            if (removeLevel[l] == false && IS_EQUAL(level, globs->LevelRequest[l]))
840
              {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
841
                if (globs->Verbose) fprintf(stdout, "  Level %2d = %13.4f double request\n", l + 1, globs->LevelRequest[l]);
842
                removeLevel[l] = true;
843
844
              }
        }
845
846

      l = 0;
847
      for (k = 0; k < globs->NumLevelRequest; k++)
848
        if (removeLevel[k] == false) globs->LevelRequest[l++] = globs->LevelRequest[k];
849
850
851

      globs->NumLevelRequest = l;

852
853
      if (globs->AnalysisData || globs->Type < 30)
        {
854
          for (k = 0; k < globs->NumLevelRequest; k++) removeLevel[k] = false;
855
856
857
          for (k = 0; k < globs->NumLevelRequest; k++)
            {
              level = globs->LevelRequest[k];
858
              found = false;
859
              for (l = 0; l < globs->NumLevelFound; l++)
860
                if (IS_EQUAL(level, LevelFound[l])) found = true;
861
862
863

              if (!found)
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
864
                  fprintf(stdout, "  Level %2d = %14.4f not in input\n", k + 1, globs->LevelRequest[k]);
865
                  removeLevel[k] = true;
866
867
868
869
870
                }
            }

          l = 0;
          for (k = 0; k < globs->NumLevelRequest; k++)
871
            if (removeLevel[k] == false) globs->LevelRequest[l++] = globs->LevelRequest[k];
872
873
874
875

          if (l != globs->NumLevelRequest)
            {
              if (globs->Verbose) lprintf(stdout);
876
              afterAbort("Inconsistent or invalid level list!");
877
878
879
880
881
882
883
            }

          globs->NumLevelRequest = l;
        }
    }

  if (globs->Verbose) lprintf(stdout);
884
885
}

886
static void
887
after_defineLevel(AfterControl *globs, struct Variable *vars)
888
889
890
891
892
{
  int code, i;

  /* hybrid, pressure, height */

893
894
895
896
897
898
899
900
901
902
903
904
905
  switch (globs->Type)
    {
    case 0:
    case 10:
    case 11:
    case 20:
      {
        if (iVertID == -1) break;

        if (zaxisInqType(iVertID) == ZAXIS_HYBRID)
          {
            if (oVertID == -1)
              {
906
                if (globs->NumLevelRequest > globs->NumLevelFound) afterAbort("Too much level requested");
907
908
909
910

                if (globs->NumLevelFound == globs->NumLevelRequest)
                  {
                    for (i = 0; i < globs->NumLevelRequest; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
911
                      if (IS_NOT_EQUAL(globs->LevelRequest[i], LevelFound[i])) break;
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929

                    if (i == globs->NumLevelRequest) oVertID = iVertID;
                  }

                if (oVertID == -1 && globs->NumLevelRequest > 0)
                  {
                    oVertID = zaxisCreate(ZAXIS_HYBRID, globs->NumLevelRequest);
                    zaxisDefLevels(oVertID, globs->LevelRequest);
                    zaxisDefVct(oVertID, globs->nvct, globs->vct);
                  }
              }

            for (code = 0; code < MaxCodes; code++)
              {
                if (vars[code].selected)
                  {
                    if (vars[code].izaxisID != -1)
                      if (zaxisInqType(vars[code].izaxisID) == ZAXIS_HYBRID
Uwe Schulzweida's avatar
Uwe Schulzweida committed
930
                          && zaxisInqSize(vars[code].izaxisID) >= globs->NumLevelRequest)
931
932
933
934
935
936
937
                        vars[code].ozaxisID = oVertID;
                  }
              }
          }
        else
          {
            zaxisName(zaxisInqType(iVertID), zaxistypename);
938
            afterAbort("%s level data unsupported for TYPE %d", zaxistypename, globs->Type);
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
          }
        break;
      }
    case 30:
    case 40:
    case 41:
    case 50:
    case 60:
    case 61:
    case 70:
      {
        if (iVertID == -1) break;

        if (oVertID == -1)
          {
            if (globs->unitsel == 0)
              oVertID = zaxisCreate(ZAXIS_PRESSURE, globs->NumLevelRequest);
            else
              oVertID = zaxisCreate(ZAXIS_HEIGHT, globs->NumLevelRequest);

            zaxisDefLevels(oVertID, globs->LevelRequest);
          }

        for (code = 0; code < MaxCodes; code++)
          {
            if (vars[code].selected)
              {
                if (vars[code].izaxisID != -1)
                  {
                    int nlev = zaxisInqSize(vars[code].izaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
969
970
                    if (zaxisInqType(vars[code].izaxisID) == zaxisInqType(iVertID)
                        && (nlev == globs->NumLevel || nlev == globs->NumLevel + 1) && nlev > 1)
971
972
973
974
975
976
977
                      vars[code].ozaxisID = oVertID;
                  }
              }
          }

        break;
      }
978
    default: afterAbort("TYPE %d unsupported", globs->Type);
979
980
981
    }
}

982
static void
983
after_defineGrid(AfterControl *globs, struct Variable *vars)
984
985
986
987
{
  int ogridID = -1;
  int code;

988
  // spectral, fourier, gauss, zonal mean
989

990
991
992
993
994
995
996
  switch (globs->Type)
    {
    case 0:
    case 50:
      {
        if (specGridID == -1)
          {
997
998
            if (globs->DimSP == 0) afterAbort("dim spectral undefined");
            if (globs->Truncation == 0) afterAbort("truncation undefined");
999
1000
1001
1002

            specGridID = gridCreate(GRID_SPECTRAL, globs->DimSP);
            gridDefTrunc(specGridID, globs->Truncation);
          }
1003

1004
1005
1006
1007
1008
1009
1010
1011
1012
        ogridID = specGridID;
        break;
      }
    case 20:
    case 30:
    case 70:
      {
        if (gaussGridID == -1)
          {
1013
1014
            if (globs->Longitudes == 0) afterAbort("number of longitudes undefined");
            if (globs->Latitudes == 0) afterAbort("number of latitudes undefined");
1015

Uwe Schulzweida's avatar
Uwe Schulzweida committed
1016
            gaussGridID = gridCreate(GRID_GAUSSIAN, globs->Longitudes * globs->Latitudes);
1017
1018
            gridDefXsize(gaussGridID, globs->Longitudes);
            gridDefYsize(gaussGridID, globs->Latitudes);