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-2026 Forschungszentrum Juelich GmbH
*/
 
#ifndef LIBTRAC_H
#define LIBTRAC_H
 
/* ------------------------------------------------------------
   Includes...
   ------------------------------------------------------------ */
 
#include <ctype.h>
#include <float.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort.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 DD
#include <netcdf_par.h>
#endif
 
#ifdef _OPENACC
#include "openacc.h"
#endif
 
#ifdef CURAND
#include "curand.h"
#endif
 
#ifdef THRUST
#include "thrustsort.h"
#endif
 
#ifdef ZFP
#include "zfp.h"
#endif
 
#ifdef ZSTD
#include "zstd.h"
#endif
 
#ifdef LZ4
#include "lz4.h"
#endif
 
#ifdef SZ3
#include "SZ3c/sz3c.h"
#endif
 
#ifdef CMS
#include "cmultiscale.h"
#endif
 
#ifdef KPP
#include "chem_Parameters.h"
#include "chem_Global.h"
#include "chem_Sparse.h"
#endif
 
#ifdef ECCODES
#include "eccodes.h"
#else
#define codes_handle void*
#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 KARMAN
#define KARMAN 0.40
#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 KAPPA
#define KAPPA 0.286
#endif
 
#ifndef OMEGA_EARTH
#define OMEGA_EARTH 7.2921e-5
#endif
 
#ifndef DOBSON_UNIT
#define DOBSON_UNIT 2.1415e-5
#endif
 
#ifndef M_AIR_MOLECULE
#define M_AIR_MOLECULE 4.8096e-26
#endif
 
#ifndef CHEM_REF_TEMP
#define CHEM_REF_TEMP 298.15
#endif
 
#ifndef H2O2_SO2_RATE_REF
#define H2O2_SO2_RATE_REF 9.1e7
#endif
 
#ifndef H2O2_SO2_RATE_TEMP
#define H2O2_SO2_RATE_TEMP 29700.
#endif
 
#ifndef SO2_HENRY_REF
#define SO2_HENRY_REF 1.3e-2
#endif
 
#ifndef SO2_HENRY_TEMP
#define SO2_HENRY_TEMP 2900.
#endif
 
#ifndef H2O2_HENRY_REF
#define H2O2_HENRY_REF 8.3e2
#endif
 
#ifndef H2O2_HENRY_TEMP
#define H2O2_HENRY_TEMP 7600.
#endif
 
#ifndef SO2_CORR_A
#define SO2_CORR_A 3.12541941e-06
#endif
 
#ifndef SO2_CORR_B
#define SO2_CORR_B -5.72532259e-01
#endif
 
#ifndef O1D_RATE_CCL4_A
#define O1D_RATE_CCL4_A 3.30e-10
#endif
 
#ifndef O1D_RATE_CCL4_B
#define O1D_RATE_CCL4_B 0.
#endif
 
#ifndef O1D_RATE_CFC11_A
#define O1D_RATE_CFC11_A 2.30e-10
#endif
 
#ifndef O1D_RATE_CFC11_B
#define O1D_RATE_CFC11_B 0.
#endif
 
#ifndef O1D_RATE_CFC12_A
#define O1D_RATE_CFC12_A 1.40e-10
#endif
 
#ifndef O1D_RATE_CFC12_B
#define O1D_RATE_CFC12_B -25.
#endif
 
#ifndef O1D_RATE_N2O_A
#define O1D_RATE_N2O_A 1.19e-10
#endif
 
#ifndef O1D_RATE_N2O_B
#define O1D_RATE_N2O_B -20.
#endif
 
#ifndef SO2_DISS_K1_REF
#define SO2_DISS_K1_REF 1.23e-2
#endif
 
#ifndef SO2_DISS_K1_TEMP
#define SO2_DISS_K1_TEMP 2.01e3
#endif
 
#ifndef SO2_DISS_K2_REF
#define SO2_DISS_K2_REF 6e-8
#endif
 
#ifndef SO2_DISS_K2_TEMP
#define SO2_DISS_K2_TEMP 1.12e3
#endif
 
#ifndef WET_DEPO_T_LIQUID
#define WET_DEPO_T_LIQUID T0
#endif
 
#ifndef WET_DEPO_T_ICE
#define WET_DEPO_T_ICE 238.15
#endif
 
#ifndef WET_DEPO_T_LIQUID_BC
#define WET_DEPO_T_LIQUID_BC 270.
#endif
 
#ifndef RADIO_HALF_LIFE_RN222
#define RADIO_HALF_LIFE_RN222 (3.8235 * 86400.0)
#endif
 
#ifndef RADIO_HALF_LIFE_PB210
#define RADIO_HALF_LIFE_PB210 (22.3 * 365.25 * 86400.0)
#endif
 
#ifndef RADIO_HALF_LIFE_BE7
#define RADIO_HALF_LIFE_BE7 (53.22 * 86400.0)
#endif
 
#ifndef RADIO_HALF_LIFE_CS137
#define RADIO_HALF_LIFE_CS137 (30.05 * 365.25 * 86400.0)
#endif
 
#ifndef RADIO_HALF_LIFE_I131
#define RADIO_HALF_LIFE_I131 (8.02 * 86400.0)
#endif
 
#ifndef RADIO_HALF_LIFE_XE133
#define RADIO_HALF_LIFE_XE133 (5.2474 * 86400.0)
#endif
 
#ifndef RADIO_DRY_VDEP_PB210
#define RADIO_DRY_VDEP_PB210 0.002
#endif
 
#ifndef RADIO_DRY_VDEP_BE7
#define RADIO_DRY_VDEP_BE7 0.001
#endif
 
#ifndef RADIO_DRY_VDEP_CS137
#define RADIO_DRY_VDEP_CS137 0.002
#endif
 
#ifndef RADIO_DRY_VDEP_I131
#define RADIO_DRY_VDEP_I131 0.0005
#endif
 
#ifndef RADIO_WET_COEFF_PB210
#define RADIO_WET_COEFF_PB210 1e-4
#endif
 
#ifndef RADIO_WET_COEFF_BE7
#define RADIO_WET_COEFF_BE7 5e-5
#endif
 
#ifndef RADIO_WET_COEFF_CS137
#define RADIO_WET_COEFF_CS137 1e-4
#endif
 
#ifndef RADIO_WET_COEFF_I131
#define RADIO_WET_COEFF_I131 5e-5
#endif
 
#ifndef T0
#define T0 273.15
#endif
 
/* ------------------------------------------------------------
   Dimensions...
   ------------------------------------------------------------ */
 
#ifndef EP
#define EP 140
#endif
 
#ifndef EX
#define EX 1444
#endif
 
#ifndef EY
#define EY 724
#endif
 
#ifndef LEN
#define LEN 5000
#endif
 
#ifndef METVAR
#define METVAR 13
#endif
 
#ifndef NP
#define NP 10000000
#endif
 
#ifndef NQ
#define NQ 15
#endif
 
#ifndef NCSI
#define NCSI 1000000
#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
 
#ifndef DD_EX_GLOB
#define DD_EX_GLOB (EX * 16)
#endif
 
#ifndef DD_EY_GLOB
#define DD_EY_GLOB (EY * 16)
#endif
 
/* ------------------------------------------------------------
   Macros...
   ------------------------------------------------------------ */
 
#ifdef _OPENACC
#define ALLOC(ptr, type, n) do {                                       \
    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!");                                        \
  } while (0)
#else
#define ALLOC(ptr, type, n) do {                                       \
    if (((ptr) = calloc((size_t)(n), sizeof(type))) == NULL)            \
      ERRMSG("Out of memory!");                                        \
  } while (0)
#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 CLAMP(v, lo, hi)                                \
  (((v) < (lo)) ? (lo) : (((v) > (hi)) ? (hi) : (v)))
 
#define COMPRESS_SPEED(nbytes, dt) \
  ((dt) > 0 ? ((double)(nbytes)) / ((dt) * 1024. * 1024.) : NAN)
 
#define COMPRESS_RATIO(raw_size, stored_size) \
  ((double) (raw_size) / (double) (stored_size))
 
#define COMPRESS_BPV(n, stored_size) \
  ((8.0 * (double) (stored_size)) / (double) (n))
 
#define DEG2DX(dlon, lat)                       \
  (RE * DEG2RAD(dlon) * cos(DEG2RAD(lat)))
 
#define DEG2DY(dlat)                            \
  (RE * DEG2RAD(dlat))
 
#define DEG2RAD(deg)                            \
  ((deg) * (M_PI / 180.0))
 
#define DP2DZ(dp, p)                            \
  (- (dp) * H0 / (p))
 
#define DX2DEG(dx, lat)                                         \
  (((lat) < -89.999 || (lat) > 89.999) ? 0                      \
   : (dx) * 180. / (M_PI * RE * cos(DEG2RAD(lat))))
 
#define DY2DEG(dy)                              \
  ((dy) * 180. / (M_PI * RE))
 
#define DZ2DP(dz, p)                            \
  (-(dz) * (p) / H0)
 
#define DX2COORD(met, dx, lat) (((met)->coord_type == 0)? DX2DEG((dx) / 1000.0, (lat)) : ((dx) ))
 
#define DY2COORD(met, dy)  (((met)->coord_type == 0)? DY2DEG((dy) / 1000.0) : ((dy) ))
 
#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 ECC(cmd) {                                                      \
    int ecc_result=(cmd);                                               \
    if(ecc_result!=0)                                                   \
      ERRMSG("ECCODES error: %s", codes_get_error_message(ecc_result)); \
  }
 
#define ECC_READ_2D(variable, target, scaling_factor, found_flag) {     \
    if(strcmp(short_name, variable) == 0) {                             \
      if ((found_flag) == 0)                                            \
        LOG(2, "Read 2-D variable: %s (GRIB)", variable);               \
      for (int ix = 0; ix < met->nx; ix++)                              \
        for (int iy = 0; iy < met->ny; iy++)                            \
          target[ix][iy] = (float)(values[iy * met->nx + ix] * scaling_factor); \
      found_flag = 1;                                                   \
    }                                                                   \
  }
 
#define ECC_READ_3D(variable, level, target, scaling_factor, found_flag) { \
    if(strcmp(short_name, variable) == 0) {                             \
      if ((found_flag) == 0)                                            \
        LOG(2, "Read 3-D variable: %s (GRIB)", variable);               \
      for (int ix = 0; ix < met->nx; ix++)                              \
        for (int iy = 0; iy < met->ny; iy++)                            \
          target[ix][iy][level] = (float) (values[iy * met->nx + ix] * scaling_factor); \
      found_flag += 1;                                                  \
    }                                                                   \
  }
 
#define FMOD(x, y)                              \
  ((x) - (int) ((x) / (y)) * (y))
 
/* Arguments must be side-effect-free; each may be evaluated more than once. */
#define DT_ALIGNED(dt, dt_mod) \
  (isfinite(dt) && isfinite(dt_mod) && (dt_mod) > 0 && \
   (dt) >= (dt_mod) && (dt) / (dt_mod) <= 0x1p52 && \
   fabs((dt) / (dt_mod) - nearbyint((dt) / (dt_mod))) <= \
   fmin(8 * DBL_EPSILON * ((dt) / (dt_mod)), 1e-8))
 
#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->ess, ess, lat, &ess, ci, cw, 0);      \
    intpol_met_space_2d(met, met->nss, nss, lat, &nss, ci, cw, 0);      \
    intpol_met_space_2d(met, met->shf, shf, lat, &shf, 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->ess, met1, met1->ess, time, lon, lat, &ess, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->nss, met1, met1->nss, time, lon, lat, &nss, ci, cw, 0); \
    intpol_met_time_2d(met0, met0->shf, met1, met1->shf, time, lon, lat, &shf, 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 = eastward turbulent surface stress [N/m^2]\n"         \
          "# $20 = northward turbulent surface stress [N/m^2]\n");      \
  fprintf(out,                                                          \
          "# $21 = surface sensible heat flux [W/m^2]\n"                \
          "# $22 = land-sea mask [1]\n"                                 \
          "# $23 = sea surface temperature [K]\n"                       \
          "# $24 = tropopause pressure [hPa]\n"                         \
          "# $25 = tropopause geopotential height [km]\n"               \
          "# $26 = tropopause temperature [K]\n"                        \
          "# $27 = tropopause water vapor [ppv]\n"                      \
          "# $28 = cloud liquid water content [kg/kg]\n"                \
          "# $29 = cloud rain water content [kg/kg]\n"                  \
          "# $30 = cloud ice water content [kg/kg]\n");                 \
  fprintf(out,                                                          \
          "# $31 = cloud snow water content [kg/kg]\n"                  \
          "# $32 = cloud cover [1]\n"                                   \
          "# $33 = total column cloud water [kg/m^2]\n"                 \
          "# $34 = cloud top pressure [hPa]\n"                          \
          "# $35 = cloud bottom pressure [hPa]\n"                       \
          "# $36 = pressure at lifted condensation level (LCL) [hPa]\n" \
          "# $37 = pressure at level of free convection (LFC) [hPa]\n"  \
          "# $38 = pressure at equilibrium level (EL) [hPa]\n"          \
          "# $39 = convective available potential energy (CAPE) [J/kg]\n" \
          "# $40 = convective inhibition (CIN) [J/kg]\n");              \
  fprintf(out,                                                          \
          "# $41 = relative humidity over water [%%]\n"                 \
          "# $42 = relative humidity over ice [%%]\n"                   \
          "# $43 = dew point temperature [K]\n"                         \
          "# $44 = frost point temperature [K]\n"                       \
          "# $45 = NAT temperature [K]\n"                               \
          "# $46 = HNO3 volume mixing ratio [ppv]\n"                    \
          "# $47 = OH volume mixing ratio [ppv]\n"                      \
          "# $48 = H2O2 volume mixing ratio [ppv]\n"                    \
          "# $49 = HO2 volume mixing ratio [ppv]\n"                     \
          "# $50 = O(1D) volume mixing ratio [ppv]\n");                 \
  fprintf(out,                                                          \
          "# $51 = boundary layer pressure [hPa]\n"                     \
          "# $52 = total column ozone [DU]\n"                           \
          "# $53 = number of data points\n"                             \
          "# $54 = number of tropopause data points\n"                  \
          "# $55 = 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, level, quant) { \
    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)); \
    NC_DEF_VAR_QUANTIZE(quant); \
    if((level) != 0) {                                                  \
      NC(nc_def_var_deflate(ncid, varid, 1, 1, level));                 \
      /* unsigned int ulevel = (unsigned int)level; */                  \
      /* NC(nc_def_var_filter(ncid, varid, 32015, 1, (unsigned int[]){ulevel})); */ \
    }                                                                   \
  }
 
#ifdef NC_QUANTIZE_GRANULARBR
#define NC_DEF_VAR_QUANTIZE(quant) do { \
    if((quant) > 0) \
      NC(nc_def_var_quantize(ncid, varid, NC_QUANTIZE_GRANULARBR, quant)); \
  } while(0)
#else
#define NC_DEF_VAR_QUANTIZE(quant) do { \
    if((quant) > 0) \
      ERRMSG("NetCDF quantization is not supported by this NetCDF version"); \
  } while(0)
#endif
 
#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, check) {       \
    int dimid; size_t naux;                               \
    NC(nc_inq_dimid(ncid, dimname, &dimid));              \
    NC(nc_inq_dimlen(ncid, dimid, &naux));                \
    *ptr = (int)naux;                                     \
    if (check)                                            \
      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))
 
#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 || cache->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 || cache->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 RAD2DEG(rad)                            \
  ((rad) * (180.0 / M_PI))
 
#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 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), KAPPA))
 
#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 USAGE                                                           \
  do {                                                                  \
    int iusage;                                                         \
    for (iusage = 1; iusage < argc; iusage++)                           \
      if (!strcmp(argv[iusage], "-h")                                   \
          || !strcmp(argv[iusage], "--help")) {                         \
        usage();                                                        \
        return EXIT_SUCCESS;                                            \
      }                                                                 \
  } while (0)
 
#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)                 \
  timer(id, group, 0);
 
/* ------------------------------------------------------------
   Structs...
   ------------------------------------------------------------ */
 
typedef struct {
 
  /* ------------------------------------------------------------
     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_ess;
 
  int qnt_nss;
 
  int qnt_shf;
 
  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_iso_var;
 
  int qnt_zeta;
 
  int qnt_zeta_d;
 
  int qnt_zeta_dot;
 
  int qnt_eta;
 
  int qnt_eta_d;
 
  int qnt_eta_dot;
 
  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 qnt_Arn222;
 
  int qnt_Apb210;
 
  int qnt_Abe7;
 
  int qnt_Acs137;
 
  int qnt_Ai131;
 
  int qnt_Axe133;
 
  int qnt_current_subdomain;
 
  int qnt_target_subdomain;
 
  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_coord_type;
 
  double met_utm_ref_lat;
 
  double met_utm_ref_lon;
 
  int met_vert_coord;
 
  int met_gp2z;
 
  int met_type;
 
  int met_clams;
 
  int met_nc_scale;
 
  int met_nc_read_mode;
 
  int met_nc_level;
 
  int met_nc_quant;
 
  int met_zstd_level;
 
  int met_zstd_nworkers;
 
  int met_lz4_accel;
 
  int met_pck_zstd;
 
  int met_zfp_prec[METVAR];
 
  double met_zfp_tol[METVAR];
 
  int met_sz3_prec[METVAR];
 
  double met_sz3_tol[METVAR];
 
  int met_lossy_scale[METVAR];
 
  char met_comp_logfile[LEN];
 
  int met_cms_batch;
 
  int met_cms_zstd;
 
  int met_cms_nd0x;
 
  int met_cms_nd0y;
 
  int met_cms_maxlev;
 
  double met_cms_eps[METVAR];
 
  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_press_level_def;
 
  int met_nlev;
 
  double met_lev_hyam[EP];
 
  double met_lev_hybm[EP];
 
  int met_geopot_sx;
 
  int met_geopot_sy;
 
  int met_relhum;
 
  int met_cape;
 
  int met_pbl;
 
  double met_pbl_min;
 
  double met_pbl_max;
 
  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 advect_vert_coord;
 
  int rng_type;
 
  int diffusion;
 
  int turb_pbl_scheme;
 
  double turb_dx_pbl;
 
  double turb_dx_trop;
 
  double turb_dx_strat;
 
  double turb_dz_pbl;
 
  double turb_dz_trop;
 
  double turb_dz_strat;
 
  double turb_mesox;
 
  double turb_mesoz;
 
  double turb_pbl_trans;
 
  int conv_mix_pbl;
 
  double conv_pbl_trans;
 
  double conv_cape;
 
  double conv_cin;
 
  double conv_dt;
 
  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_tau_trop;
 
  double mixing_tau_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;
 
  int radio_decay;
 
  int radio_depo;
 
  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[2];
 
  double wet_depo_bc_h[2];
 
  double wet_depo_so2_ph;
 
  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 depo_basename[LEN];
 
  double depo_dt_out;
 
  int depo_type;
 
  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 atm_nc_level;
 
  int atm_nc_quant[NQ];
 
  int obs_type;
 
  char budget_basename[LEN];
 
  double budget_dt_out;
 
  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;
 
  int nens;
 
  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_nc_level;
 
  int grid_nc_quant[NQ];
 
  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;
 
  /* ------------------------------------------------------------
     Domain decomposition...
     ------------------------------------------------------------ */
 
  int dd;
 
  int dd_subdomains_zonal;
 
  int dd_subdomains_meridional;
 
  int dd_halos_size;
 
  double dd_sort_dt;
 
  int dd_exchange_mode;
 
  int dd_supported;
 
} 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 time;
 
  double p;
 
  double lon;
 
  double lat;
 
  double dt;
 
  double q[NQ];
 
  float uvwp[3];
 
} particle_t;
 
 
typedef struct {
 
  double iso_ps[NP];
 
  double iso_ts[NP];
 
  int iso_n;
 
  float uvwp[NP][3];
 
  double rs[3 * NP + 1];
 
  double dt[NP];
 
} cache_t;
 
typedef struct {
 
  double Apb210[EX * EY];
 
  double Abe7[EX * EY];
 
  double Acs137[EX * EY];
 
  double Ai131[EX * EY];
 
} depo_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;
 
  /* Grid type: 0=lat/lon [deg], 1=UTM [m]. */
  int coord_type;
 
  int nx;
 
  int ny;
 
  int np;
 
  int npl;
 
  double lon[EX];
 
  double lat[EY];
 
  double p[EP];
 
  double hybrid[EP];
 
  double hyam[EP];
 
  double hybm[EP];
 
  double eta[EP];
 
  float ps[EX][EY];
 
  float ts[EX][EY];
 
  float zs[EX][EY];
 
  float us[EX][EY];
 
  float vs[EX][EY];
 
  float ess[EX][EY];
 
  float nss[EX][EY];
 
  float shf[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;
 
typedef struct {
 
  /* ------------------------------------------------------------
     Global grid...
     ------------------------------------------------------------ */
 
  int nx_glob;
 
  int ny_glob;
 
  double lon_glob[DD_EX_GLOB];
 
  double lat_glob[DD_EY_GLOB];
 
  /* ------------------------------------------------------------
     Subdomains...
     ------------------------------------------------------------ */
 
  size_t subdomain_start[4];
 
  size_t subdomain_count[4];
 
  size_t halo_bnd_start[4];
 
  size_t halo_bnd_count[4];
 
  int halo_offset_start;
 
  int halo_offset_end;
 
  /* ------------------------------------------------------------
     Helpers...
     ------------------------------------------------------------ */
 
#ifdef DD
 
  MPI_Datatype MPI_Particle;
 
  double sort_key[NP];
 
  int perm[NP];
 
  double tmp[NP];
 
#endif
 
} dd_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(
  const ctl_t * ctl,
  clim_t * clim);
 
double clim_photo(
  const double rate[CP][CSZA][CO3],
  const clim_photo_t * photo,
  const double p,
  const double sza,
  const 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_log_level(
  FILE * out,
  const char *codec,
  const char *varname,
  const size_t lev,
  const double plev,
  const double ratio,
  const double bpv,
  const double t_comp,
  const double t_decomp,
  const size_t n,
  const size_t nbytes,
  const float *org,
  const float *cmp);
 
void compress_log_levels_3d(
  FILE * out,
  const char *codec,
  const char *varname,
  const met_t * met,
  const float *org_all,
  const float *cmp_all,
  const size_t nxy,
  const size_t nz,
  const double ratio,
  const double bpv,
  const double t_comp,
  const double t_decomp,
  const size_t nbytes);
 
void compress_cms(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
void compress_scale_to_unit(
  float *array,
  const size_t nxy,
  const size_t nz,
  double *off,
  double *scl);
 
void compress_unscale_from_unit(
  float *array,
  const size_t nxy,
  const size_t nz,
  const double *off,
  const double *scl);
 
int compress_read_lossy_scale(
  FILE * in,
  const size_t nz,
  double **off,
  double **scl);
 
void compress_write_lossy_scale(
  FILE * out,
  const int enabled,
  float *array,
  const size_t nxy,
  const size_t nz,
  double **off,
  double **scl);
 
#ifdef ZSTD
ZSTD_CCtx *compress_zstd_create_cctx(
  const int level,
  const int nworkers);
#endif
 
int compress_metvar_index(
  const char *varname);
 
void compress_pck(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
void compress_sz3(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
void compress_zfp(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
void compress_zstd(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
void compress_lz4(
  const ctl_t * ctl,
  const met_t * met,
  const char *varname,
  float *array,
  const int decompress,
  FILE * level_log,
  FILE * inout);
 
double cos_sza(
  const double sec,
  const double lon,
  const double lat);
 
void day2doy(
  const int year,
  const int mon,
  const int day,
  int *doy);
 
void dd_assign_subdomains(
  const ctl_t * ctl,
  const dd_t * dd,
  atm_t * atm,
  const int init);
 
void dd_atm2particles(
  const ctl_t * ctl,
  cache_t * cache,
  atm_t * atm,
  particle_t * particles,
  const int npart);
 
void dd_pack_compact(
  const ctl_t * ctl,
  cache_t * cache,
  atm_t * atm,
  dd_t * dd,
  particle_t * particles,
  const int nkeep,
  const int npart);
 
int dd_calc_subdomain_from_coords(
  const ctl_t * ctl,
  const dd_t * dd,
  const double lon,
  const double lat);
 
void dd_communicate_particles(
  const ctl_t * ctl,
  const dd_t * dd,
  particle_t ** particles,
  int *npart,
  int *capacity);
 
void dd_init(
  const ctl_t * ctl,
  dd_t * dd,
  atm_t * atm);
 
void dd_normalize_lon_lat(
  const dd_t * dd,
  double *lon,
  double *lat);
 
void dd_particles2atm(
  const ctl_t * ctl,
  cache_t * cache,
  const particle_t * particles,
  const int npart,
  atm_t * atm);
 
 
void dd_push(
  const ctl_t * ctl,
  atm_t * atm,
  cache_t * cache,
  int *npart);
 
void dd_count_outgoing(
  const ctl_t * ctl,
  const atm_t * atm,
  dd_t * dd,
  int *nkeep,
  int *npart);
 
 
void dd_read_met_nc_grid(
  dd_t * dd,
  const ctl_t * ctl,
  met_t * met,
  const int ncid);
 
void dd_sort(
  const ctl_t * ctl,
  const met_t * met0,
  cache_t * cache,
  atm_t * atm,
  dd_t * dd,
  int *npart);
 
void doy2day(
  const int year,
  const int doy,
  int *mon,
  int *day);
 
void fft_help(
  double *fcReal,
  double *fcImag,
  const int n);
 
void geo2cart(
  const double z,
  const double lon,
  const double lat,
  double *x);
 
void get_met_filename(
  const ctl_t * ctl,
  const double t,
  const int direct,
  const char *metbase,
  const double dt_met,
  char *filename);
 
void get_met_replace(
  char *orig,
  const char *search,
  const char *repl);
 
void get_tropo(
  const int met_tropo,
  ctl_t * ctl,
  const clim_t * clim,
  met_t * met,
  const double *lons,
  const int nx,
  const double *lats,
  const int ny,
  double *pt,
  double *zt,
  double *tt,
  double *qt,
  double *o3t,
  double *ps,
  double *zs);
 
void intpol_check_lon_lat(
  const double *lons,
  const int nlon,
  const double *lats,
  const int nlat,
  const double lon,
  const double lat,
  double *lon2,
  double *lat2);
 
void intpol_check_cartesian(
  const double *lons,
  const int nlon,
  const double *lats,
  const int nlat,
  const double lon,
  const double lat,
  double *lon2,
  double *lat2);
 
void intpol_met_4d_zeta(
  const met_t * met0,
  float height0[EX][EY][EP],
  float array0[EX][EY][EP],
  const met_t * met1,
  float height1[EX][EY][EP],
  float array1[EX][EY][EP],
  const double ts,
  const double height,
  const double lon,
  const double lat,
  double *var,
  int *ci,
  double *cw,
  const int init);
 
void intpol_met_space_3d(
  const met_t * met,
  float array[EX][EY][EP],
  const double p,
  const double lon,
  const double lat,
  double *var,
  int *ci,
  double *cw,
  const int init);
 
void intpol_met_space_2d(
  const met_t * met,
  float array[EX][EY],
  const double lon,
  const double lat,
  double *var,
  int *ci,
  double *cw,
  const int init);
 
void intpol_met_time_3d(
  const met_t * met0,
  float array0[EX][EY][EP],
  const met_t * met1,
  float array1[EX][EY][EP],
  const double ts,
  const double p,
  const double lon,
  const double lat,
  double *var,
  int *ci,
  double *cw,
  const int init);
 
void intpol_met_time_2d(
  const met_t * met0,
  float array0[EX][EY],
  const met_t * met1,
  float array1[EX][EY],
  const double ts,
  const double lon,
  const double lat,
  double *var,
  int *ci,
  double *cw,
  const int init);
 
void intpol_tropo_3d(
  const double time0,
  float array0[EX][EY],
  const double time1,
  float array1[EX][EY],
  const double lons[EX],
  const double lats[EY],
  const int nlon,
  const int nlat,
  const double time,
  const double lon,
  const double lat,
  const 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_reg(
  const double *xx,
  const int n,
  const double x);
 
void locate_vert(
  float profiles[EX][EY][EP],
  const int np,
  const int lon_ap_ind,
  const int lat_ap_ind,
  const double alt_ap,
  int *ind);
 
void module_advect(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_advect_init(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_bound_cond(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_chem_grid(
  const ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  const double t);
 
void module_chem_init(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_convection(
  const ctl_t * ctl,
  cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_dd(
  double t,
  const ctl_t * ctl,
  cache_t * cache,
  dd_t * dd,
  atm_t * atm,
  met_t ** met);
 
void module_decay(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  atm_t * atm);
 
void module_diff_meso(
  const ctl_t * ctl,
  cache_t * cache,
  const met_t * met0,
  const met_t * met1,
  atm_t * atm);
 
void module_diff_pbl(
  const ctl_t * ctl,
  cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_diff_turb(
  const ctl_t * ctl,
  cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_dry_depo(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_h2o2_chem(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_isosurf_init(
  const ctl_t * ctl,
  cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_isosurf(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_kpp_chem(
  ctl_t * ctl,
  cache_t * cache,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_meteo(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_mixing(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  atm_t * atm,
  const double t);
 
void module_mixing_help(
  const ctl_t * ctl,
  atm_t * atm,
  const int *ixs,
  const int *iys,
  const int *izs,
  const double *mixparam,
  const int qnt_idx,
  const int use_ensemble);
 
void module_oh_chem(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_position(
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  const int reflect);
 
void module_rng_init(
  const int ntask);
 
void module_rng(
  const ctl_t * ctl,
  double *rs,
  const size_t n,
  const int method);
 
void module_radio_decay(
  const ctl_t * ctl,
  const cache_t * cache,
  atm_t * atm);
 
void module_radio_depo(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  depo_t * depo);
 
void module_sedi(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void sort_array(
  double *array,
  const int *perm,
  double *tmp,
  const int np);
 
void sort_uvwp(
  cache_t * cache,
  const int *perm,
  double *tmp,
  const int np);
 
void module_sort(
  const ctl_t * ctl,
  const met_t * met0,
  cache_t * cache,
  atm_t * atm);
 
void module_timesteps(
  const ctl_t * ctl,
  cache_t * cache,
  const met_t * met0,
  atm_t * atm,
  const double t);
 
void module_timesteps_init(
  ctl_t * ctl,
  const atm_t * atm);
 
void module_tracer_chem(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void module_wet_depo(
  const ctl_t * ctl,
  const cache_t * cache,
  met_t * met0,
  met_t * met1,
  atm_t * atm);
 
void mptrac_alloc(
  ctl_t ** ctl,
  cache_t ** cache,
  clim_t ** clim,
  met_t ** met0,
  met_t ** met1,
  atm_t ** atm,
  depo_t ** depo,
  dd_t ** dd);
 
void mptrac_free(
  ctl_t * ctl,
  cache_t * cache,
  clim_t * clim,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  depo_t * depo,
  dd_t * dd);
 
void mptrac_get_met(
  ctl_t * ctl,
  clim_t * clim,
  const double t,
  met_t ** met0,
  met_t ** met1,
  dd_t * dd);
 
void mptrac_init(
  ctl_t * ctl,
  cache_t * cache,
  clim_t * clim,
  atm_t * atm,
  depo_t * depo,
  const int ntask);
 
int mptrac_read_atm(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
void mptrac_read_clim(
  const ctl_t * ctl,
  clim_t * clim);
 
void mptrac_read_ctl(
  const char *filename,
  int argc,
  char *argv[],
  ctl_t * ctl);
 
int mptrac_read_met(
  const char *filename,
  const ctl_t * ctl,
  const clim_t * clim,
  met_t * met,
  dd_t * dd);
 
void mptrac_run_timestep(
  ctl_t * ctl,
  cache_t * cache,
  clim_t * clim,
  met_t ** met0,
  met_t ** met1,
  atm_t * atm,
  depo_t * depo,
  double t,
  dd_t * dd);
 
void mptrac_write_atm(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void mptrac_write_met(
  const char *filename,
  const ctl_t * ctl,
  met_t * met);
 
void mptrac_write_output(
  const char *dirname,
  const ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  atm_t * atm,
  depo_t * depo,
  const double t);
 
void mptrac_update_device(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t ** met0,
  met_t ** met1,
  const atm_t * atm);
 
void mptrac_update_host(
  const ctl_t * ctl,
  const cache_t * cache,
  const clim_t * clim,
  met_t ** met0,
  met_t ** met1,
  const atm_t * atm);
 
double nat_temperature(
  const double p,
  const double h2o,
  const double hno3);
 
double pbl_weight(
  const ctl_t * ctl,
  const atm_t * atm,
  const int ip,
  const double pbl,
  const double ps);
 
int read_atm_asc(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
int read_atm_bin(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
int read_atm_clams(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
int read_atm_nc(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
void read_clim_photo(
  const char *filename,
  clim_photo_t * photo);
 
void read_clim_photo_help(
  const int ncid,
  const char *varname,
  const clim_photo_t * photo,
  double var[CP][CSZA][CO3]);
 
int read_clim_ts(
  const char *filename,
  clim_ts_t * ts);
 
void read_clim_zm(
  const char *filename,
  const char *varname,
  clim_zm_t * zm);
 
void read_kernel(
  const char *filename,
  double kz[EP],
  double kw[EP],
  int *nk);
 
int read_met_bin(
  const char *filename,
  const ctl_t * ctl,
  met_t * met);
 
void read_met_bin_2d(
  FILE * in,
  const met_t * met,
  float var[EX][EY],
  const char *varname);
 
void read_met_bin_3d(
  FILE * in,
  const ctl_t * ctl,
  const met_t * met,
  float var[EX][EY][EP],
  const char *varname,
  const float bound_min,
  const float bound_max);
 
void read_met_cape(
  const ctl_t * ctl,
  const clim_t * clim,
  met_t * met);
 
void read_met_cloud(
  met_t * met);
 
void read_met_detrend(
  const ctl_t * ctl,
  met_t * met);
 
void read_met_extrapolate(
  met_t * met);
 
void read_met_geopot(
  const ctl_t * ctl,
  met_t * met);
 
int read_met_grib(
  const char *filename,
  const ctl_t * ctl,
  met_t * met);
 
void read_met_grib_grid(
  codes_handle ** handles,
  int count_handles,
  met_t * met);
 
void read_met_grib_levels(
  codes_handle ** handles,
  const int num_messages,
  const ctl_t * ctl,
  met_t * met);
 
void read_met_grib_surface(
  codes_handle ** handles,
  const int num_messages,
  const ctl_t * ctl,
  met_t * met);
 
void read_met_ml2pl(
  const ctl_t * ctl,
  const met_t * met,
  float var[EX][EY][EP],
  const char *varname);
 
void read_met_monotonize(
  const ctl_t * ctl,
  met_t * met);
 
int read_met_nc(
  const char *filename,
  const ctl_t * ctl,
  met_t * met,
  dd_t * dd);
 
void read_met_nc_grid(
  const char *filename,
  const int ncid,
  const ctl_t * ctl,
  met_t * met,
  dd_t * dd);
 
void read_met_nc_levels(
  const int ncid,
  const ctl_t * ctl,
  met_t * met,
  dd_t * dd);
 
void read_met_nc_surface(
  const int ncid,
  const ctl_t * ctl,
  met_t * met,
  dd_t * dd);
 
int read_met_nc_2d(
  const int ncid,
  const char *const varnames[],
  const size_t nvarnames,
  const ctl_t * ctl,
  const met_t * met,
  dd_t * dd,
  float dest[EX][EY],
  const float scl,
  const int init);
 
int read_met_nc_3d(
  const int ncid,
  const char *const varnames[],
  const size_t nvarnames,
  const ctl_t * ctl,
  const met_t * met,
  dd_t * dd,
  float dest[EX][EY][EP],
  const float scl);
 
void read_met_pbl(
  const ctl_t * ctl,
  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(
  const ctl_t * ctl,
  met_t * met);
 
void read_met_tropo(
  const ctl_t * ctl,
  const clim_t * clim,
  met_t * met);
 
void read_obs(
  const char *filename,
  const 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,
  const 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);
 
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,
  const int output);
 
double time_from_filename(
  const char *filename,
  const int offset,
  const int with_seconds);
 
double tropo_weight(
  const ctl_t * ctl,
  const clim_t * clim,
  const atm_t * atm,
  const int ip);
 
void write_atm_asc(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_atm_bin(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_atm_clams(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_atm_clams_traj(
  const char *dirname,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_atm_nc(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_budget(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_csi(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_csi_ens(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_ens(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const double t);
 
void write_depo(
  const char *filename,
  const ctl_t * ctl,
  const depo_t * depo,
  const double t);
 
void write_depo_asc(
  const char *filename,
  const ctl_t * ctl,
  const double *data,
  const double t,
  const double *lon,
  const double *lat,
  const double *area);
 
void write_depo_nc(
  const char *filename,
  const ctl_t * ctl,
  const double *data,
  const double t,
  const double *lon,
  const double *lat,
  const double *area);
 
void write_grid(
  const char *filename,
  const ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  const atm_t * atm,
  const double t);
 
void write_grid_asc(
  const char *filename,
  const ctl_t * ctl,
  const double *cd,
  double *mean[NQ],
  double *sigma[NQ],
  const double *vmr_impl,
  const double t,
  const double *z,
  const double *lon,
  const double *lat,
  const double *area,
  const double dz,
  const int *np);
 
void write_grid_nc(
  const char *filename,
  const ctl_t * ctl,
  const double *cd,
  double *mean[NQ],
  double *sigma[NQ],
  const double *vmr_impl,
  const double t,
  const double *z,
  const double *lon,
  const double *lat,
  const double *area,
  const double dz,
  const int *np);
 
void write_met_bin(
  const char *filename,
  const ctl_t * ctl,
  met_t * met);
 
void write_met_bin_2d(
  FILE * out,
  met_t * met,
  float var[EX][EY],
  const char *varname);
 
void write_met_bin_3d(
  FILE * out,
  const ctl_t * ctl,
  met_t * met,
  float var[EX][EY][EP],
  const char *varname,
  const int metvar,
  FILE * level_log);
 
void write_met_nc(
  const char *filename,
  const ctl_t * ctl,
  met_t * met);
 
void write_met_nc_2d(
  const int ncid,
  const char *varname,
  met_t * met,
  float var[EX][EY],
  const float scl);
 
void write_met_nc_3d(
  const int ncid,
  const char *varname,
  met_t * met,
  float var[EX][EY][EP],
  const float scl);
 
void write_prof(
  const char *filename,
  const ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  const atm_t * atm,
  const double t);
 
void write_sample(
  const char *filename,
  const ctl_t * ctl,
  met_t * met0,
  met_t * met1,
  const atm_t * atm,
  const double t);
 
void write_station(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm,
  const double t);
 
void write_vtk(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  const 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 (cos_sza)
#pragma acc routine (dd_calc_subdomain_from_coords)
#pragma acc routine (dd_normalize_lon_lat)
#pragma acc routine (intpol_check_lon_lat)
#pragma acc routine (intpol_met_4d_zeta)
#pragma acc routine (intpol_met_space_3d)
#pragma acc routine (intpol_met_space_2d)
#pragma acc routine (intpol_met_time_3d)
#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 (nat_temperature)
#pragma acc routine (pbl_weight)
#pragma acc routine (sedi)
#pragma acc routine (stddev)
#pragma acc routine (tropo_weight)
#endif
 
#endif /* LIBTRAC_H */