tropo.c

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
*/
 
#include "mptrac.h"
 
/* ------------------------------------------------------------
   Functions...
   ------------------------------------------------------------ */
 
void usage(
  void);
 
/* ------------------------------------------------------------
   Main...
   ------------------------------------------------------------ */
 
int main(
  int argc,
  char *argv[]) {
 
  ctl_t ctl;
 
  clim_t *clim;
 
  met_t *met;
 
  dd_t *dd;
 
  static double ps[EX * EY], pt[EX * EY], qt[EX * EY], o3t[EX * EY],
    zs[EX * EY], zt[EX * EY], tt[EX * EY], lon, lons[EX], lat, lats[EY];
 
  static int init, nx, ny, nt, ncid, varid, dims[3];
 
  static size_t count[10], start[10];
 
  /* Allocate... */
  ALLOC(clim, clim_t, 1);
  ALLOC(met, met_t, 1);
  ALLOC(dd, dd_t, 1);
 
  /* Print usage information... */
  USAGE;
 
  /* Check arguments... */
  if (argc < 4)
    ERRMSG("Missing or invalid command-line arguments.\n\n"
           "Usage: tropo <ctl> <tropo.nc> <met0> [<met1> ...]\n\n"
           "Use -h for full help.");
 
  /* Read control parameters... */
  mptrac_read_ctl(argv[1], argc, argv, &ctl);
  double lon0 = scan_ctl(argv[1], argc, argv, "TROPO_LON0", -1, "-180", NULL);
  double lon1 = scan_ctl(argv[1], argc, argv, "TROPO_LON1", -1, "180", NULL);
  double dlon = scan_ctl(argv[1], argc, argv, "TROPO_DLON", -1, "-999", NULL);
  double lat0 = scan_ctl(argv[1], argc, argv, "TROPO_LAT0", -1, "-90", NULL);
  double lat1 = scan_ctl(argv[1], argc, argv, "TROPO_LAT1", -1, "90", NULL);
  double dlat = scan_ctl(argv[1], argc, argv, "TROPO_DLAT", -1, "-999", NULL);
  int h2o = (int) scan_ctl(argv[1], argc, argv, "TROPO_H2O", -1, "1", NULL);
  int o3 = (int) scan_ctl(argv[1], argc, argv, "TROPO_O3", -1, "1", NULL);
 
  /* Read climatological data... */
  mptrac_read_clim(&ctl, clim);
 
  /* Loop over files... */
  for (int i = 3; i < argc; i++) {
 
    /* Set control parameters... */
    ctl.met_tropo = 0;
 
    /* Read meteorological data... */
    if (!mptrac_read_met(argv[i], &ctl, clim, met, dd))
      continue;
 
    /* Set horizontal grid... */
    if (!init) {
      init = 1;
 
      /* Get grid... */
      if (dlon <= 0)
        dlon = fabs(met->lon[1] - met->lon[0]);
      if (dlat <= 0)
        dlat = fabs(met->lat[1] - met->lat[0]);
      if (lon0 < -360 && lon1 > 360) {
        lon0 = gsl_stats_min(met->lon, 1, (size_t) met->nx);
        lon1 = gsl_stats_max(met->lon, 1, (size_t) met->nx);
      }
      nx = ny = 0;
      for (lon = lon0; lon <= lon1 + 0.001; lon += dlon) {
        lons[nx] = round(lon * 1e3) / 1e3;
        if ((++nx) >= EX)
          ERRMSG("Too many longitudes!");
      }
      if (lat0 < -90 && lat1 > 90) {
        lat0 = gsl_stats_min(met->lat, 1, (size_t) met->ny);
        lat1 = gsl_stats_max(met->lat, 1, (size_t) met->ny);
      }
      for (lat = lat0; lat <= lat1 + 0.001; lat += dlat) {
        lats[ny] = round(lat * 1e3) / 1e3;
        if ((++ny) >= EY)
          ERRMSG("Too many latitudes!");
      }
 
      /* Create netCDF file... */
      LOG(1, "Write tropopause data file: %s", argv[2]);
      NC(nc_create(argv[2], NC_NETCDF4, &ncid));
 
      /* Create dimensions... */
      NC(nc_def_dim(ncid, "time", (size_t) NC_UNLIMITED, &dims[0]));
      NC(nc_def_dim(ncid, "lat", (size_t) ny, &dims[1]));
      NC(nc_def_dim(ncid, "lon", (size_t) nx, &dims[2]));
 
      /* Create variables... */
      NC_DEF_VAR("time", NC_DOUBLE, 1, &dims[0], "time",
                 "seconds since 2000-01-01 00:00:00 UTC", 0, 0);
      NC_DEF_VAR("lat", NC_DOUBLE, 1, &dims[1], "latitude", "degrees_north",
                 0, 0);
      NC_DEF_VAR("lon", NC_DOUBLE, 1, &dims[2], "longitude", "degrees_east",
                 0, 0);
 
      NC_DEF_VAR("clp_z", NC_FLOAT, 3, &dims[0], "cold point height", "km", 0,
                 0);
      NC_DEF_VAR("clp_p", NC_FLOAT, 3, &dims[0], "cold point pressure", "hPa",
                 0, 0);
      NC_DEF_VAR("clp_t", NC_FLOAT, 3, &dims[0], "cold point temperature",
                 "K", 0, 0);
      if (h2o)
        NC_DEF_VAR("clp_q", NC_FLOAT, 3, &dims[0], "cold point water vapor",
                   "ppv", 0, 0);
      if (o3)
        NC_DEF_VAR("clp_o3", NC_FLOAT, 3, &dims[0], "cold point ozone",
                   "ppv", 0, 0);
 
      NC_DEF_VAR("dyn_z", NC_FLOAT, 3, &dims[0],
                 "dynamical tropopause height", "km", 0, 0);
      NC_DEF_VAR("dyn_p", NC_FLOAT, 3, &dims[0],
                 "dynamical tropopause pressure", "hPa", 0, 0);
      NC_DEF_VAR("dyn_t", NC_FLOAT, 3, &dims[0],
                 "dynamical tropopause temperature", "K", 0, 0);
      if (h2o)
        NC_DEF_VAR("dyn_q", NC_FLOAT, 3, &dims[0],
                   "dynamical tropopause water vapor", "ppv", 0, 0);
      if (o3)
        NC_DEF_VAR("dyn_o3", NC_FLOAT, 3, &dims[0],
                   "dynamical tropopause ozone", "ppv", 0, 0);
 
      NC_DEF_VAR("wmo_1st_z", NC_FLOAT, 3, &dims[0],
                 "WMO 1st tropopause height", "km", 0, 0);
      NC_DEF_VAR("wmo_1st_p", NC_FLOAT, 3, &dims[0],
                 "WMO 1st tropopause pressure", "hPa", 0, 0);
      NC_DEF_VAR("wmo_1st_t", NC_FLOAT, 3, &dims[0],
                 "WMO 1st tropopause temperature", "K", 0, 0);
      if (h2o)
        NC_DEF_VAR("wmo_1st_q", NC_FLOAT, 3, &dims[0],
                   "WMO 1st tropopause water vapor", "ppv", 0, 0);
      if (o3)
        NC_DEF_VAR("wmo_1st_o3", NC_FLOAT, 3, &dims[0],
                   "WMO 1st tropopause ozone", "ppv", 0, 0);
 
      NC_DEF_VAR("wmo_2nd_z", NC_FLOAT, 3, &dims[0],
                 "WMO 2nd tropopause height", "km", 0, 0);
      NC_DEF_VAR("wmo_2nd_p", NC_FLOAT, 3, &dims[0],
                 "WMO 2nd tropopause pressure", "hPa", 0, 0);
      NC_DEF_VAR("wmo_2nd_t", NC_FLOAT, 3, &dims[0],
                 "WMO 2nd tropopause temperature", "K", 0, 0);
      if (h2o)
        NC_DEF_VAR("wmo_2nd_q", NC_FLOAT, 3, &dims[0],
                   "WMO 2nd tropopause water vapor", "ppv", 0, 0);
      if (o3)
        NC_DEF_VAR("wmo_2nd_o3", NC_FLOAT, 3, &dims[0],
                   "WMO 2nd tropopause ozone", "ppv", 0, 0);
 
      NC_DEF_VAR("ps", NC_FLOAT, 3, &dims[0], "surface pressure", "hPa", 0,
                 0);
      NC_DEF_VAR("zs", NC_FLOAT, 3, &dims[0], "surface height", "km", 0, 0);
 
      /* End definition... */
      NC(nc_enddef(ncid));
 
      /* Write longitude and latitude... */
      NC_PUT_DOUBLE("lat", lats, 0);
      NC_PUT_DOUBLE("lon", lons, 0);
    }
 
    /* Write time... */
    start[0] = (size_t) nt;
    count[0] = 1;
    start[1] = 0;
    count[1] = (size_t) ny;
    start[2] = 0;
    count[2] = (size_t) nx;
    NC_PUT_DOUBLE("time", &met->time, 1);
 
    /* Get cold point... */
    get_tropo(2, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
              zs);
    NC_PUT_DOUBLE("clp_z", zt, 1);
    NC_PUT_DOUBLE("clp_p", pt, 1);
    NC_PUT_DOUBLE("clp_t", tt, 1);
    if (h2o)
      NC_PUT_DOUBLE("clp_q", qt, 1);
    if (o3)
      NC_PUT_DOUBLE("clp_o3", o3t, 1);
 
    /* Get dynamical tropopause... */
    get_tropo(5, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
              zs);
    NC_PUT_DOUBLE("dyn_z", zt, 1);
    NC_PUT_DOUBLE("dyn_p", pt, 1);
    NC_PUT_DOUBLE("dyn_t", tt, 1);
    if (h2o)
      NC_PUT_DOUBLE("dyn_q", qt, 1);
    if (o3)
      NC_PUT_DOUBLE("dyn_o3", o3t, 1);
 
    /* Get WMO 1st tropopause... */
    get_tropo(3, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
              zs);
    NC_PUT_DOUBLE("wmo_1st_z", zt, 1);
    NC_PUT_DOUBLE("wmo_1st_p", pt, 1);
    NC_PUT_DOUBLE("wmo_1st_t", tt, 1);
    if (h2o)
      NC_PUT_DOUBLE("wmo_1st_q", qt, 1);
    if (o3)
      NC_PUT_DOUBLE("wmo_1st_o3", o3t, 1);
 
    /* Get WMO 2nd tropopause... */
    get_tropo(4, &ctl, clim, met, lons, nx, lats, ny, pt, zt, tt, qt, o3t, ps,
              zs);
    NC_PUT_DOUBLE("wmo_2nd_z", zt, 1);
    NC_PUT_DOUBLE("wmo_2nd_p", pt, 1);
    NC_PUT_DOUBLE("wmo_2nd_t", tt, 1);
    if (h2o)
      NC_PUT_DOUBLE("wmo_2nd_q", qt, 1);
    if (o3)
      NC_PUT_DOUBLE("wmo_2nd_o3", o3t, 1);
 
    /* Write surface data... */
    NC_PUT_DOUBLE("ps", ps, 1);
    NC_PUT_DOUBLE("zs", zs, 1);
 
    /* Increment time step counter... */
    nt++;
  }
 
  /* Close file... */
  NC(nc_close(ncid));
 
  /* Free... */
  free(clim);
  free(met);
  free(dd);
 
  return EXIT_SUCCESS;
}
 
/*****************************************************************************/
 
void usage(
  void) {
 
  printf("\nMPTRAC tropo tool.\n\n");
  printf("Create tropopause data sets from meteorological data.\n");
  printf("\n");
  printf("Usage:\n");
  printf("  tropo <ctl> <tropo.nc> <met0> [<met1> ...]\n");
  printf("\n");
  printf("Arguments:\n");
  printf("  <ctl>       Control file.\n");
  printf("  <tropo.nc>  Output tropopause netCDF file.\n");
  printf("  <met*>      Meteorological input files.\n");
  printf("\nFurther information:\n");
  printf("  Manual: https://slcs-jsc.github.io/mptrac/\n");
}