trac.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"
 
#ifdef KPP
#include "kpp_chem.h"
#endif
 
/* ------------------------------------------------------------
   Functions...
   ------------------------------------------------------------ */
 
void usage(
  void);
 
/* ------------------------------------------------------------
   Main...
   ------------------------------------------------------------ */
 
int main(
  int argc,
  char *argv[]) {
 
  ctl_t *ctl;
 
  atm_t *atm;
 
  cache_t *cache;
 
  clim_t *clim;
 
  depo_t *depo;
 
  met_t *met0, *met1;
 
  dd_t *dd;
 
  FILE *dirlist;
 
  char dirname[LEN], filename[2 * LEN];
 
  int ntask = -1, rank = 0, size = 1;
 
  /* Print usage information... */
  USAGE;
 
  /* Initialize MPI... */
#ifdef MPI
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
#ifdef _OPENACC
  const int num_devices = acc_get_num_devices(acc_device_nvidia);
  if (num_devices <= 0)
    ERRMSG("Not running on a GPU device!");
  acc_set_device_num(rank % num_devices, acc_device_nvidia);
#endif
#endif
 
  /* Check arguments... */
  if (argc < 4)
    ERRMSG("Missing or invalid command-line arguments.\n\n"
           "Usage: trac <dirlist> <ctl> <atm_in> [KEY VALUE ...]\n\n"
           "Use -h for full help.");
 
  /* Open directory list... */
  if (!(dirlist = fopen(argv[1], "r")))
    ERRMSG("Cannot open directory list!");
 
  /* Loop over directories... */
  while (fscanf(dirlist, "%4999s", dirname) != EOF) {
 
    /* MPI parallelization... */
    if ((++ntask) % size != rank)
      continue;
 
    /* Write info... */
    LOG(1, "Parallelization: ntask= %d | rank= %d | size= %d",
        ntask, rank, size);
 
    /* ------------------------------------------------------------
       Initialize model run...
       ------------------------------------------------------------ */
 
    /* Allocate memory... */
    mptrac_alloc(&ctl, &cache, &clim, &met0, &met1, &atm, &depo, &dd);
 
    /* Read control parameters... */
    sprintf(filename, "%s/%s", dirname, argv[2]);
    mptrac_read_ctl(filename, argc, argv, ctl);
 
    /* This tool implements the domain decomposition... */
    ctl->dd_supported = 1;
 
    /* Read climatological data... */
    mptrac_read_clim(ctl, clim);
 
    /* Read atmospheric data... */
    sprintf(filename, "%s/%s", dirname, argv[3]);
    if (!mptrac_read_atm(filename, ctl, atm))
      ERRMSG("Cannot open file!");
 
    /* Initialize MPTRAC... */
    mptrac_init(ctl, cache, clim, atm, depo, ntask);
#ifdef DD
    int peak_np = atm->np;
#endif
 
    /* ------------------------------------------------------------
       Loop over timesteps...
       ------------------------------------------------------------ */
 
    /* Loop over timesteps... */
    for (double t = ctl->t_start;
         ctl->direction * (t - ctl->t_stop) < ctl->dt_mod;
         t += ctl->direction * ctl->dt_mod) {
 
      /* Adjust length of final time step... */
      if (ctl->direction * (t - ctl->t_stop) > 0)
        t = ctl->t_stop;
 
      /* Get meteo data... */
      mptrac_get_met(ctl, clim, t, &met0, &met1, dd);
 
      /* Recommend time step from horizontal grid spacing... */
      if (t == ctl->t_start) {
        double dt_rec;
        if (met0->coord_type == 0)
          dt_rec = fabs(met0->lon[1] - met0->lon[0]) * 111132. / 150.;
        else {
          double spacing = INFINITY;
          for (int ix = 1; ix < met0->nx; ix++) {
            const double dx = fabs(met0->lon[ix] - met0->lon[ix - 1]);
            if (!isfinite(dx) || dx <= 0) {
              spacing = NAN;
              break;
            }
            spacing = MIN(spacing, dx);
          }
          if (isfinite(spacing))
            for (int iy = 1; iy < met0->ny; iy++) {
              const double dy = fabs(met0->lat[iy] - met0->lat[iy - 1]);
              if (!isfinite(dy) || dy <= 0) {
                spacing = NAN;
                break;
              }
              spacing = MIN(spacing, dy);
            }
          dt_rec = spacing / 50.;
        }
        if (!isfinite(dt_rec) || dt_rec <= 0) {
          WARN
            ("Cannot recommend DT_MOD: invalid horizontal meteo grid spacing!");
        } else if (ctl->dt_mod > dt_rec) {
          WARN
            ("DT_MOD=%g s may be too large for accurate trajectory integration; "
             "recommended maximum from horizontal meteo grid resolution is %g s.",
             ctl->dt_mod, dt_rec);
        }
      }
#ifdef DD
      /* Set-up domain decomposition... */
      if ((t == ctl->t_start) && ctl->dd)
        dd_init(ctl, dd, atm);
#else
      /* Catch misconfiguration... */
      if ((ctl->dd)
          || (ctl->dd_subdomains_zonal * ctl->dd_subdomains_meridional > 1))
        ERRMSG
          ("Model configured for domain decomposition, but not compiled for it!");
#endif
 
      /* Run a single time step... */
      mptrac_run_timestep(ctl, cache, clim, &met0, &met1, atm, depo, t, dd);
#ifdef DD
      if (ctl->dd && atm->np > peak_np)
        peak_np = atm->np;
#endif
 
      /* Write output... */
      mptrac_write_output(dirname, ctl, met0, met1, atm, depo, t);
    }
 
    /* ------------------------------------------------------------
       Finalize model run...
       ------------------------------------------------------------ */
 
    /* Flush output buffer... */
    fflush(NULL);
 
    /* Report problem size... */
    LOG(1, "SIZE_NP = %d", atm->np);
#ifdef DD
    if (ctl->dd)
      LOG(1, "SIZE_NP_PEAK = %d", peak_np);
#endif
    LOG(1, "SIZE_MPI_TASKS = %d", size);
    LOG(1, "SIZE_OMP_THREADS = %d", omp_get_max_threads());
 
    /* Report memory usage... */
    LOG(1, "MEMORY_ATM = %g MByte", sizeof(atm_t) / 1024. / 1024.);
    LOG(1, "MEMORY_CACHE = %g MByte", sizeof(cache_t) / 1024. / 1024.);
    LOG(1, "MEMORY_DEPO = %g MByte", sizeof(depo_t) / 1024. / 1024.);
    LOG(1, "MEMORY_CLIM = %g MByte", sizeof(clim_t) / 1024. / 1024.);
    LOG(1, "MEMORY_METEO = %g MByte", sizeof(met_t) / 1024. / 1024.);
 
    /* Free memory... */
    mptrac_free(ctl, cache, clim, met0, met1, atm, depo, dd);
 
    /* Report timers... */
    PRINT_TIMERS;
  }
 
  /* Finalize MPI... */
#ifdef MPI
  MPI_Finalize();
#endif
 
  return EXIT_SUCCESS;
}
 
/*****************************************************************************/
 
void usage(
  void) {
 
  printf("\nMPTRAC trac tool.\n\n");
  printf("Run forward or backward trajectory calculations.\n");
  printf("\n");
  printf("Usage:\n");
  printf("  trac <dirlist> <ctl> <atm_in> [KEY VALUE ...]\n");
  printf("\n");
  printf("Arguments:\n");
  printf("  <dirlist>  Text file containing work directories to process.\n");
  printf("  <ctl>      Control file name relative to each work directory.\n");
  printf
    ("  <atm_in>   Atmospheric input file name relative to each work directory.\n");
  printf("  [KEY VALUE]  Optional control parameters.\n");
  printf("\nFurther information:\n");
  printf("  Manual: https://slcs-jsc.github.io/mptrac/\n");
}