Afterburner.cc

/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  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.
*/
/* ============================================================= */
/*                                                               */
/* 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                */
/*                                                               */
/* ============================================================= */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef _OPENMP
#include <omp.h>
#endif

#include <cdi.h>

#include "cdo_default_values.h"

#define streamOpenWrite cdoOpenWrite
#define streamDefVlist cdoDefVlist
#define streamDefTimestep cdoDefTimestep

#include "afterburner.h"
#include "constants.h"
#include "compare.h"
#include "cdo_options.h"
#include "mpmo_color.h"
#include "cdi_lockedIO.h"

int scan_par_obsolete(char *namelist, const char *name, int def);
int scan_par(int verbose, char *namelist, const char *name, int def);
void scan_code(char *namelist, struct Variable *vars, int maxCodes, int *numCodes);
int scan_time(int verbose, char *namelist, int *hours, int max_hours);
void scan_darray(char *namelist, const char *name, double *values, int maxValues, int *numValues);

char zaxistypename[CDI_MAX_NAME];

struct RARG
{
  int lana, nrecs;
  struct Variable *vars;
  AfterControl *globs;
};

cdo::Task *afterReadTask = nullptr;

static bool lstdout = true;

static int Source = 0;

static int ofiletype = -1;

static int DataType = -1;

static char *filename;
static const char **ifiles;
static char *ifile = nullptr;
static const char *ofile = nullptr;

static int ofileidx = 0;

static int specGridID = -1;
static int gaussGridID = -1;
static int iVertID = -1;
static int oVertID = -1;

static bool Lhybrid2pressure = false;

static int TsID;
bool afterReadAsync = true;

#define TIMESTEP_INTERVAL -1
#define MONTHLY_INTERVAL 0
#define DAILY_INTERVAL 1
#define UNLIM_INTERVAL 2

#define MaxHours 24
static int nrqh;
static int hours[MaxHours + 1];

static std::vector<double> LevelFound;

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;
}

static void
lprintf(FILE *fp)
{
  int num = 67;
  int cval = '-';

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

static void
FreeMean(struct Variable *vars)
{
  for (int code = 0; code < MaxCodes; code++)
    if (vars[code].mean)
      {
        Free(vars[code].mean);
        vars[code].mean = nullptr;
      }
}

static void
after_PostProcess(AfterControl *globs)
{
  if (globs->EndOfInterval)
    {
      if (lstdout)
        {
          if (globs->OutputInterval == DAILY_INTERVAL)
            fprintf(stdout, " Processed Day %2d  Month %2d  Year %04d", globs->OldDate.dy, globs->OldDate.mo, globs->OldDate.yr);
          else if (globs->OutputInterval == MONTHLY_INTERVAL)
            fprintf(stdout, " Processed Month %2d  Year %04d", globs->OldDate.mo, globs->OldDate.yr);
          else if (globs->OutputInterval == UNLIM_INTERVAL)
            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);

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

      globs->EndOfInterval = false;
      globs->MeanCount = 0;
    }
}

/* ================= */
/* switch input file */
/* ================= */
static void
after_SwitchFile(AfterControl *globs)
{
  bool echam4 = false;
  int n;
  char y3, y2, y1, y0;
  char m1, m0;
  char d1, d0;

  streamClose(globs->istreamID);

  if (globs->Multi > 0)
    {
      int i = strlen(ifile);
      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];
        }
      else
        {
          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;
        }
      else
        {
          ifile[i - 6] = y3;
          ifile[i - 5] = y2;
          ifile[i - 4] = y1;
          ifile[i - 3] = y0;
          ifile[i - 2] = m1;
          ifile[i - 1] = m0;
        }

      globs->Multi--;
    }

  if (globs->Nfiles > 0) ifile = (char *) ifiles[--globs->Nfiles];

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

  globs->istreamID = streamOpenReadLocked(ifile);

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

static int64_t
after_getDate(struct Date datetime)
{
  return cdiEncodeDate(datetime.yr, datetime.mo, datetime.dy);
}

static int
after_getTime(struct Date datetime)
{
  return cdiEncodeTime(datetime.hr, datetime.mn, 0);
}

static void
after_setDateTime(struct Date *datetime, int64_t date, int time)
{
  int sec;
  cdiDecodeDate(date, &datetime->yr, &datetime->mo, &datetime->dy);
  cdiDecodeTime(time, &datetime->hr, &datetime->mn, &sec);
}

static void
after_printProcessStatus(int tsID)
{
  static bool counthead = false;

  if (tsID == -1)
    {
      if (stdout_is_tty)
        {
          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");
          fflush(stdout);
        }

      counthead = false;
    }
  else
    {
      if (!counthead)
        {
          if (stdout_is_tty) fprintf(stdout, " Process timestep :       ");

          counthead = true;
        }

      if (stdout_is_tty)
        {
          fprintf(stdout, "\b\b\b\b\b\b%6d", tsID);
          fflush(stdout);
        }
    }
}

static int
after_setNextDate(AfterControl *globs)
{
  int nrecs = 0;

  bool righttime = false;
  while (true)
    {
      nrecs = streamInqTimestep(globs->istreamID, TsID);
      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;

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

      for (int i = 0; i < nrqh; i++)
        if (hours[i] < 0 || hours[i] == globs->NextDate.hr)
          {
            righttime = true;
            break;
          }

      if (righttime)
        break;
      else
        TsID += 1;
    }

  return nrecs;
}

static int num_recs = 0;

static void *
after_readTimestep(void *arg)
{
  int varID, gridID, zaxisID, levelID, timeID;
  size_t nmiss;
  auto rarg = (RARG *) arg;

  auto nrecs = rarg->nrecs;
  auto analysisData = rarg->lana;
  auto vars = rarg->vars;
  auto globs = rarg->globs;

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

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

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

      /* Skip records containing unneeded codes */

      if (!vars[code].needed0) continue;

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

      const auto leveltype = zaxisInqType(zaxisID);

      /* Skip records with unselected levels */

      int levelOffset = -1;
      /*
        if ( vars[code].ozaxisID != vars[code].izaxisID && ! Lhybrid2pressure )
      */
      if ((vars[code].ozaxisID != vars[code].izaxisID) && (leveltype == ZAXIS_PRESSURE))
        {
          const auto level = (int) zaxisInqLevel(zaxisID, levelID);
          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);

          fprintf(stderr, "  Code %3d   Level%6d   %6.4d-%2.2d-%2.2d  %2.2d:%2.2d:00\n", code,
                  (int) zaxisInqLevel(zaxisID, levelID), globs->OldDate.yr, globs->OldDate.mo, globs->OldDate.dy, globs->OldDate.hr,
                  globs->OldDate.mn);
        }

      if (analysisData)
        {
          streamReadRecord(globs->istreamID, globs->Field, &nmiss);
          after_AnalysisAddRecord(globs, vars, code, gridID, zaxisID, levelID, nmiss);
        }
      else
        {
          double *dataptr = after_get_dataptr(vars, code, gridID, zaxisID, levelID);
          streamReadRecord(globs->istreamID, dataptr, &nmiss);
          after_EchamAddRecord(globs, vars, code, gridID, zaxisID, levelID, nmiss);
        }

      if (iVertID != -1 && oVertID != -1 && (vars[code].izaxisID == iVertID)) vars[code].ozaxisID = oVertID;
    }

  TsID++;

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

  num_recs = after_setNextDate(globs);

  return (void *) &num_recs;
}

static void
after_defineNextTimestep(AfterControl *globs)
{
  static int otsID = 0;
  const auto vdate = after_getDate(globs->OldDate);
  const auto vtime = after_getTime(globs->OldDate);
  taxisDefVdate(globs->taxisID2, vdate);
  taxisDefVtime(globs->taxisID2, vtime);

  if (globs->Mean != 2)
    {
      if (otsID == 0)
        {
          const auto nvars = vlistNvars(globs->ovlistID);
          if (nvars == 0) afterAbort("No variable selected!");
          vlistDefTaxis(globs->ovlistID, globs->taxisID2);
          streamDefVlist(globs->ostreamID, globs->ovlistID);
        }
      taxisDefNumavg(globs->taxisID2, globs->MeanCount + 1);
      streamDefTimestep(globs->ostreamID, otsID);
    }

  otsID++;
}

static void
after_setEndOfInterval(AfterControl *globs, int nrecs)
{
  if (nrecs == 0)
    {
      globs->EndOfInterval = true;
    }
  else
    {
      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)
        globs->EndOfInterval = false;
      else
        afterAbort("output interval %d not implemented!", globs->OutputInterval);
    }
}

static void
after_moveTimestep(struct Variable *vars)
{
  int code;

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

  for (code = 0; code < MaxCodes; code++)
    if (vars[code].hybrid0)
      {
        vars[code].hybrid = vars[code].hybrid0;
        vars[code].hybrid0 = nullptr;
      }

  for (code = 0; code < MaxCodes; code++)
    if (vars[code].spectral0)
      {
        vars[code].spectral = vars[code].spectral0;
        vars[code].spectral0 = nullptr;
      }

  for (code = 0; code < MaxCodes; code++)
    if (vars[code].grid0)
      {
        vars[code].grid = vars[code].grid0;
        vars[code].grid0 = nullptr;
      }
}

static void
after_check_content(struct Variable *vars, int timestep)
{
  for (int code = 0; code < 272; code++)
    {
      /*  if ( code == GEOPOTENTIAL ) continue; */
      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)
        {
          if (vars[code].needed && !vars[code].selected && vars[code].spectral == nullptr && vars[code].hybrid == nullptr)
            {
              static bool lwarn = true;
              if (lwarn) cdoWarning("No humidity in data file, set to zero !");
              lwarn = false;
              vars[code].needed = false;
            }
        }
      else
        {
          if (vars[code].needed && !vars[code].comp && vars[code].spectral == nullptr && vars[code].hybrid == nullptr)
            {
              afterAbort("Code  %3d not found at timestep %d!", code, timestep);
            }
        }
    }
}

static void
after_control(AfterControl *globs, struct Variable *vars)
{
  int i;
  int nrecs = 0;
  int code;
  RARG rarg;

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

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

  TsID = 0;

  bool righttime = false;
  while (true)
    {
      nrecs = streamInqTimestep(globs->istreamID, TsID);
      if (nrecs <= 0) break;

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

      for (i = 0; i < nrqh; i++)
        if (hours[i] < 0 || hours[i] == globs->NewDate.hr)
          {
            righttime = true;
            break;
          }

      if (righttime) break;

      TsID++;
    }

  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))
    {
      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;

  bool tsFirst = true;

  while (nrecs > 0)
    {
      rarg.nrecs = nrecs;
      rarg.lana = globs->AnalysisData;
      rarg.vars = vars;
      rarg.globs = globs;

      int status = -1;

      if (tsFirst || !afterReadAsync)
        {
          if (afterReadAsync)
            {
              afterReadTask->start(after_readTimestep, &rarg);
            }
          else
            {
              status = *(int *) after_readTimestep(&rarg);
            }

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

          if (afterReadAsync)
            {
              status = *(int *) afterReadTask->wait();
              if (status < 0) cdoAbort("after_readTimestep error! (status = %d)", status);
            }
          tsFirst = false;
        }
      else
        {
          status = *(int *) afterReadTask->wait();
          if (status < 0) cdoAbort("after_readTimestep error! (status = %d)", status);
        }

      nrecs = status;

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

      after_moveTimestep(vars);

      if (nrecs && afterReadAsync)
        {
          afterReadTask->start(after_readTimestep, &rarg);
        }

      after_setEndOfInterval(globs, nrecs);

      if (lstdout) after_printProcessStatus(TsID);

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

      if (globs->Mean == 0 || globs->EndOfInterval)
        {
          if (!globs->AnalysisData) after_check_content(vars, globs->TermCount + 1);
          after_defineNextTimestep(globs);
        }

      if (globs->AnalysisData)
        after_processPL(globs, vars);
      else
        after_processML(globs, vars);

      after_PostProcess(globs);

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

      globs->OldDate = globs->NewDate;
    }

  if (afterReadTask) delete afterReadTask;
}

static void
after_setLevel(AfterControl *globs)
{
  int k, l, found;
  bool removeLevel[MaxLevel];
  double level;
  bool checkLevel = true;
  // default pressure level
  long plevelDefault[]
      = { 100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 25000, 20000, 15000, 10000, 7000, 5000, 3000, 2000, 1000 };
  // default height level
  long hlevelDefault[] = { 0, 1000, 2000, 5000, 10000, 15000, 20000, 25000, 30000 };

  const int numplevelDefault = sizeof(plevelDefault) / sizeof(plevelDefault[0]);
  const int numhlevelDefault = sizeof(hlevelDefault) / sizeof(hlevelDefault[0]);

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

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

  if (globs->NumLevelRequest == 0)
    {
      if (iVertID == -1)
        {
          if (globs->Verbose) fprintf(stdout, " No level detected\n");
        }
      else
        {
          if (Lhybrid2pressure)
            {
              if (globs->unitsel == 0)
                {
                  if (globs->Verbose) fprintf(stdout, " Default pressure level selected:\n");
                  globs->NumLevelRequest = numplevelDefault;
                  for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = plevelDefault[l];
                  oVertID = zaxisCreate(ZAXIS_PRESSURE, globs->NumLevelRequest);
                  zaxisDefLevels(oVertID, globs->LevelRequest);
                }
              else
                {
                  if (globs->Verbose) fprintf(stdout, " Default height level selected:\n");
                  globs->NumLevelRequest = numhlevelDefault;
                  for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = hlevelDefault[l];
                  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;
              for (l = 0; l < globs->NumLevelRequest; l++) globs->LevelRequest[l] = LevelFound[l];
              oVertID = iVertID;
            }
        }
      checkLevel = false;
    }
  else
    {
      if (iVertID == -1)
        {
          if (globs->Verbose) fprintf(stdout, " No level detected\n");
          checkLevel = false;
        }
      else if (globs->NumLevelRequest == 1 && IS_EQUAL(globs->LevelRequest[0], 0))
        {
          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)
    for (l = 0; l < globs->NumLevelRequest; l++) fprintf(stdout, "  Level %2d = %13.4f\n", l + 1, globs->LevelRequest[l]);

  if (checkLevel)
    {
      for (k = 0; k < globs->NumLevelRequest; k++) removeLevel[k] = false;
      for (k = 0; k < globs->NumLevelRequest; k++)
        {
          level = globs->LevelRequest[k];
          for (l = k + 1; l < globs->NumLevelRequest; l++)
            if (removeLevel[l] == false && IS_EQUAL(level, globs->LevelRequest[l]))
              {
                if (globs->Verbose) fprintf(stdout, "  Level %2d = %13.4f double request\n", l + 1, globs->LevelRequest[l]);
                removeLevel[l] = true;
              }
        }

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

      globs->NumLevelRequest = l;

      if (globs->AnalysisData || globs->Type < 30)
        {
          for (k = 0; k < globs->NumLevelRequest; k++) removeLevel[k] = false;
          for (k = 0; k < globs->NumLevelRequest; k++)
            {
              level = globs->LevelRequest[k];
              found = false;
              for (l = 0; l < globs->NumLevelFound; l++)
                if (IS_EQUAL(level, LevelFound[l])) found = true;

              if (!found)
                {
                  fprintf(stdout, "  Level %2d = %14.4f not in input\n", k + 1, globs->LevelRequest[k]);
                  removeLevel[k] = true;
                }
            }

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

          if (l != globs->NumLevelRequest)
            {
              if (globs->Verbose) lprintf(stdout);
              afterAbort("Inconsistent or invalid level list!");
            }

          globs->NumLevelRequest = l;
        }
    }

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

static void
after_defineLevel(AfterControl *globs, struct Variable *vars)
{
  int code, i;

  /* hybrid, pressure, height */

  switch (globs->Type)
    {
    case 0:
    case 10:
    case 11:
    case 20:
      {
        if (iVertID == -1) break;

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

                if (globs->NumLevelFound == globs->NumLevelRequest)
                  {
                    for (i = 0; i < globs->NumLevelRequest; i++)
                      if (IS_NOT_EQUAL(globs->LevelRequest[i], LevelFound[i])) break;

                    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
                          && zaxisInqSize(vars[code].izaxisID) >= globs->NumLevelRequest)
                        vars[code].ozaxisID = oVertID;
                  }
              }
          }
        else
          {
            zaxisName(zaxisInqType(iVertID), zaxistypename);
            afterAbort("%s level data unsupported for TYPE %d", zaxistypename, globs->Type);
          }
        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);
                    if (zaxisInqType(vars[code].izaxisID) == zaxisInqType(iVertID)
                        && (nlev == globs->NumLevel || nlev == globs->NumLevel + 1) && nlev > 1)
                      vars[code].ozaxisID = oVertID;
                  }
              }
          }

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

static void
after_defineGrid(AfterControl *globs, struct Variable *vars)
{
  int ogridID = -1;
  int code;

  // spectral, fourier, gauss, zonal mean

  switch (globs->Type)
    {
    case 0:
    case 50:
      {
        if (specGridID == -1)
          {
            if (globs->DimSP == 0) afterAbort("dim spectral undefined");
            if (globs->Truncation == 0) afterAbort("truncation undefined");

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

        ogridID = specGridID;
        break;
      }
    case 20:
    case 30:
    case 70:
      {
        if (gaussGridID == -1)
          {
            if (globs->Longitudes == 0) afterAbort("number of longitudes undefined");
            if (globs->Latitudes == 0) afterAbort("number of latitudes undefined");

            gaussGridID = gridCreate(GRID_GAUSSIAN, globs->Longitudes * globs->Latitudes);
            gridDefXsize(gaussGridID, globs->Longitudes);
            gridDefYsize(gaussGridID, globs->Latitudes);