mptrac.h

Go to the documentation of this file.

/*
  This file is part of MPTRAC.
  
  MPTRAC 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, either version 3 of the License, or
  (at your option) any later version.
  
  MPTRAC 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.
  
  You should have received a copy of the GNU General Public License
  along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
  
  Copyright (C) 2013-2024 Forschungszentrum Juelich GmbH
*/
 
#ifndef LIBTRAC_H
#define LIBTRAC_H
 
/* ------------------------------------------------------------
   Includes...
   ------------------------------------------------------------ */
 
#include <ctype.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_statistics.h>
#include <math.h>
#include <netcdf.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
 
#ifdef MPI
#include "mpi.h"
#endif
 
#ifdef _OPENACC
#include "openacc.h"
#endif
 
#ifdef CURAND
#include "curand.h"
#endif
 
#ifdef ZFP
#include "zfp.h"
#endif
 
#ifdef ZSTD
#include "zstd.h"
#endif
 
#ifdef CMS
#include "cmultiscale.h"
#endif
 
#ifdef KPP
#include "chem_Parameters.h"
#include "chem_Global.h"
#include "chem_Sparse.h"
#endif
 
/* ------------------------------------------------------------
   Constants...
   ------------------------------------------------------------ */
 
#ifndef AVO
#define AVO 6.02214076e23
#endif
 
#ifndef CPD
#define CPD 1003.5
#endif
 
#ifndef EPS
#define EPS (MH2O / MA)
#endif
 
#ifndef G0
#define G0 9.80665
#endif
 
#ifndef H0
#define H0 7.0
#endif
 
#ifndef LV
#define LV 2501000.
#endif
 
#ifndef KB
#define KB 1.3806504e-23
#endif
 
#ifndef MA
#define MA 28.9644
#endif
 
#ifndef MH2O
#define MH2O 18.01528
#endif
 
#ifndef MO3
#define MO3 48.00
#endif
 
#ifndef P0
#define P0 1013.25
#endif
 
#ifndef RA
#define RA (1e3 * RI / MA)
#endif
 
#ifndef RE
#define RE 6367.421
#endif
 
#ifndef RI
#define RI 8.3144598
#endif
 
#ifndef T0
#define T0 273.15
#endif
 
/* ------------------------------------------------------------
   Dimensions...
   ------------------------------------------------------------ */
 
#ifndef LEN
#define LEN 5000
#endif
 
#ifndef NP
#define NP 10000000
#endif
 
#ifndef NQ
#define NQ 15
#endif
 
#ifndef NCSI
#define NCSI 1000000
#endif
 
#ifndef EP
#define EP 140
#endif
 
#ifndef EX
#define EX 1202
#endif
 
#ifndef EY
#define EY 602
#endif
 
#ifndef NENS
#define NENS 2000
#endif
 
#ifndef NOBS
#define NOBS 10000000
#endif
 
#ifndef NTHREADS
#define NTHREADS 512
#endif
 
#ifndef CY
#define CY 250
#endif
 
#ifndef CO3
#define CO3 30
#endif
 
#ifndef CP
#define CP 70
#endif
 
#ifndef CSZA
#define CSZA 50
#endif
 
#ifndef CT
#define CT 12
#endif
 
#ifndef CTS
#define CTS 1000
#endif
 
/* ------------------------------------------------------------
   Macros...
   ------------------------------------------------------------ */
 
#ifdef _OPENACC
#define ALLOC(ptr, type, n)                             \
  if(acc_get_num_devices(acc_device_nvidia) <= 0)       \
    ERRMSG("Not running on a GPU device!");             \
  if((ptr=calloc((size_t)(n), sizeof(type)))==NULL)     \
    ERRMSG("Out of memory!");
#else
#define ALLOC(ptr, type, n)                              \
  if((ptr=calloc((size_t)(n), sizeof(type)))==NULL)      \
    ERRMSG("Out of memory!");
#endif
 
#define ARRAY_2D(ix, iy, ny)                    \
  ((ix) * (ny) + (iy))
 
#define ARRAY_3D(ix, iy, ny, iz, nz)            \
  (((ix)*(ny) + (iy)) * (nz) + (iz))
 
#define ARRHENIUS(a, b, t)                      \
  ((a) * exp( -(b) / (t)))
 
#define DEG2DX(dlon, lat)                                       \
  ((dlon) * M_PI * RE / 180. * cos((lat) / 180. * M_PI))
 
#define DEG2DY(dlat)                            \
  ((dlat) * M_PI * RE / 180.)
 
#define DP2DZ(dp, p)                            \
  (- (dp) * H0 / (p))
 
#define DX2DEG(dx, lat)                                         \
  (((lat) < -89.999 || (lat) > 89.999) ? 0                      \
   : (dx) * 180. / (M_PI * RE * cos((lat) / 180. * M_PI)))
 
#define DY2DEG(dy)                              \
  ((dy) * 180. / (M_PI * RE))
 
#define DZ2DP(dz, p)                            \
  (-(dz) * (p) / H0)
 
#define DIST(a, b)                              \
  sqrt(DIST2(a, b))
 
#define DIST2(a, b)                                                     \
  ((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])+(a[2]-b[2])*(a[2]-b[2]))
 
#define DOTP(a, b)                              \
  (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
 
#define FMOD(x, y)                              \
  ((x) - (int) ((x) / (y)) * (y))
 
#define FREAD(ptr, type, size, in) {                                    \
    if(fread(ptr, sizeof(type), size, in)!=size)                        \
      ERRMSG("Error while reading!");                                   \
  }
 
#define FWRITE(ptr, type, size, out) {                                  \
    if(fwrite(ptr, sizeof(type), size, out)!=size)                      \
      ERRMSG("Error while writing!");                                   \
  }
 
#define INTPOL_INIT                                             \
  double cw[4] = {0.0, 0.0, 0.0, 0.0}; int ci[3] = {0, 0, 0};
 
#define INTPOL_2D(var, init)                                            \
  intpol_met_time_2d(met0, met0->var, met1, met1->var,                  \
                     atm->time[ip], atm->lon[ip], atm->lat[ip],         \
                     &var, ci, cw, init);
 
#define INTPOL_3D(var, init)                                            \
  intpol_met_time_3d(met0, met0->var, met1, met1->var,                  \
                     atm->time[ip], atm->p[ip],                         \
                     atm->lon[ip], atm->lat[ip],                        \
                     &var, ci, cw, init);
 
#define INTPOL_SPACE_ALL(p, lon, lat) {                                 \
    intpol_met_space_3d(met, met->z, p, lon, lat, &z, ci, cw, 1);       \
    intpol_met_space_3d(met, met->t, p, lon, lat, &t, ci, cw, 0);       \
    intpol_met_space_3d(met, met->u, p, lon, lat, &u, ci, cw, 0);       \
    intpol_met_space_3d(met, met->v, p, lon, lat, &v, ci, cw, 0);       \
    intpol_met_space_3d(met, met->w, p, lon, lat, &w, ci, cw, 0);       \
    intpol_met_space_3d(met, met->pv, p, lon, lat, &pv, ci, cw, 0);     \
    intpol_met_space_3d(met, met->h2o, p, lon, lat, &h2o, ci, cw, 0);   \
    intpol_met_space_3d(met, met->o3, p, lon, lat, &o3, ci, cw, 0);     \
    intpol_met_space_3d(met, met->lwc, p, lon, lat, &lwc, ci, cw, 0);   \
    intpol_met_space_3d(met, met->rwc, p, lon, lat, &rwc, ci, cw, 0);   \
    intpol_met_space_3d(met, met->iwc, p, lon, lat, &iwc, ci, cw, 0);   \
    intpol_met_space_3d(met, met->swc, p, lon, lat, &swc, ci, cw, 0);   \
    intpol_met_space_3d(met, met->cc, p, lon, lat, &cc, ci, cw, 0);     \
    intpol_met_space_2d(met, met->ps, lon, lat, &ps, ci, cw, 0);        \
    intpol_met_space_2d(met, met->ts, lon, lat, &ts, ci, cw, 0);        \
    intpol_met_space_2d(met, met->zs, lon, lat, &zs, ci, cw, 0);        \
    intpol_met_space_2d(met, met->us, lon, lat, &us, ci, cw, 0);        \
    intpol_met_space_2d(met, met->vs, lon, lat, &vs, ci, cw, 0);        \
    intpol_met_space_2d(met, met->lsm, lon, lat, &lsm, ci, cw, 0);      \
    intpol_met_space_2d(met, met->sst, lon, lat, &sst, ci, cw, 0);      \
    intpol_met_space_2d(met, met->pbl, lon, lat, &pbl, ci, cw, 0);      \
    intpol_met_space_2d(met, met->pt, lon, lat, &pt, ci, cw, 0);        \
    intpol_met_space_2d(met, met->tt, lon, lat, &tt, ci, cw, 0);        \
    intpol_met_space_2d(met, met->zt, lon, lat, &zt, ci, cw, 0);        \
    intpol_met_space_2d(met, met->h2ot, lon, lat, &h2ot, ci, cw, 0);    \
    intpol_met_space_2d(met, met->pct, lon, lat, &pct, ci, cw, 0);      \
    intpol_met_space_2d(met, met->pcb, lon, lat, &pcb, ci, cw, 0);      \
    intpol_met_space_2d(met, met->cl, lon, lat, &cl, ci, cw, 0);        \
    intpol_met_space_2d(met, met->plcl, lon, lat, &plcl, ci, cw, 0);    \
    intpol_met_space_2d(met, met->plfc, lon, lat, &plfc, ci, cw, 0);    \
    intpol_met_space_2d(met, met->pel, lon, lat, &pel, ci, cw, 0);      \
    intpol_met_space_2d(met, met->cape, lon, lat, &cape, ci, cw, 0);    \
    intpol_met_space_2d(met, met->cin, lon, lat, &cin, ci, cw, 0);      \
    intpol_met_space_2d(met, met->o3c, lon, lat, &o3c, ci, cw, 0);      \
  }
 
#define INTPOL_TIME_ALL(time, p, lon, lat) {                            \
    intpol_met_time_3d(met0, met0->z, met1, met1->z, time, p, lon, lat, &z, ci, cw, 1); \
    intpol_met_time_3d(met0, met0->t, met1, met1->t, time, p, lon, lat, &t, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->u, met1, met1->u, time, p, lon, lat, &u, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->v, met1, met1->v, time, p, lon, lat, &v, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->w, met1, met1->w, time, p, lon, lat, &w, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->pv, met1, met1->pv, time, p, lon, lat, &pv, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->h2o, met1, met1->h2o, time, p, lon, lat, &h2o, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->o3, met1, met1->o3, time, p, lon, lat, &o3, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->lwc, met1, met1->lwc, time, p, lon, lat, &lwc, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->rwc, met1, met1->rwc, time, p, lon, lat, &rwc, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->iwc, met1, met1->iwc, time, p, lon, lat, &iwc, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->swc, met1, met1->swc, time, p, lon, lat, &swc, ci, cw, 0); \
    intpol_met_time_3d(met0, met0->cc, met1, met1->cc, time, p, lon, lat, &cc, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->ps, met1, met1->ps, time, lon, lat, &ps, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->ts, met1, met1->ts, time, lon, lat, &ts, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->zs, met1, met1->zs, time, lon, lat, &zs, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->us, met1, met1->us, time, lon, lat, &us, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->vs, met1, met1->vs, time, lon, lat, &vs, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->lsm, met1, met1->lsm, time, lon, lat, &lsm, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->sst, met1, met1->sst, time, lon, lat, &sst, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->pbl, met1, met1->pbl, time, lon, lat, &pbl, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->pt, met1, met1->pt, time, lon, lat, &pt, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->tt, met1, met1->tt, time, lon, lat, &tt, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->zt, met1, met1->zt, time, lon, lat, &zt, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->h2ot, met1, met1->h2ot, time, lon, lat, &h2ot, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->pct, met1, met1->pct, time, lon, lat, &pct, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->pcb, met1, met1->pcb, time, lon, lat, &pcb, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->cl, met1, met1->cl, time, lon, lat, &cl, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->plcl, met1, met1->plcl, time, lon, lat, &plcl, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->plfc, met1, met1->plfc, time, lon, lat, &plfc, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->pel, met1, met1->pel, time, lon, lat, &pel, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->cape, met1, met1->cape, time, lon, lat, &cape, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->cin, met1, met1->cin, time, lon, lat, &cin, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->o3c, met1, met1->o3c, time, lon, lat, &o3c, ci, cw, 0); \
  }
 
#define LAPSE(p1, t1, p2, t2)                                           \
  (1e3 * G0 / RA * ((t2) - (t1)) / ((t2) + (t1))                        \
   * ((p2) + (p1)) / ((p2) - (p1)))
 
#define LIN(x0, y0, x1, y1, x)                  \
  ((y0)+((y1)-(y0))/((x1)-(x0))*((x)-(x0)))
 
#define MAX(a,b)                                \
  (((a)>(b))?(a):(b))
 
#define MET_HEADER                                                      \
  fprintf(out,                                                          \
          "# $1 = time [s]\n"                                           \
          "# $2 = altitude [km]\n"                                      \
          "# $3 = longitude [deg]\n"                                    \
          "# $4 = latitude [deg]\n"                                     \
          "# $5 = pressure [hPa]\n"                                     \
          "# $6 = temperature [K]\n"                                    \
          "# $7 = zonal wind [m/s]\n"                                   \
          "# $8 = meridional wind [m/s]\n"                              \
          "# $9 = vertical velocity [hPa/s]\n"                          \
          "# $10 = H2O volume mixing ratio [ppv]\n");                   \
  fprintf(out,                                                          \
          "# $11 = O3 volume mixing ratio [ppv]\n"                      \
          "# $12 = geopotential height [km]\n"                          \
          "# $13 = potential vorticity [PVU]\n"                         \
          "# $14 = surface pressure [hPa]\n"                            \
          "# $15 = surface temperature [K]\n"                           \
          "# $16 = surface geopotential height [km]\n"                  \
          "# $17 = surface zonal wind [m/s]\n"                          \
          "# $18 = surface meridional wind [m/s]\n"                     \
          "# $19 = land-sea mask [1]\n"                                 \
          "# $20 = sea surface temperature [K]\n");                     \
  fprintf(out,                                                          \
          "# $21 = tropopause pressure [hPa]\n"                         \
          "# $22 = tropopause geopotential height [km]\n"               \
          "# $23 = tropopause temperature [K]\n"                        \
          "# $24 = tropopause water vapor [ppv]\n"                      \
          "# $25 = cloud liquid water content [kg/kg]\n"                \
          "# $26 = cloud rain water content [kg/kg]\n"                  \
          "# $27 = cloud ice water content [kg/kg]\n"                   \
          "# $28 = cloud snow water content [kg/kg]\n"                  \
          "# $29 = cloud cover [1]\n"                                   \
          "# $30 = total column cloud water [kg/m^2]\n");               \
  fprintf(out,                                                          \
          "# $31 = cloud top pressure [hPa]\n"                          \
          "# $32 = cloud bottom pressure [hPa]\n"                       \
          "# $33 = pressure at lifted condensation level (LCL) [hPa]\n" \
          "# $34 = pressure at level of free convection (LFC) [hPa]\n"  \
          "# $35 = pressure at equilibrium level (EL) [hPa]\n"          \
          "# $36 = convective available potential energy (CAPE) [J/kg]\n" \
          "# $37 = convective inhibition (CIN) [J/kg]\n"                \
          "# $38 = relative humidity over water [%%]\n"                 \
          "# $39 = relative humidity over ice [%%]\n"                   \
          "# $40 = dew point temperature [K]\n");                       \
  fprintf(out,                                                          \
          "# $41 = frost point temperature [K]\n"                       \
          "# $42 = NAT temperature [K]\n"                               \
          "# $43 = HNO3 volume mixing ratio [ppv]\n"                    \
          "# $44 = OH volume mixing ratio [ppv]\n"                      \
          "# $45 = H2O2 volume mixing ratio [ppv]\n"                    \
          "# $46 = HO2 volume mixing ratio [ppv]\n"                     \
          "# $47 = O(1D) volume mixing ratio [ppv]\n"                   \
          "# $48 = boundary layer pressure [hPa]\n"                     \
          "# $49 = total column ozone [DU]\n"                           \
          "# $50 = number of data points\n");                           \
  fprintf(out,                                                          \
          "# $51 = number of tropopause data points\n"                  \
          "# $52 = number of CAPE data points\n");
 
#define MIN(a,b)                                \
  (((a)<(b))?(a):(b))
 
#define MOLEC_DENS(p,t)                 \
  (AVO * 1e-6 * ((p) * 100) / (RI * (t)))
 
#define NC(cmd) {                                    \
    int nc_result=(cmd);                             \
    if(nc_result!=NC_NOERR)                          \
      ERRMSG("%s", nc_strerror(nc_result));          \
  }
 
#define NC_DEF_VAR(varname, type, ndims, dims, long_name, units) {      \
    NC(nc_def_var(ncid, varname, type, ndims, dims, &varid));           \
    NC(nc_put_att_text(ncid, varid, "long_name", strnlen(long_name, LEN), long_name)); \
    NC(nc_put_att_text(ncid, varid, "units", strnlen(units, LEN), units)); \
  }
 
#define NC_GET_DOUBLE(varname, ptr, force) {                    \
    if(force) {                                                 \
      NC(nc_inq_varid(ncid, varname, &varid));                  \
      NC(nc_get_var_double(ncid, varid, ptr));                  \
    } else {                                                    \
      if(nc_inq_varid(ncid, varname, &varid) == NC_NOERR) {     \
        NC(nc_get_var_double(ncid, varid, ptr));                \
      } else                                                    \
        WARN("netCDF variable %s is missing!", varname);        \
    }                                                           \
  }
 
#define NC_INQ_DIM(dimname, ptr, min, max) {            \
    int dimid; size_t naux;                             \
    NC(nc_inq_dimid(ncid, dimname, &dimid));            \
    NC(nc_inq_dimlen(ncid, dimid, &naux));              \
    *ptr = (int)naux;                                   \
    if ((*ptr) < (min) || (*ptr) > (max))               \
      ERRMSG("Dimension %s is out of range!", dimname); \
  }
 
#define NC_PUT_DOUBLE(varname, ptr, hyperslab) {                \
    NC(nc_inq_varid(ncid, varname, &varid));                    \
    if(hyperslab) {                                             \
      NC(nc_put_vara_double(ncid, varid, start, count, ptr));   \
    } else {                                                    \
      NC(nc_put_var_double(ncid, varid, ptr));                  \
    }                                                           \
  }
 
#define NC_PUT_FLOAT(varname, ptr, hyperslab) {                 \
    NC(nc_inq_varid(ncid, varname, &varid));                    \
    if(hyperslab) {                                             \
      NC(nc_put_vara_float(ncid, varid, start, count, ptr));    \
    } else {                                                    \
      NC(nc_put_var_float(ncid, varid, ptr));                   \
    }                                                           \
  }
 
#define NC_PUT_INT(varname, ptr, hyperslab) {                   \
    NC(nc_inq_varid(ncid, varname, &varid));                    \
    if(hyperslab) {                                             \
      NC(nc_put_vara_int(ncid, varid, start, count, ptr));      \
    } else {                                                    \
      NC(nc_put_var_int(ncid, varid, ptr));                     \
    }                                                           \
  }
 
#define NC_PUT_ATT(varname, attname, text) {                            \
    NC(nc_inq_varid(ncid, varname, &varid));                            \
    NC(nc_put_att_text(ncid, varid, attname, strnlen(text, LEN), text)); \
  }
 
#define NC_PUT_ATT_GLOBAL(attname, text)                                \
  NC(nc_put_att_text(ncid, NC_GLOBAL, attname, strnlen(text, LEN), text));
 
#define NN(x0, y0, x1, y1, x)                           \
  (fabs((x) - (x0)) <= fabs((x) - (x1)) ? (y0) : (y1))
 
#define NORM(a)                                 \
  sqrt(DOTP(a, a))
 
#ifdef _OPENACC
#define PARTICLE_LOOP(ip0, ip1, check_dt, ...)          \
  const int ip0_const = ip0;                            \
  const int ip1_const = ip1;                            \
  _Pragma(__VA_ARGS__)                                  \
  _Pragma("acc parallel loop independent gang vector")  \
  for (int ip = ip0_const; ip < ip1_const; ip++)        \
    if (!check_dt || dt[ip] != 0)
#else
#define PARTICLE_LOOP(ip0, ip1, check_dt, ...)          \
  const int ip0_const = ip0;                            \
  const int ip1_const = ip1;                            \
  _Pragma("omp parallel for default(shared)")           \
  for (int ip = ip0_const; ip < ip1_const; ip++)        \
    if (!check_dt || dt[ip] != 0)
#endif
 
#define P(z)                                    \
  (P0 * exp(-(z) / H0))
 
#define PSAT(t)                                                 \
  (6.112 * exp(17.62 * ((t) - T0) / (243.12 + (t) - T0)))
 
#define PSICE(t)                                                \
  (6.112 * exp(22.46 * ((t) - T0) / (272.62 + (t) - T0)))
 
#define PW(p, h2o)                                                      \
  ((p) * MAX((h2o), 0.1e-6) / (1. + (1. - EPS) * MAX((h2o), 0.1e-6)))
 
#define RH(p, t, h2o)                           \
  (PW(p, h2o) / PSAT(t) * 100.)
 
#define RHICE(p, t, h2o)                        \
  (PW(p, h2o) / PSICE(t) * 100.)
 
#define RHO(p, t)                               \
  (100. * (p) / (RA * (t)))
 
#define ROETH_PHOTOL(a, b, c, sza)                                      \
  ((c)*(sza) < M_PI/2. ? (a) * exp((b) * (1 - 1/cos((c) * (sza)))) : 0)
 
#define SET_ATM(qnt, val)                       \
  if (ctl->qnt >= 0)                            \
    atm->q[ctl->qnt][ip] = val;
 
#define SET_QNT(qnt, name, longname, unit)              \
  if (strcasecmp(ctl->qnt_name[iq], name) == 0) {       \
    ctl->qnt = iq;                                      \
    sprintf(ctl->qnt_longname[iq], longname);           \
    sprintf(ctl->qnt_unit[iq], unit);                   \
  } else
 
#define SH(h2o)                                 \
  (EPS * MAX((h2o), 0.1e-6))
 
#define SQR(x)                                  \
  ((x)*(x))
 
#define SWAP(x, y, type)                        \
  do {type tmp = x; x = y; y = tmp;} while(0);
 
#define TDEW(p, h2o)                            \
  (T0 + 243.12 * log(PW((p), (h2o)) / 6.112)    \
   / (17.62 - log(PW((p), (h2o)) / 6.112)))
 
#define TICE(p, h2o)                            \
  (T0 + 272.62 * log(PW((p), (h2o)) / 6.112)    \
   / (22.46 - log(PW((p), (h2o)) / 6.112)))
 
#define THETA(p, t)                             \
  ((t) * pow(1000. / (p), 0.286))
 
#define THETAVIRT(p, t, h2o)                    \
  (TVIRT(THETA((p), (t)), MAX((h2o), 0.1e-6)))
 
#define TOK(line, tok, format, var) {                                   \
    if(((tok)=strtok((line), " \t"))) {                                 \
      if(sscanf(tok, format, &(var))!=1) continue;                      \
    } else ERRMSG("Error while reading!");                              \
  }
 
#define TVIRT(t, h2o)                                   \
  ((t) * (1. + (1. - EPS) * MAX((h2o), 0.1e-6)))
 
#define Z(p)                                    \
  (H0 * log(P0 / (p)))
 
#define ZDIFF(lnp0, t0, h2o0, lnp1, t1, h2o1)                           \
  (RI / MA / G0 * 0.5 * (TVIRT((t0), (h2o0)) + TVIRT((t1), (h2o1)))     \
   * ((lnp0) - (lnp1)))
 
#define ZETA(ps, p, t)                                                  \
  (((p) / (ps) <= 0.3 ? 1. :                                            \
    sin(M_PI / 2. * (1. - (p) / (ps)) / (1. - 0.3)))                    \
   * THETA((p), (t)))
 
/* ------------------------------------------------------------
   Log messages...
   ------------------------------------------------------------ */
 
#ifndef LOGLEV
#define LOGLEV 2
#endif
 
#define LOG(level, ...) {                                               \
    if(level >= 2)                                                      \
      printf("  ");                                                     \
    if(level <= LOGLEV) {                                               \
      printf(__VA_ARGS__);                                              \
      printf("\n");                                                     \
    }                                                                   \
  }
 
#define WARN(...) {                                                     \
    printf("\nWarning (%s, %s, l%d): ", __FILE__, __func__, __LINE__);  \
    LOG(0, __VA_ARGS__);                                                \
  }
 
#define ERRMSG(...) {                                                   \
    printf("\nError (%s, %s, l%d): ", __FILE__, __func__, __LINE__);    \
    LOG(0, __VA_ARGS__);                                                \
    exit(EXIT_FAILURE);                                                 \
  }
 
#define PRINT(format, var)                                              \
  printf("Print (%s, %s, l%d): %s= "format"\n",                         \
         __FILE__, __func__, __LINE__, #var, var);
 
/* ------------------------------------------------------------
   Timers...
   ------------------------------------------------------------ */
 
#define NTIMER 100
 
#define PRINT_TIMERS                            \
  timer("END", "END", 1);
 
#define SELECT_TIMER(id, group, color) {                                \
    NVTX_POP;                                                           \
    NVTX_PUSH(id, color);                                               \
    timer(id, group, 0);                                                \
  }
 
#define START_TIMERS                            \
  NVTX_PUSH("START", NVTX_CPU);
 
#define STOP_TIMERS                             \
  NVTX_POP;
 
/* ------------------------------------------------------------
   NVIDIA Tools Extension (NVTX)...
   ------------------------------------------------------------ */
 
#ifdef NVTX
#include "nvToolsExt.h"
 
#define NVTX_CPU 0xFFADD8E6
 
#define NVTX_GPU 0xFF00008B
 
#define NVTX_H2D 0xFFFFFF00
 
#define NVTX_D2H 0xFFFF8800
 
#define NVTX_READ 0xFFFFCCCB
 
#define NVTX_WRITE 0xFF8B0000
 
#define NVTX_RECV 0xFFCCFFCB
 
#define NVTX_SEND 0xFF008B00
 
#define NVTX_PUSH(range_title, range_color) {           \
    nvtxEventAttributes_t eventAttrib = {0};            \
    eventAttrib.version = NVTX_VERSION;                 \
    eventAttrib.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;   \
    eventAttrib.messageType = NVTX_MESSAGE_TYPE_ASCII;  \
    eventAttrib.colorType = NVTX_COLOR_ARGB;            \
    eventAttrib.color = range_color;                    \
    eventAttrib.message.ascii = range_title;            \
    nvtxRangePushEx(&eventAttrib);                      \
  }
 
#define NVTX_POP {                              \
    nvtxRangePop();                             \
  }
#else
 
/* Empty definitions of NVTX_PUSH and NVTX_POP... */
#define NVTX_PUSH(range_title, range_color) {}
#define NVTX_POP {}
#endif
 
/* ------------------------------------------------------------
   Thrust...
   ------------------------------------------------------------ */
 
void thrustSortWrapper(
  double *__restrict__ c,
  int n,
  int *__restrict__ index);
 
/* ------------------------------------------------------------
   Structs...
   ------------------------------------------------------------ */
 
typedef struct {
 
  /* TODO: finally sort ctl parameters once Fortran wrapper is working! */
 
  int advect_cpl_zeta_and_press_modules;
 
  int met_press_level_def;
 
  int advect_vert_coord;
 
  int met_vert_coord;
 
  int met_clams;
 
  /* ------------------------------------------------------------
     Quantity parameters...
     ------------------------------------------------------------ */
 
  int nq;
 
  char qnt_name[NQ][LEN];
 
  char qnt_longname[NQ][LEN];
 
  char qnt_unit[NQ][LEN];
 
  char qnt_format[NQ][LEN];
 
  int qnt_idx;
 
  int qnt_ens;
 
  int qnt_stat;
 
  int qnt_m;
 
  int qnt_vmr;
 
  int qnt_rp;
 
  int qnt_rhop;
 
  int qnt_ps;
 
  int qnt_ts;
 
  int qnt_zs;
 
  int qnt_us;
 
  int qnt_vs;
 
  int qnt_lsm;
 
  int qnt_sst;
 
  int qnt_pbl;
 
  int qnt_pt;
 
  int qnt_tt;
 
  int qnt_zt;
 
  int qnt_h2ot;
 
  int qnt_zg;
 
  int qnt_p;
 
  int qnt_t;
 
  int qnt_rho;
 
  int qnt_u;
 
  int qnt_v;
 
  int qnt_w;
 
  int qnt_h2o;
 
  int qnt_o3;
 
  int qnt_lwc;
 
  int qnt_rwc;
 
  int qnt_iwc;
 
  int qnt_swc;
 
  int qnt_cc;
 
  int qnt_pct;
 
  int qnt_pcb;
 
  int qnt_cl;
 
  int qnt_plcl;
 
  int qnt_plfc;
 
  int qnt_pel;
 
  int qnt_cape;
 
  int qnt_cin;
 
  int qnt_o3c;
 
  int qnt_hno3;
 
  int qnt_oh;
 
  int qnt_h2o2;
 
  int qnt_ho2;
 
  int qnt_o1d;
 
  int qnt_mloss_oh;
 
  int qnt_mloss_h2o2;
 
  int qnt_mloss_kpp;
 
  int qnt_mloss_wet;
 
  int qnt_mloss_dry;
 
  int qnt_mloss_decay;
 
  int qnt_loss_rate;
 
  int qnt_psat;
 
  int qnt_psice;
 
  int qnt_pw;
 
  int qnt_sh;
 
  int qnt_rh;
 
  int qnt_rhice;
 
  int qnt_theta;
 
  int qnt_zeta;
 
  int qnt_zeta_d;
 
  int qnt_tvirt;
 
  int qnt_lapse;
 
  int qnt_vh;
 
  int qnt_vz;
 
  int qnt_pv;
 
  int qnt_tdew;
 
  int qnt_tice;
 
  int qnt_tsts;
 
  int qnt_tnat;
 
  int qnt_Cx;
 
  int qnt_Ch2o;
 
  int qnt_Co3;
 
  int qnt_Cco;
 
  int qnt_Coh;
 
  int qnt_Ch;
 
  int qnt_Cho2;
 
  int qnt_Ch2o2;
 
  int qnt_Co1d;
 
  int qnt_Co3p;
 
  int qnt_Cccl4;
 
  int qnt_Cccl3f;
 
  int qnt_Cccl2f2;
 
  int qnt_Cn2o;
 
  int qnt_Csf6;
 
  int qnt_aoa;
 
  int direction;
 
  double t_start;
 
  double t_stop;
 
  double dt_mod;
 
  /* ------------------------------------------------------------
     Meteo data parameters...
     ------------------------------------------------------------ */
 
  char metbase[LEN];
 
  double dt_met;
 
  int met_convention;
 
  int met_type;
 
  int met_nc_scale;
 
  int met_zfp_prec;
 
  double met_zfp_tol_t;
 
  double met_zfp_tol_z;
 
  int met_cms_heur;
 
  double met_cms_eps_z;
 
  double met_cms_eps_t;
 
  double met_cms_eps_u;
 
  double met_cms_eps_v;
 
  double met_cms_eps_w;
 
  double met_cms_eps_pv;
 
  double met_cms_eps_h2o;
 
  double met_cms_eps_o3;
 
  double met_cms_eps_lwc;
 
  double met_cms_eps_rwc;
 
  double met_cms_eps_iwc;
 
  double met_cms_eps_swc;
 
  double met_cms_eps_cc;
 
  int met_dx;
 
  int met_dy;
 
  int met_dp;
 
  int met_sx;
 
  int met_sy;
 
  int met_sp;
 
  double met_detrend;
 
  int met_np;
 
  double met_p[EP];
 
  int met_geopot_sx;
 
  int met_geopot_sy;
 
  int met_relhum;
 
  int met_tropo;
 
  double met_tropo_pv;
 
  double met_tropo_theta;
 
  int met_tropo_spline;
 
  double met_dt_out;
 
  int met_cache;
 
  int met_mpi_share;
 
  /* ------------------------------------------------------------
     Geophysical module parameters...
     ------------------------------------------------------------ */
 
  double sort_dt;
 
  int isosurf;
 
  char balloon[LEN];
 
  int advect;
 
  int reflect;
 
  int rng_type;
 
  double turb_dx_trop;
 
  double turb_dx_strat;
 
  double turb_dz_trop;
 
  double turb_dz_strat;
 
  double turb_mesox;
 
  double turb_mesoz;
 
  double conv_cape;
 
  double conv_cin;
 
  double conv_dt;
 
  int conv_mix_bot;
 
  int conv_mix_top;
 
  double bound_mass;
 
  double bound_mass_trend;
 
  double bound_vmr;
 
  double bound_vmr_trend;
 
  double bound_lat0;
 
  double bound_lat1;
 
  double bound_p0;
 
  double bound_p1;
 
  double bound_dps;
 
  double bound_dzs;
 
  double bound_zetas;
 
  int bound_pbl;
 
  char species[LEN];
 
  double molmass;
 
  double tdec_trop;
 
  double tdec_strat;
 
  char clim_photo[LEN];
 
  char clim_hno3_filename[LEN];
 
  char clim_oh_filename[LEN];
 
  char clim_h2o2_filename[LEN];
 
  char clim_ho2_filename[LEN];
 
  char clim_o1d_filename[LEN];
 
  char clim_o3_filename[LEN];
 
  char clim_ccl4_timeseries[LEN];
 
  char clim_ccl3f_timeseries[LEN];
 
  char clim_ccl2f2_timeseries[LEN];
 
  char clim_n2o_timeseries[LEN];
 
  char clim_sf6_timeseries[LEN];
 
  double mixing_dt;
 
  double mixing_trop;
 
  double mixing_strat;
 
  int mixing_nz;
 
  double mixing_z0;
 
  double mixing_z1;
 
  int mixing_nx;
 
  double mixing_lon0;
 
  double mixing_lon1;
 
  int mixing_ny;
 
  double mixing_lat0;
 
  double mixing_lat1;
 
  int chemgrid_nz;
 
  double chemgrid_z0;
 
  double chemgrid_z1;
 
  int chemgrid_nx;
 
  double chemgrid_lon0;
 
  double chemgrid_lon1;
 
  int chemgrid_ny;
 
  double chemgrid_lat0;
 
  double chemgrid_lat1;
 
  int oh_chem_reaction;
 
  double oh_chem[4];
 
  double oh_chem_beta;
 
  int h2o2_chem_reaction;
 
  int kpp_chem;
 
  double dt_kpp;
 
  int tracer_chem;
 
  double wet_depo_pre[2];
 
  double wet_depo_bc_a;
 
  double wet_depo_bc_b;
 
  double wet_depo_ic_a;
 
  double wet_depo_ic_b;
 
  double wet_depo_ic_h[3];
 
  double wet_depo_bc_h[2];
 
  double wet_depo_ic_ret_ratio;
 
  double wet_depo_bc_ret_ratio;
 
  double dry_depo_dp;
 
  double dry_depo_vdep;
 
  double psc_h2o;
 
  double psc_hno3;
 
  /* ------------------------------------------------------------
     Output parameters...
     ------------------------------------------------------------ */
 
  char atm_basename[LEN];
 
  char atm_gpfile[LEN];
 
  double atm_dt_out;
 
  int atm_filter;
 
  int atm_stride;
 
  int atm_type;
 
  int atm_type_out;
 
  int obs_type;
 
  char csi_basename[LEN];
 
  char csi_kernel[LEN];
 
  double csi_dt_out;
 
  char csi_obsfile[LEN];
 
  double csi_obsmin;
 
  double csi_modmin;
 
  int csi_nz;
 
  double csi_z0;
 
  double csi_z1;
 
  int csi_nx;
 
  double csi_lon0;
 
  double csi_lon1;
 
  int csi_ny;
 
  double csi_lat0;
 
  double csi_lat1;
 
  char ens_basename[LEN];
 
  double ens_dt_out;
 
  char grid_basename[LEN];
 
  char grid_kernel[LEN];
 
  char grid_gpfile[LEN];
 
  double grid_dt_out;
 
  int grid_sparse;
 
  int grid_stddev;
 
  int grid_nz;
 
  double grid_z0;
 
  double grid_z1;
 
  int grid_nx;
 
  double grid_lon0;
 
  double grid_lon1;
 
  int grid_ny;
 
  double grid_lat0;
 
  double grid_lat1;
 
  int grid_type;
 
  char prof_basename[LEN];
 
  char prof_obsfile[LEN];
 
  int prof_nz;
 
  double prof_z0;
 
  double prof_z1;
 
  int prof_nx;
 
  double prof_lon0;
 
  double prof_lon1;
 
  int prof_ny;
 
  double prof_lat0;
 
  double prof_lat1;
 
  char sample_basename[LEN];
 
  char sample_kernel[LEN];
 
  char sample_obsfile[LEN];
 
  double sample_dx;
 
  double sample_dz;
 
  char stat_basename[LEN];
 
  double stat_lon;
 
  double stat_lat;
 
  double stat_r;
 
  double stat_t0;
 
  double stat_t1;
 
  char vtk_basename[LEN];
 
  double vtk_dt_out;
 
  int vtk_stride;
 
  double vtk_scale;
 
  double vtk_offset;
 
  int vtk_sphere;
 
} ctl_t;
 
typedef struct {
 
  int np;
 
  double time[NP];
 
  double p[NP];
 
  double lon[NP];
 
  double lat[NP];
 
  double q[NQ][NP];
 
} atm_t;
 
typedef struct {
 
  double iso_var[NP];
 
  double iso_ps[NP];
 
  double iso_ts[NP];
 
  int iso_n;
 
  float uvwp[NP][3];
 
} cache_t;
 
typedef struct {
 
  int np;
 
  int nsza;
 
  int no3c;
 
  double p[CP];
 
  double sza[CSZA];
 
  double o3c[CO3];
 
  double n2o[CP][CSZA][CO3];
 
  double ccl4[CP][CSZA][CO3];
 
  double ccl3f[CP][CSZA][CO3];
 
  double ccl2f2[CP][CSZA][CO3];
 
  double o2[CP][CSZA][CO3];
 
  double o3_1[CP][CSZA][CO3];
 
  double o3_2[CP][CSZA][CO3];
 
  double h2o2[CP][CSZA][CO3];
 
  double h2o[CP][CSZA][CO3];
 
} clim_photo_t;
 
typedef struct {
 
  int ntime;
 
  double time[CTS];
 
  double vmr[CTS];
 
} clim_ts_t;
 
typedef struct {
 
  int ntime;
 
  int nlat;
 
  int np;
 
  double time[CT];
 
  double lat[CY];
 
  double p[CP];
 
  double vmr[CT][CP][CY];
 
} clim_zm_t;
 
typedef struct {
 
  int tropo_ntime;
 
  int tropo_nlat;
 
  double tropo_time[12];
 
  double tropo_lat[73];
 
  double tropo[12][73];
 
  clim_photo_t photo;
 
  clim_zm_t hno3;
 
  clim_zm_t oh;
 
  clim_zm_t h2o2;
 
  clim_zm_t ho2;
 
  clim_zm_t o1d;
 
  clim_ts_t ccl4;
 
  clim_ts_t ccl3f;
 
  clim_ts_t ccl2f2;
 
  clim_ts_t n2o;
 
  clim_ts_t sf6;
 
} clim_t;
 
typedef struct {
 
  double time;
 
  int nx;
 
  int ny;
 
  int np;
 
  int npl;
 
  double lon[EX];
 
  double lat[EY];
 
  double p[EP];
 
  double hybrid[EP];
 
  float ps[EX][EY];
 
  float ts[EX][EY];
 
  float zs[EX][EY];
 
  float us[EX][EY];
 
  float vs[EX][EY];
 
  float lsm[EX][EY];
 
  float sst[EX][EY];
 
  float pbl[EX][EY];
 
  float pt[EX][EY];
 
  float tt[EX][EY];
 
  float zt[EX][EY];
 
  float h2ot[EX][EY];
 
  float pct[EX][EY];
 
  float pcb[EX][EY];
 
  float cl[EX][EY];
 
  float plcl[EX][EY];
 
  float plfc[EX][EY];
 
  float pel[EX][EY];
 
  float cape[EX][EY];
 
  float cin[EX][EY];
 
  float o3c[EX][EY];
 
  float z[EX][EY][EP];
 
  float t[EX][EY][EP];
 
  float u[EX][EY][EP];
 
  float v[EX][EY][EP];
 
  float w[EX][EY][EP];
 
  float pv[EX][EY][EP];
 
  float h2o[EX][EY][EP];
 
  float o3[EX][EY][EP];
 
  float lwc[EX][EY][EP];
 
  float rwc[EX][EY][EP];
 
  float iwc[EX][EY][EP];
 
  float swc[EX][EY][EP];
 
  float cc[EX][EY][EP];
 
  float pl[EX][EY][EP];
 
  float ul[EX][EY][EP];
 
  float vl[EX][EY][EP];
 
  float wl[EX][EY][EP];
 
  float zetal[EX][EY][EP];
 
  float zeta_dotl[EX][EY][EP];
 
} met_t;
 
/* ------------------------------------------------------------
   Functions...
   ------------------------------------------------------------ */
 
void broadcast_large_data(
  void *data,
  size_t N);
 
void cart2geo(
  const double *x,
  double *z,
  double *lon,
  double *lat);
 
double clim_oh(
  const ctl_t * ctl,
  const clim_t * clim,
  const double t,
  const double lon,
  const double lat,
  const double p);
 
void clim_oh_diurnal_correction(
  ctl_t * ctl,
  clim_t * clim);
 
double clim_photo(
  double rate[CP][CSZA][CO3],
  clim_photo_t * photo,
  double p,
  double sza,
  double o3c);
 
double clim_tropo(
  const clim_t * clim,
  const double t,
  const double lat);
 
void clim_tropo_init(
  clim_t * clim);
 
double clim_ts(
  const clim_ts_t * ts,
  const double t);
 
double clim_zm(
  const clim_zm_t * zm,
  const double t,
  const double lat,
  const double p);
 
void compress_cms(
  ctl_t * ctl,
  char *varname,
  float *array,
  size_t nx,
  size_t ny,
  size_t np,
  int decompress,
  FILE * inout);
 
void compress_pck(
  char *varname,
  float *array,
  size_t nxy,
  size_t nz,
  int decompress,
  FILE * inout);
 
void compress_zfp(
  char *varname,
  float *array,
  int nx,
  int ny,
  int nz,
  int precision,
  double tolerance,
  int decompress,
  FILE * inout);
 
void compress_zstd(
  char *varname,
  float *array,
  size_t n,
  int decompress,
  FILE * inout);
 
void day2doy(
  const int year,
  const int mon,
  const int day,
  int *doy);
 
void doy2day(
  const int year,
  const int doy,
  int *mon,
  int *day);
 
void fft_help(
  double *fcReal,
  double *fcImag,
  int n);
 
void geo2cart(
  const double z,
  const double lon,
  const double lat,
  double *x);
 
void get_met(
  ctl_t * ctl,
  clim_t * clim,
  double t,
  met_t ** met0,
  met_t ** met1);
 
void get_met_help(
  ctl_t * ctl,
  double t,
  int direct,
  char *metbase,
  double dt_met,
  char *filename);
 
void get_met_replace(
  char *orig,
  char *search,
  char *repl);
 
void get_tropo(
  int met_tropo,
  ctl_t * ctl,
  clim_t * clim,
  met_t * met,
  double *lons,
  int nx,
  double *lats,
  int ny,
  double *pt,
  double *zt,
  double *tt,
  double *qt,
  double *o3t,
  double *ps,
  double *zs);
 
void intpol_met_4d_coord(
  met_t * met0,
  float height0[EX][EY][EP],
  float array0[EX][EY][EP],
  met_t * met1,
  float height1[EX][EY][EP],
  float array1[EX][EY][EP],
  double ts,
  double height,
  double lon,
  double lat,
  double *var,
  int *ci,
  double *cw,
  int init);
 
void intpol_met_space_3d(
  met_t * met,
  float array[EX][EY][EP],
  double p,
  double lon,
  double lat,
  double *var,
  int *ci,
  double *cw,
  int init);
 
void intpol_met_space_3d_ml(
  met_t * met,
  float array[EX][EY][EP],
  double p,
  double lon,
  double lat,
  double *var);
 
void intpol_met_space_2d(
  met_t * met,
  float array[EX][EY],
  double lon,
  double lat,
  double *var,
  int *ci,
  double *cw,
  int init);
 
void intpol_met_time_3d(
  met_t * met0,
  float array0[EX][EY][EP],
  met_t * met1,
  float array1[EX][EY][EP],
  double ts,
  double p,
  double lon,
  double lat,
  double *var,
  int *ci,
  double *cw,
  int init);
 
void intpol_met_time_3d_ml(
  met_t * met0,
  float array0[EX][EY][EP],
  met_t * met1,
  float array1[EX][EY][EP],
  double ts,
  double p,
  double lon,
  double lat,
  double *var);
 
void intpol_met_time_2d(
  met_t * met0,
  float array0[EX][EY],
  met_t * met1,
  float array1[EX][EY],
  double ts,
  double lon,
  double lat,
  double *var,
  int *ci,
  double *cw,
  int init);
 
void intpol_tropo_3d(
  double time0,
  float array0[EX][EY],
  double time1,
  float array1[EX][EY],
  double lons[EX],
  double lats[EY],
  int nlon,
  int nlat,
  double time,
  double lon,
  double lat,
  int method,
  double *var,
  double *sigma);
 
void jsec2time(
  const double jsec,
  int *year,
  int *mon,
  int *day,
  int *hour,
  int *min,
  int *sec,
  double *remain);
 
double kernel_weight(
  const double kz[EP],
  const double kw[EP],
  const int nk,
  const double p);
 
double lapse_rate(
  const double t,
  const double h2o);
 
void level_definitions(
  ctl_t * ctl);
 
int locate_irr(
  const double *xx,
  const int n,
  const double x);
 
int locate_irr_float(
  const float *xx,
  const int n,
  const double x,
  const int ig);
 
int locate_irr_3d(
  float profiles[EX][EY][EP],
  int np,
  int ind_lon,
  int ind_lat,
  double x);
 
int locate_reg(
  const double *xx,
  const int n,
  const double x);
 
void locate_vert(
  float profiles[EX][EY][EP],
  int np,
  int lon_ap_ind,
  int lat_ap_ind,
  double alt_ap,
  int *ind);
 
void module_advect(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_advect_init(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_bound_cond(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_chemgrid(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double t);
 
void module_chem_init(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_convection(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt,
  double *rs);
 
void module_decay(
  ctl_t * ctl,
  clim_t * clim,
  atm_t * atm,
  double *dt);
 
void module_diffusion_meso(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  cache_t * cache,
  double *dt,
  double *rs);
 
void module_diffusion_turb(
  ctl_t * ctl,
  clim_t * clim,
  atm_t * atm,
  double *dt,
  double *rs);
 
void module_dry_deposition(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_h2o2_chem(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_isosurf_init(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  cache_t * cache);
 
void module_isosurf(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  cache_t * cache,
  double *dt);
 
void module_kpp_chem(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_meteo(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_mixing(
  ctl_t * ctl,
  clim_t * clim,
  atm_t * atm,
  double t);
 
void module_mixing_help(
  ctl_t * ctl,
  clim_t * clim,
  atm_t * atm,
  const int *ixs,
  const int *iys,
  const int *izs,
  int qnt_idx);
 
void module_oh_chem(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_position(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_rng_init(
  int ntask);
 
void module_rng(
  ctl_t * ctl,
  double *rs,
  size_t n,
  int method);
 
void module_sedi(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_sort(
  ctl_t * ctl,
  met_t * met0,
  atm_t * atm);
 
void module_sort_help(
  double *a,
  const int *p,
  const int np);
 
void module_timesteps(
  ctl_t * ctl,
  met_t * met0,
  atm_t * atm,
  double *dt,
  double t);
 
void module_timesteps_init(
  ctl_t * ctl,
  atm_t * atm);
 
void module_tracer_chem(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
void module_wet_deposition(
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double *dt);
 
double nat_temperature(
  const double p,
  const double h2o,
  const double hno3);
 
int read_atm(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
int read_atm_asc(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
int read_atm_bin(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
int read_atm_clams(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
int read_atm_nc(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
void read_clim(
  ctl_t * ctl,
  clim_t * clim);
 
void read_clim_photo(
  const char *filename,
  clim_photo_t * photo);
 
int read_clim_ts(
  const char *filename,
  clim_ts_t * ts);
 
void read_clim_zm(
  const char *filename,
  char *varname,
  clim_zm_t * zm);
 
void read_ctl(
  const char *filename,
  int argc,
  char *argv[],
  ctl_t * ctl);
 
void read_kernel(
  const char *filename,
  double kz[EP],
  double kw[EP],
  int *nk);
 
int read_met(
  const char *filename,
  ctl_t * ctl,
  clim_t * clim,
  met_t * met);
 
void read_met_bin_2d(
  FILE * in,
  met_t * met,
  float var[EX][EY],
  char *varname);
 
void read_met_bin_3d(
  FILE * in,
  ctl_t * ctl,
  met_t * met,
  float var[EX][EY][EP],
  char *varname,
  float bound_min,
  float bound_max);
 
void read_met_cape(
  clim_t * clim,
  met_t * met);
 
void read_met_cloud(
  met_t * met);
 
void read_met_detrend(
  ctl_t * ctl,
  met_t * met);
 
void read_met_extrapolate(
  met_t * met);
 
void read_met_geopot(
  ctl_t * ctl,
  met_t * met);
 
void read_met_grid(
  const char *filename,
  int ncid,
  ctl_t * ctl,
  met_t * met);
 
void read_met_levels(
  int ncid,
  ctl_t * ctl,
  met_t * met);
 
void read_met_ml2pl(
  ctl_t * ctl,
  met_t * met,
  float var[EX][EY][EP],
  char *varname);
 
void read_met_monotonize(
  met_t * met);
 
int read_met_nc_2d(
  int ncid,
  char *varname,
  char *varname2,
  char *varname3,
  char *varname4,
  ctl_t * ctl,
  met_t * met,
  float dest[EX][EY],
  float scl,
  int init);
 
int read_met_nc_3d(
  int ncid,
  char *varname,
  char *varname2,
  char *varname3,
  char *varname4,
  ctl_t * ctl,
  met_t * met,
  float dest[EX][EY][EP],
  float scl);
 
void read_met_pbl(
  met_t * met);
 
void read_met_periodic(
  met_t * met);
 
void read_met_polar_winds(
  met_t * met);
 
void read_met_pv(
  met_t * met);
 
void read_met_ozone(
  met_t * met);
 
void read_met_sample(
  ctl_t * ctl,
  met_t * met);
 
void read_met_surface(
  int ncid,
  met_t * met,
  ctl_t * ctl);
 
void read_met_tropo(
  ctl_t * ctl,
  clim_t * clim,
  met_t * met);
 
void read_obs(
  const char *filename,
  ctl_t * ctl,
  double *rt,
  double *rz,
  double *rlon,
  double *rlat,
  double *robs,
  int *nobs);
 
void read_obs_asc(
  const char *filename,
  double *rt,
  double *rz,
  double *rlon,
  double *rlat,
  double *robs,
  int *nobs);
 
void read_obs_nc(
  const char *filename,
  double *rt,
  double *rz,
  double *rlon,
  double *rlat,
  double *robs,
  int *nobs);
 
double scan_ctl(
  const char *filename,
  int argc,
  char *argv[],
  const char *varname,
  int arridx,
  const char *defvalue,
  char *value);
 
double sedi(
  const double p,
  const double T,
  const double rp,
  const double rhop);
 
void spline(
  const double *x,
  const double *y,
  const int n,
  const double *x2,
  double *y2,
  const int n2,
  const int method);
 
float stddev(
  const float *data,
  const int n);
 
double sza_calc(
  const double sec,
  const double lon,
  const double lat);
 
void time2jsec(
  const int year,
  const int mon,
  const int day,
  const int hour,
  const int min,
  const int sec,
  const double remain,
  double *jsec);
 
void timer(
  const char *name,
  const char *group,
  int output);
 
double time_from_filename(
  const char *filename,
  int offset);
 
double tropo_weight(
  const clim_t * clim,
  const double t,
  const double lat,
  const double p);
 
void write_atm(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_atm_asc(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_atm_bin(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
void write_atm_clams(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
void write_atm_clams_traj(
  const char *dirname,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_atm_nc(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm);
 
void write_csi(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_ens(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_grid(
  const char *filename,
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double t);
 
void write_grid_asc(
  const char *filename,
  ctl_t * ctl,
  double *cd,
  double *mean[NQ],
  double *sigma[NQ],
  double *vmr_impl,
  double t,
  double *z,
  double *lon,
  double *lat,
  double *area,
  double dz,
  int *np);
 
void write_grid_nc(
  const char *filename,
  ctl_t * ctl,
  double *cd,
  double *mean[NQ],
  double *sigma[NQ],
  double *vmr_impl,
  double t,
  double *z,
  double *lon,
  double *lat,
  double *area,
  double dz,
  int *np);
 
int write_met(
  const char *filename,
  ctl_t * ctl,
  met_t * met);
 
void write_met_bin_2d(
  FILE * out,
  met_t * met,
  float var[EX][EY],
  char *varname);
 
void write_met_bin_3d(
  FILE * out,
  ctl_t * ctl,
  met_t * met,
  float var[EX][EY][EP],
  char *varname,
  int precision,
  double tolerance);
 
void write_output(
  const char *dirname,
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double t);
 
void write_prof(
  const char *filename,
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double t);
 
void write_sample(
  const char *filename,
  ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  double t);
 
void write_station(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
void write_vtk(
  const char *filename,
  ctl_t * ctl,
  atm_t * atm,
  double t);
 
/* ------------------------------------------------------------
   OpenACC routines...
   ------------------------------------------------------------ */
 
#ifdef _OPENACC
#pragma acc routine (clim_oh)
#pragma acc routine (clim_photo)
#pragma acc routine (clim_tropo)
#pragma acc routine (clim_ts)
#pragma acc routine (clim_zm)
#pragma acc routine (intpol_met_4d_coord)
#pragma acc routine (intpol_met_space_3d)
#pragma acc routine (intpol_met_space_3d_ml)
#pragma acc routine (intpol_met_space_2d)
#pragma acc routine (intpol_met_time_3d)
#pragma acc routine (intpol_met_time_3d_ml)
#pragma acc routine (intpol_met_time_2d)
#pragma acc routine (kernel_weight)
#pragma acc routine (lapse_rate)
#pragma acc routine (locate_irr)
#pragma acc routine (locate_irr_float)
#pragma acc routine (locate_reg)
#pragma acc routine (locate_vert)
#pragma acc routine (locate_irr_3d)
#pragma acc routine (nat_temperature)
#pragma acc routine (sedi)
#pragma acc routine (stddev)
#pragma acc routine (sza_calc)
#pragma acc routine (tropo_weight)
#endif
 
#endif /* LIBTRAC_H */