jurassic.h

Go to the documentation of this file.

/*
  This file is part of JURASSIC.
  
  JURASSIC 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.
  
  JURASSIC 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 JURASSIC. If not, see <http://www.gnu.org/licenses/>.
  
  Copyright (C) 2003-2026 Forschungszentrum Juelich GmbH
*/
 
#ifndef JURASSIC_H
#define JURASSIC_H
 
/* ------------------------------------------------------------
   Includes...
   ------------------------------------------------------------ */
 
#include <errno.h>
#include <fcntl.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.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 <unistd.h>
 
enum {
  FORMOD_STATUS_OK = 0,
  FORMOD_STATUS_RFM_UNSUPPORTED = 1,
  FORMOD_STATUS_FOV_DATA_MISSING = 2,
  FORMOD_STATUS_OBSERVER_BELOW_SURFACE = 3,
  FORMOD_STATUS_TOO_MANY_LOS_POINTS = 4
};
 
/* ------------------------------------------------------------
   Constants...
   ------------------------------------------------------------ */
 
#ifndef C1
#define C1 1.19104259e-8
#endif
 
#ifndef C2
#define C2 1.43877506
#endif
 
#ifndef EPSMIN
#define EPSMIN 0
#endif
 
#ifndef EPSMAX
#define EPSMAX 1
#endif
 
#ifndef G0
#define G0 9.80665
#endif
 
#ifndef H0
#define H0 7.0
#endif
 
#ifndef KB
#define KB 1.3806504e-23
#endif
 
#ifndef ME
#define ME 5.976e24
#endif
 
#ifndef NA
#define NA 6.02214199e23
#endif
 
#ifndef N2
#define N2 0.78084
#endif
 
#ifndef O2
#define O2 0.20946
#endif
 
#ifndef P0
#define P0 1013.25
#endif
 
#ifndef RE
#define RE 6367.421
#endif
 
#ifndef RI
#define RI 8.3144598
#endif
 
#ifndef T0
#define T0 273.15
#endif
 
#ifndef TMAX
#define TMAX 400.
#endif
 
#ifndef TMIN
#define TMIN 100.
#endif
 
#ifndef TSUN
#define TSUN 5780.
#endif
 
#ifndef UMIN
#define UMIN 0
#endif
 
#ifndef UMAX
#define UMAX 1e30
#endif
 
/* ------------------------------------------------------------
   Kernel finite-difference perturbations...
   ------------------------------------------------------------ */
 
#ifndef KERNEL_P_REL
#define KERNEL_P_REL 0.01
#endif
 
#ifndef KERNEL_P_MIN
#define KERNEL_P_MIN 1e-7
#endif
 
#ifndef KERNEL_T_ABS
#define KERNEL_T_ABS 1.0
#endif
 
#ifndef KERNEL_Q_REL
#define KERNEL_Q_REL 0.01
#endif
 
#ifndef KERNEL_Q_MIN
#define KERNEL_Q_MIN 1e-15
#endif
 
#ifndef KERNEL_K_ABS
#define KERNEL_K_ABS 1e-4
#endif
 
#ifndef KERNEL_CLZ_ABS
#define KERNEL_CLZ_ABS 0.1
#endif
 
#ifndef KERNEL_CLDZ_ABS
#define KERNEL_CLDZ_ABS 0.1
#endif
 
#ifndef KERNEL_CLK_ABS
#define KERNEL_CLK_ABS 1e-4
#endif
 
#ifndef KERNEL_SFT_ABS
#define KERNEL_SFT_ABS 1.0
#endif
 
#ifndef KERNEL_SFEPS_ABS
#define KERNEL_SFEPS_ABS 1e-2
#endif
 
/* ------------------------------------------------------------
   Dimensions...
   ------------------------------------------------------------ */
 
#ifndef LEN
#define LEN 10000
#endif
 
#ifndef M
#define M (NR*ND)
#endif
 
#ifndef N
#define N ((2 + NG + NW) * NP + NCL + NSF + 3)
#endif
 
#ifndef NCL
#define NCL 8
#endif
 
#ifndef ND
#define ND 128
#endif
 
#ifndef NFOV
#define NFOV 5
#endif
 
#ifndef NG
#define NG 8
#endif
 
#ifndef CT
#define CT 12
#endif
 
#ifndef CZ
#define CZ 151
#endif
 
#ifndef CLAT
#define CLAT 36
#endif
 
#ifndef NLOS
#define NLOS 4096
#endif
 
#ifndef NP
#define NP 256
#endif
 
#ifndef NQ
#define NQ (5 + NG + NW + NCL + NSF)
#endif
 
#ifndef NR
#define NR 256
#endif
 
#ifndef NSF
#define NSF 8
#endif
 
#ifndef NSHAPE
#define NSHAPE 20000
#endif
 
#ifndef NTIMER
#define NTIMER 100
#endif
 
#ifndef NW
#define NW 4
#endif
 
#ifndef RFMNPTS
#define RFMNPTS 10000000
#endif
 
#ifndef TBLNP
#define TBLNP 41
#endif
 
#ifndef TBLNS
#define TBLNS 1200
#endif
 
#ifndef TBLNT
#define TBLNT 30
#endif
 
#ifndef TBLNU
#define TBLNU 512
#endif
 
/* ------------------------------------------------------------
   Quantity indices...
   ------------------------------------------------------------ */
 
#define IDXP 0
 
#define IDXT 1
 
#define IDXQ(ig) (2 + (ig))
 
#define IDXK(iw) (2 + (ctl->ng) + (iw))
 
#define IDXCLZ (2 + (ctl->ng) + (ctl->nw))
 
#define IDXCLDZ (3 + (ctl->ng) + (ctl->nw))
 
#define IDXCLK(icl) (4 + (ctl->ng) + (ctl->nw) + (icl))
 
#define IDXSFT (4 + (ctl->ng) + (ctl->nw) + (ctl->ncl))
 
#define IDXSFEPS(isf) (5 + (ctl->ng) + (ctl->nw) + (ctl->ncl) + (isf))
 
/* ------------------------------------------------------------
   Macros...
   ------------------------------------------------------------ */
 
#define ALLOC(ptr, type, n)                              \
  if((ptr=calloc((size_t)(n), sizeof(type)))==NULL)      \
    ERRMSG("Out of memory!");
 
#define BRIGHT(rad, nu)                                 \
  (C2 * (nu) / gsl_log1p(C1 * POW3(nu) / (rad)))
 
#define CLAMP(v, lo, hi)                                \
  (((v) < (lo)) ? (lo) : (((v) > (hi)) ? (hi) : (v)))
 
#define CO2_FIT(TIME) \
  (373.87884731310215e-6 \
   + 1.9021498612395775e-6 * (((TIME) - 94694400.) / 31557600.) \
   + 0.016490044598136553e-6 \
   * POW2((((TIME) - 94694400.) / 31557600.)))
 
#define DEG2RAD(deg) ((deg) * (M_PI / 180.0))
 
#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, out) { \
    if(fread(ptr, sizeof(type), size, out)!=size) \
      ERRMSG("Error while reading!"); \
  }
 
#define FWRITE(ptr, type, size, out) { \
    if(fwrite(ptr, sizeof(type), size, out)!=size) \
      ERRMSG("Error while writing!"); \
  }
 
#define LIN(x0, y0, x1, y1, x) \
  ((y0)+((y1)-(y0))/((x1)-(x0))*((x)-(x0)))
 
#define LOGY(x0, y0, x1, y1, x) \
  (((y1)/(y0)>0) \
   ? ((y0)*exp(log((y1)/(y0))/((x1)-(x0))*((x)-(x0)))) \
   : LIN(x0, y0, x1, y1, x))
 
#define MAX(a,b) (((a)>(b))?(a):(b))
 
#define MIN(a,b) (((a)<(b))?(a):(b))
 
#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 NEDT(t_bg, nesr, nu)                                    \
  (BRIGHT(PLANCK((t_bg), (nu)) + (nesr), (nu)) - (t_bg))
 
#define NESR(t_bg, nedt, nu)                                    \
  (PLANCK((t_bg) + (nedt), (nu)) - PLANCK((t_bg), (nu)))
 
#define NORM(a) sqrt(DOTP(a, a))
 
#define PLANCK(T, nu) \
  (C1 * POW3(nu) / gsl_expm1(C2 * (nu) / (T)))
 
#define POW2(x) ((x)*(x))
 
#define POW3(x) ((x)*(x)*(x))
 
#define RAD2DEG(rad) ((rad) * (180.0 / M_PI))
 
#define REFRAC_SW53(p, T, h2o) \
  (1e-6 * (77.6 * (p) / (T) + 3.73e5 * (h2o) * (p) / ((T) * (T))))
 
#define REFRAC_EDLEN(p, T)                                                      \
  ((1e-8                                                                \
    * (8342.54 + 2406147.0 / (130.0 - 1.0 / (10.0 * 10.0))              \
       + 15998.0 / (38.9 - 1.0 / (10.0 * 10.0))))                       \
   * ((p) / 1013.25) * (288.15 / (T)))
 
#define TBL_LOG(tbl, id, ig)                                            \
  do {                                                                  \
    for (int ip = 0; ip < (tbl)->np[(id)][(ig)]; ip++)                  \
      LOG(2,                                                            \
          "p[%2d]= %.5e hPa | T[0:%2d]= %.2f ... %.2f K | "             \
          "u[0:%3d]= %.5e ... %.5e molec/cm^2 | "                       \
          "eps[0:%3d]= %.5e ... %.5e",                                  \
          (ip),                                                         \
          (tbl)->p[(id)][(ig)][(ip)],                                   \
          (tbl)->nt[(id)][(ig)][(ip)] - 1,                              \
          (tbl)->t[(id)][(ig)][(ip)][0],                                \
          (tbl)->t[(id)][(ig)][(ip)]                                    \
          [(tbl)->nt[(id)][(ig)][(ip)] - 1],                            \
          (tbl)->nu[(id)][(ig)][(ip)][0] - 1,                           \
          exp((tbl)->logu[(id)][(ig)][(ip)][0][0]),                     \
          exp((tbl)->logu[(id)][(ig)][(ip)][0]                          \
              [(tbl)->nu[(id)][(ig)][(ip)][0] - 1]),                    \
          (tbl)->nu[(id)][(ig)][(ip)][0] - 1,                           \
          exp((tbl)->logeps[(id)][(ig)][(ip)][0][0]),                   \
          exp((tbl)->logeps[(id)][(ig)][(ip)][0]                        \
              [(tbl)->nu[(id)][(ig)][(ip)][0] - 1]));                   \
  } while (0)
 
#define SELECT_TIMER(id, group) \
  {timer_group(id, group, 0);}
 
#define PRINT_TIMERS \
  {timer_group("END", "END", 1);}
 
#define TOK(line, tok, format, var) { \
    if(((tok)=strtok((line), " \t"))) { \
      if(sscanf(tok, format, &(var))!=1) continue; \
    } else ERRMSG("Error while reading!"); \
  }
 
#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)
 
/* ------------------------------------------------------------
   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);
 
/* ------------------------------------------------------------
   Structs...
   ------------------------------------------------------------ */
 
typedef struct {
 
  int np;
 
  double time[NP];
 
  double z[NP];
 
  double lon[NP];
 
  double lat[NP];
 
  double p[NP];
 
  double t[NP];
 
  double q[NG][NP];
 
  double k[NW][NP];
 
  double clz;
 
  double cldz;
 
  double clk[NCL];
 
  double sft;
 
  double sfeps[NSF];
 
} atm_t;
 
typedef struct {
  int nalt;
 
  int nlat;
 
  int ntime;
 
  double alt[CZ];
 
  double lat[CLAT];
 
  double time[CT];
 
  double t[CT][CZ][CLAT];
 
  double q[NG][CT][CZ][CLAT];
 
  int has_t;
 
  int has_q[NG];
 
} clim_t;
 
typedef struct {
 
  int ng;
 
  char emitter[NG][LEN];
 
  int ig_co2;
 
  int ig_h2o;
 
  int ig_n2;
 
  int ig_o2;
 
  int nd;
 
  double nu[ND];
 
  int nw;
 
  int window[ND];
 
  int ncl;
 
  double clnu[NCL];
 
  int nsf;
 
  double sfnu[NSF];
 
  int sftype;
 
  double sfsza;
 
  char tblbase[LEN];
 
  int tblfmt;
 
  char climfile[LEN];
 
  int atmfmt;
 
  int obsfmt;
 
  int matrixfmt;
 
  double hydz;
 
  int ctm_co2;
 
  int ctm_h2o;
 
  int ctm_n2;
 
  int ctm_o2;
 
  int refrac;
 
  double rayds;
 
  double raydz;
 
  char fov[LEN];
 
  double fov_dz[NSHAPE];
 
  double fov_w[NSHAPE];
 
  int fov_n;
 
  double retp_zmin;
 
  double retp_zmax;
 
  double rett_zmin;
 
  double rett_zmax;
 
  double retq_zmin[NG];
 
  double retq_zmax[NG];
 
  double retk_zmin[NW];
 
  double retk_zmax[NW];
 
  int ret_clz;
 
  int ret_cldz;
 
  int ret_clk;
 
  int ret_sft;
 
  int ret_sfeps;
 
  int write_bbt;
 
  int write_matrix;
 
  int formod;
 
  char rfmbin[LEN];
 
  char rfmhit[LEN];
 
  char rfmxsc[NG][LEN];
 
} ctl_t;
 
typedef struct {
 
  int np;
 
  double z[NLOS];
 
  double lon[NLOS];
 
  double lat[NLOS];
 
  double p[NLOS];
 
  double t[NLOS];
 
  double q[NLOS][NG];
 
  double k[NLOS][ND];
 
  double sft;
 
  double sfeps[ND];
 
  double ds[NLOS];
 
  double u[NLOS][NG];
 
  double cgp[NLOS][NG];
 
  double cgt[NLOS][NG];
 
  double cgu[NLOS][NG];
 
  double eps[NLOS][ND];
 
  double src[NLOS][ND];
 
} los_t;
 
typedef struct {
 
  int nr;
 
  double time[NR];
 
  double obsz[NR];
 
  double obslon[NR];
 
  double obslat[NR];
 
  double vpz[NR];
 
  double vplon[NR];
 
  double vplat[NR];
 
  double tpz[NR];
 
  double tplon[NR];
 
  double tplat[NR];
 
  double tau[ND][NR];
 
  double rad[ND][NR];
 
  int mask[ND][NR];
 
} obs_t;
 
typedef struct {
 
  char dir[LEN];
 
  int kernel_recomp;
 
  int conv_itmax;
 
  double conv_dmin;
 
  int err_ana;
 
  double err_formod[ND];
 
  double err_noise[ND];
 
  double err_press;
 
  double err_press_cz;
 
  double err_press_ch;
 
  double err_temp;
 
  double err_temp_cz;
 
  double err_temp_ch;
 
  double err_q[NG];
 
  double err_q_cz[NG];
 
  double err_q_ch[NG];
 
  double err_k[NW];
 
  double err_k_cz[NW];
 
  double err_k_ch[NW];
 
  double err_clz;
 
  double err_cldz;
 
  double err_clk[NCL];
 
  double err_sft;
 
  double err_sfeps[NSF];
 
  int shared_io_profile;
 
  char shared_io_proflist[LEN];
 
  char shared_io_atm_apr_file[LEN];
 
  char shared_io_obs_meas_file[LEN];
 
  char shared_io_atm_final_file[LEN];
 
  char shared_io_obs_final_file[LEN];
 
  char shared_io_matrix_cov_apr_file[LEN];
 
  char shared_io_matrix_kernel_file[LEN];
 
  char shared_io_matrix_cov_ret_file[LEN];
 
  char shared_io_matrix_corr_file[LEN];
 
  char shared_io_matrix_gain_file[LEN];
 
  char shared_io_matrix_avk_file[LEN];
 
  char shared_io_atm_err_total_file[LEN];
 
  char shared_io_atm_err_noise_file[LEN];
 
  char shared_io_atm_err_formod_file[LEN];
 
  char shared_io_atm_cont_file[LEN];
 
  char shared_io_atm_res_file[LEN];
 
} ret_t;
 
typedef struct {
 
  int np[ND][NG];
 
  int nt[ND][NG][TBLNP];
 
  int nu[ND][NG][TBLNP][TBLNT];
 
  double p[ND][NG][TBLNP];
 
  double lnp[ND][NG][TBLNP];
 
  double t[ND][NG][TBLNP][TBLNT];
 
  float *logu[ND][NG][TBLNP][TBLNT];
 
  float *logeps[ND][NG][TBLNP][TBLNT];
 
  int filt_n[ND];
 
  double filt_nu[ND][NSHAPE];
 
  double filt_f[ND][NSHAPE];
 
  double st[TBLNS];
 
  double sr[TBLNS][ND];
 
} tbl_t;
 
/* ------------------------------------------------------------
   Functions...
   ------------------------------------------------------------ */
 
void analyze_avk(
  const ret_t * ret,
  const ctl_t * ctl,
  const atm_t * atm,
  const int *iqa,
  const int *ipa,
  const gsl_matrix * avk);
 
void analyze_avk_quantity(
  const gsl_matrix * avk,
  const int iq,
  const int *ipa,
  const size_t *n0,
  const size_t *n1,
  double *cont,
  double *res);
 
size_t atm2x(
  const ctl_t * ctl,
  const atm_t * atm,
  gsl_vector * x,
  int *iqa,
  int *ipa);
 
void atm2x_help(
  const double value,
  const int value_iqa,
  const int value_ip,
  gsl_vector * x,
  int *iqa,
  int *ipa,
  size_t *n);
 
void cart2geo(
  const double *x,
  double *z,
  double *lon,
  double *lat);
 
void climatology(
  const ctl_t * ctl,
  const clim_t * clim,
  atm_t * atm);
 
int read_clim_field(
  const char *filename,
  const char *varname,
  const int ncid,
  const int dim_month,
  const int dim_alt,
  const int dim_lat,
  const int nalt,
  const int nlat,
  double field[CT][CZ][CLAT],
  double *vmin,
  double *vmax);
 
void climatology_external_interp(
  const clim_t * clim,
  const atm_t * atm,
  const double field[CT][CZ][CLAT],
  const int top_policy,
  double *out);
 
clim_t *read_clim(
  const ctl_t * ctl);
 
void climatology_external(
  const ctl_t * ctl,
  const clim_t * clim,
  atm_t * atm);
 
double cos_sza(
  const double sec,
  const double lon,
  const double lat);
 
double cost_function(
  const gsl_vector * dx,
  const gsl_vector * dy,
  const gsl_matrix * s_a_inv,
  const gsl_vector * sig_eps_inv);
 
double ctmco2(
  const double nu,
  const double p,
  const double t,
  const double u);
 
double ctmh2o(
  const double nu,
  const double p,
  const double t,
  const double q,
  const double u);
 
double ctmn2(
  const double nu,
  const double p,
  const double t);
 
double ctmo2(
  const double nu,
  const double p,
  const double t);
 
void copy_atm(
  const ctl_t * ctl,
  atm_t * atm_dest,
  const atm_t * atm_src,
  const int init);
 
void copy_obs(
  const ctl_t * ctl,
  obs_t * obs_dest,
  const obs_t * obs_src,
  const int init);
 
void day2doy(
  int year,
  int mon,
  int day,
  int *doy);
 
void doy2day(
  int year,
  int doy,
  int *mon,
  int *day);
 
int find_emitter(
  const ctl_t * ctl,
  const char *emitter);
 
int formod(
  const ctl_t * ctl,
  const tbl_t * tbl,
  atm_t * atm,
  obs_t * obs,
  los_t * los_scratch,
  obs_t * obs_scratch);
 
void formod_batch(
  const ctl_t * ctl,
  const tbl_t * tbl,
  atm_t * atm,
  obs_t * obs,
  const int nbatch,
  int *status,
  los_t * los_scratch,
  obs_t * obs_scratch);
 
void formod_continua(
  const ctl_t * ctl,
  const los_t * los,
  const int ip,
  double *beta);
 
int formod_fov(
  const ctl_t * ctl,
  obs_t * obs,
  obs_t * obs_scratch);
 
int formod_pencil(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const atm_t * atm,
  obs_t * obs,
  const int ir,
  los_t * los);
 
void formod_rfm(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const atm_t * atm,
  obs_t * obs);
 
void formod_srcfunc(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const double t,
  double *src);
 
void geo2cart(
  const double z,
  const double lon,
  const double lat,
  double *x);
 
void hydrostatic(
  const ctl_t * ctl,
  atm_t * atm);
 
void idx2name(
  const ctl_t * ctl,
  const int idx,
  char *quantity);
 
void init_srcfunc(
  const ctl_t * ctl,
  tbl_t * tbl);
 
void intpol_atm(
  const ctl_t * ctl,
  const atm_t * atm,
  const double z,
  double *p,
  double *t,
  double *q,
  double *k);
 
void intpol_tbl_cga(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const los_t * los,
  const int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
void intpol_tbl_ega(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const los_t * los,
  const int ip,
  double tau_path[ND][NG],
  double tau_seg[ND]);
 
double intpol_tbl_eps(
  const tbl_t * tbl,
  const int ig,
  const int id,
  const int ip,
  const int it,
  const double logu);
 
double intpol_tbl_u(
  const tbl_t * tbl,
  const int ig,
  const int id,
  const int ip,
  const int it,
  const double logeps);
 
void jsec2time(
  const double jsec,
  int *year,
  int *mon,
  int *day,
  int *hour,
  int *min,
  int *sec,
  double *remain);
 
void kernel(
  const ctl_t * ctl,
  const tbl_t * tbl,
  atm_t * atm,
  obs_t * obs,
  gsl_matrix * k);
 
int locate_irr(
  const double *xx,
  const int n,
  const double x);
 
int locate_reg(
  const double *xx,
  const int n,
  const double x);
 
int locate_tbl(
  const float *xx,
  const int n,
  const double x);
 
void matrix_invert(
  gsl_matrix * a);
 
void matrix_product(
  const gsl_matrix * a,
  const gsl_vector * b,
  const int transpose,
  gsl_matrix * c);
 
size_t obs2y(
  const ctl_t * ctl,
  const obs_t * obs,
  gsl_vector * y,
  int *ida,
  int *ira);
 
void optimal_estimation(
  ret_t * ret,
  ctl_t * ctl,
  tbl_t * tbl,
  obs_t * obs_meas,
  obs_t * obs_i,
  atm_t * atm_apr,
  atm_t * atm_i,
  double *chisq);
 
int raytrace(
  const ctl_t * ctl,
  const atm_t * atm,
  obs_t * obs,
  los_t * los,
  const int ir);
 
void read_atm(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm,
  int profile);
 
void read_atm_asc(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
void read_atm_bin(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm);
 
void read_atm_nc(
  const char *filename,
  const ctl_t * ctl,
  atm_t * atm,
  int profile);
 
void read_ctl(
  int argc,
  char *argv[],
  ctl_t * ctl);
 
void read_matrix(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  gsl_matrix * matrix,
  int dataset);
 
void read_matrix_asc(
  const char *dirname,
  const char *filename,
  gsl_matrix * matrix);
 
void read_matrix_bin(
  const char *dirname,
  const char *filename,
  gsl_matrix * matrix);
 
void read_matrix_nc(
  const char *dirname,
  const char *filename,
  gsl_matrix * matrix,
  int dataset);
 
void read_obs(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  obs_t * obs,
  int profile);
 
void read_obs_asc(
  const char *filename,
  const ctl_t * ctl,
  obs_t * obs);
 
void read_obs_bin(
  const char *filename,
  const ctl_t * ctl,
  obs_t * obs);
 
void read_obs_nc(
  const char *filename,
  const ctl_t * ctl,
  obs_t * obs,
  const int profile);
 
double read_obs_rfm(
  const char *basename,
  const double z,
  const double *nu,
  const double *f,
  const int n);
 
void read_ret(
  int argc,
  char *argv[],
  const ctl_t * ctl,
  ret_t * ret);
 
void read_rfm_spec(
  const char *filename,
  double *nu,
  double *rad,
  int *npts);
 
void read_shape(
  const char *filename,
  double *x,
  double *y,
  int *n);
 
tbl_t *read_tbl(
  const ctl_t * ctl);
 
void read_tbl_asc(
  const ctl_t * ctl,
  tbl_t * tbl,
  const int id,
  const int ig);
 
void read_tbl_bin(
  const ctl_t * ctl,
  tbl_t * tbl,
  const int id,
  const int ig);
 
void read_tbl_nc_channel(
  const ctl_t * ctl,
  tbl_t * tbl,
  int id,
  int ig,
  int ncid);
 
double scan_ctl(
  int argc,
  char *argv[],
  const char *varname,
  const int arridx,
  const char *defvalue,
  char *value);
 
void set_cov_apr(
  const ret_t * ret,
  const ctl_t * ctl,
  const atm_t * atm,
  const int *iqa,
  const int *ipa,
  gsl_matrix * s_a);
 
void set_cov_meas(
  const ret_t * ret,
  const ctl_t * ctl,
  const obs_t * obs,
  gsl_vector * sig_noise,
  gsl_vector * sig_formod,
  gsl_vector * sig_eps_inv);
 
const char *shared_io_input_target(
  const ret_t * ret,
  const char *shared_file,
  const char *legacy_file,
  const char **dirname,
  int *profile);
 
const char *shared_io_output_target(
  const ret_t * ret,
  const char *shared_file,
  const char *legacy_file,
  const char **dirname,
  int *profile);
 
int shared_io_lock(
  const ret_t * ret);
 
void shared_io_unlock(
  int fd);
 
const char *shared_io_output_file(
  const ret_t * ret);
 
void tangent_point(
  const los_t * los,
  double *tpz,
  double *tplon,
  double *tplat);
 
void tbl_free(
  const ctl_t * ctl,
  tbl_t * tbl);
 
void tbl_pack(
  const tbl_t * tbl,
  int id,
  int ig,
  uint8_t * buf,
  size_t *bytes_used);
 
size_t tbl_packed_size(
  const tbl_t * tbl,
  int id,
  int ig);
 
size_t tbl_unpack(
  tbl_t * tbl,
  int id,
  int ig,
  const uint8_t * buf);
 
 
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_group(
  const char *name,
  const char *group,
  int output);
 
void write_atm(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  int profile);
 
void write_atm_asc(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_atm_bin(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_atm_nc(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm,
  int profile);
 
void write_atm_rfm(
  const char *filename,
  const ctl_t * ctl,
  const atm_t * atm);
 
void write_matrix(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const gsl_matrix * matrix,
  const atm_t * atm,
  const obs_t * obs,
  const char *rowspace,
  const char *colspace,
  const char *sort,
  int dataset);
 
void write_matrix_asc(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const gsl_matrix * matrix,
  const atm_t * atm,
  const obs_t * obs,
  const char *rowspace,
  const char *colspace,
  const char *sort);
 
void write_matrix_bin(
  const char *dirname,
  const char *filename,
  const gsl_matrix * matrix);
 
void write_matrix_nc(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const gsl_matrix * matrix,
  const atm_t * atm,
  const obs_t * obs,
  const char *rowspace,
  const char *colspace,
  const char *sort,
  int dataset);
 
void write_obs(
  const char *dirname,
  const char *filename,
  const ctl_t * ctl,
  const obs_t * obs,
  int profile);
 
void write_obs_asc(
  const char *filename,
  const ctl_t * ctl,
  const obs_t * obs);
 
void write_obs_bin(
  const char *filename,
  const ctl_t * ctl,
  const obs_t * obs);
 
void write_obs_nc(
  const char *filename,
  const ctl_t * ctl,
  const obs_t * obs,
  const int profile);
 
void write_shape(
  const char *filename,
  const double *x,
  const double *y,
  const int n);
 
void write_stddev(
  const char *quantity,
  const ret_t * ret,
  const ctl_t * ctl,
  const atm_t * atm,
  const gsl_matrix * s);
 
void write_tbl(
  const ctl_t * ctl,
  const tbl_t * tbl);
 
void write_tbl_asc(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const int id,
  const int ig);
 
void write_tbl_bin(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const int id,
  const int ig);
 
void write_tbl_nc(
  const ctl_t * ctl,
  const tbl_t * tbl,
  const int id,
  const int ig);
 
void x2atm(
  const ctl_t * ctl,
  const gsl_vector * x,
  atm_t * atm);
 
void x2atm_help(
  double *value,
  const gsl_vector * x,
  size_t *n);
 
void y2obs(
  const ctl_t * ctl,
  const gsl_vector * y,
  obs_t * obs);
 
#endif