retrieval.c File Reference

Retrieval processor for IASI. More...

#include <mpi.h>
#include <omp.h>
#include <netcdf.h>
#include "jurassic.h"

Go to the source code of this file.

Data Structures
struct	ncd_t
	Buffer for netCDF data. More...

Macros
#define	NC(cmd)
	Execute netCDF library command and check result. More...
#define	L1_NCHAN   33
	Number of IASI radiance channels (don't change). More...
#define	L1_NTRACK   1800
	Maximum along-track size of IASI radiance granule (don't change). More...
#define	L1_NXTRACK   60
	Across-track size of IASI radiance granule (don't change). More...
#define	L2_NLAY   27
	Number of IASI pressure layers (don't change). More...
#define	L2_NTRACK   1800
	Maximum along-track size of IASI retrieval granule (don't change). More...
#define	L2_NXTRACK   60
	Across-track size of IASI retrieval granule (don't change). More...

Functions
void	add_var (int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
	Create variable in netCDF file. More...
void	buffer_nc (atm_t *atm, double chisq, ncd_t *ncd, int track, int xtrack, int np0, int np1)
	Buffer netCDF data. More...
void	init_l2 (ncd_t *ncd, int track, int xtrack, ctl_t *ctl, atm_t *atm)
	Initialize with IASI Level-2 data. More...
void	read_nc (char *filename, ncd_t *ncd)
	Read netCDF file. More...
void	write_nc (char *filename, ncd_t *ncd)
	Write to netCDF file... More...
int	main (int argc, char *argv[])

Detailed Description

Retrieval processor for IASI.

Definition in file retrieval.c.

Macro Definition Documentation

NC

#define NC

(

cmd

)

Value:

  {                                  \
  int nc_result=(cmd);                               \
  if(nc_result!=NC_NOERR)                            \
    ERRMSG("%s", nc_strerror(nc_result));            \
}

Execute netCDF library command and check result.

Definition at line 36 of file retrieval.c.

L1_NCHAN

#define L1_NCHAN   33

Number of IASI radiance channels (don't change).

Definition at line 47 of file retrieval.c.

L1_NTRACK

#define L1_NTRACK   1800

Maximum along-track size of IASI radiance granule (don't change).

Definition at line 50 of file retrieval.c.

L1_NXTRACK

#define L1_NXTRACK   60

Across-track size of IASI radiance granule (don't change).

Definition at line 53 of file retrieval.c.

L2_NLAY

#define L2_NLAY   27

Number of IASI pressure layers (don't change).

Definition at line 56 of file retrieval.c.

L2_NTRACK

#define L2_NTRACK   1800

Maximum along-track size of IASI retrieval granule (don't change).

Definition at line 59 of file retrieval.c.

L2_NXTRACK

#define L2_NXTRACK   60

Across-track size of IASI retrieval granule (don't change).

Definition at line 62 of file retrieval.c.

Function Documentation

add_var()

void add_var	(	int	ncid,
		const char *	varname,
		const char *	unit,
		const char *	longname,
		int	type,
		int	dimid[],
		int *	varid,
		int	ndims
	)

Create variable in netCDF file.

Add variable to netCDF file.

Definition at line 400 of file retrieval.c.

             {
 
  /* Check if variable exists... */
  if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
 
    /* Define variable... */
    NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
 
    /* Set long name... */
    NC(nc_put_att_text
       (ncid, *varid, "long_name", strlen(longname), longname));
 
    /* Set units... */
    NC(nc_put_att_text(ncid, *varid, "units", strlen(unit), unit));
  }
}

buffer_nc()

void buffer_nc	(	atm_t *	atm,
		double	chisq,
		ncd_t *	ncd,
		int	track,
		int	xtrack,
		int	np0,
		int	np1
	)

Buffer netCDF data.

Definition at line 427 of file retrieval.c.

           {
 
  int ip;
 
  /* Set number of data points... */
  ncd->np = np1 - np0 + 1;
 
  /* Save retrieval data... */
  ncd->ret_chisq[track * L1_NXTRACK + xtrack] = (float) chisq;
  ncd->ret_p[track * L1_NXTRACK + xtrack] = (float) atm->p[np0];
  for (ip = np0; ip <= np1; ip++) {
    ncd->ret_z[ip - np0] = (float) atm->z[ip];
    ncd->ret_t[(track * L1_NXTRACK + xtrack) * ncd->np + ip - np0] =
      (gsl_finite(chisq) ? (float) atm->t[ip] : GSL_NAN);
  }
}

init_l2()

void init_l2	(	ncd_t *	ncd,
		int	track,
		int	xtrack,
		ctl_t *	ctl,
		atm_t *	atm
	)

Initialize with IASI Level-2 data.

Definition at line 453 of file retrieval.c.

              {
 
  static atm_t atm_iasi;
 
  double k[NW], p, q[NG], t, w, zmax = 0, zmin = 1000;
 
  int ip, lay;
 
  /* Store IASI data in atmospheric data struct... */
  atm_iasi.np = 0;
  for (lay = 0; lay < L2_NLAY; lay++)
    if (gsl_finite(ncd->l2_z[track][xtrack][lay])
        && ncd->l2_z[track][xtrack][lay] <= 60.) {
      atm_iasi.z[atm_iasi.np] = ncd->l2_z[track][xtrack][lay];
      atm_iasi.p[atm_iasi.np] = ncd->l2_p[lay];
      atm_iasi.t[atm_iasi.np] = ncd->l2_t[track][xtrack][lay];
      if ((++atm_iasi.np) > NP)
        ERRMSG("Too many layers!");
    }
 
  /* Check number of levels... */
  if (atm_iasi.np < 2)
    return;
 
  /* Get height range of IASI data... */
  for (ip = 0; ip < atm_iasi.np; ip++) {
    zmax = GSL_MAX(zmax, atm_iasi.z[ip]);
    zmin = GSL_MIN(zmin, atm_iasi.z[ip]);
  }
 
  /* Merge IASI data... */
  for (ip = 0; ip < atm->np; ip++) {
 
    /* Interpolate IASI data... */
    intpol_atm(ctl, &atm_iasi, atm->z[ip], &p, &t, q, k);
 
    /* Weighting factor... */
    w = 1;
    if (atm->z[ip] > zmax)
      w = GSL_MAX(1 - (atm->z[ip] - zmax) / 50, 0);
    if (atm->z[ip] < zmin)
      w = GSL_MAX(1 - (zmin - atm->z[ip]) / 50, 0);
 
    /* Merge... */
    atm->t[ip] = w * t + (1 - w) * atm->t[ip];
    atm->p[ip] = w * p + (1 - w) * atm->p[ip];
  }
}

read_nc()

void read_nc	(	char *	filename,
		ncd_t *	ncd
	)

Read netCDF file.

Definition at line 509 of file retrieval.c.

              {
 
  int dimid, varid;
 
  size_t len;
 
  /* Open netCDF file... */
  printf("Read netCDF file: %s\n", filename);
  NC(nc_open(filename, NC_WRITE, &ncd->ncid));
 
  /* Read number of tracks... */
  NC(nc_inq_dimid(ncd->ncid, "L1_NTRACK", &dimid));
  NC(nc_inq_dimlen(ncd->ncid, dimid, &len));
  ncd->ntrack = (int) len;
 
  /* Read Level-1 data... */
  NC(nc_inq_varid(ncd->ncid, "l1_time", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_time[0]));
  NC(nc_inq_varid(ncd->ncid, "l1_lon", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_lon[0]));
  NC(nc_inq_varid(ncd->ncid, "l1_lat", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_lat[0]));
  NC(nc_inq_varid(ncd->ncid, "l1_sat_z", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_z));
  NC(nc_inq_varid(ncd->ncid, "l1_sat_lon", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_lon));
  NC(nc_inq_varid(ncd->ncid, "l1_sat_lat", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_sat_lat));
  NC(nc_inq_varid(ncd->ncid, "l1_nu", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l1_nu));
  NC(nc_inq_varid(ncd->ncid, "l1_rad", &varid));
  NC(nc_get_var_float(ncd->ncid, varid, ncd->l1_rad[0][0]));
 
  /* Read Level-2 data... */
  NC(nc_inq_varid(ncd->ncid, "l2_z", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_z[0][0]));
  NC(nc_inq_varid(ncd->ncid, "l2_press", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_p));
  NC(nc_inq_varid(ncd->ncid, "l2_temp", &varid));
  NC(nc_get_var_double(ncd->ncid, varid, ncd->l2_t[0][0]));
}

write_nc()

void write_nc	(	char *	filename,
		ncd_t *	ncd
	)

Write to netCDF file...

Definition at line 555 of file retrieval.c.

              {
 
  int dimid[10], c_id, p_id, t_id, z_id;
 
  /* Create netCDF file... */
  printf("Write netCDF file: %s\n", filename);
 
  /* Read existing dimensions... */
  NC(nc_inq_dimid(ncd->ncid, "L1_NTRACK", &dimid[0]));
  NC(nc_inq_dimid(ncd->ncid, "L1_NXTRACK", &dimid[1]));
 
  /* Set define mode... */
  NC(nc_redef(ncd->ncid));
 
  /* Set new dimensions... */
  if (nc_inq_dimid(ncd->ncid, "RET_NP", &dimid[2]) != NC_NOERR)
    NC(nc_def_dim(ncd->ncid, "RET_NP", (size_t) ncd->np, &dimid[2]));
 
  /* Set new variables... */
  add_var(ncd->ncid, "ret_z", "km", "altitude", NC_FLOAT, &dimid[2], &z_id,
          1);
  add_var(ncd->ncid, "ret_press", "hPa", "pressure", NC_FLOAT, dimid, &p_id,
          2);
  add_var(ncd->ncid, "ret_temp", "K", "temperature", NC_FLOAT, dimid, &t_id,
          3);
  add_var(ncd->ncid, "ret_chisq", "1", "chi^2 value of fit", NC_FLOAT, dimid,
          &c_id, 2);
 
  /* Leave define mode... */
  NC(nc_enddef(ncd->ncid));
 
  /* Write data... */
  NC(nc_put_var_float(ncd->ncid, z_id, ncd->ret_z));
  NC(nc_put_var_float(ncd->ncid, p_id, ncd->ret_p));
  NC(nc_put_var_float(ncd->ncid, t_id, ncd->ret_t));
  NC(nc_put_var_float(ncd->ncid, c_id, ncd->ret_chisq));
 
  /* Close netCDF file... */
  NC(nc_close(ncd->ncid));
}

Here is the call graph for this function:

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 174 of file retrieval.c.

                {
 
  static ctl_t ctl;
  static atm_t atm_apr, atm_clim, atm_i;
  static obs_t obs_i, obs_meas;
  static ncd_t ncd;
  static ret_t ret;
 
  FILE *in;
 
  char filename[LEN], filename2[LEN];
 
  double chisq, sza_thresh, z[NP], t0;
 
  int channel[ND], i, id, ip, iz, nz, ntask = -1, rank, size,
    np0, np1, track, track0, track1, xtrack, xtrack0, xtrack1,
    task0, task1, debug;
 
  /* ------------------------------------------------------------
     Init...
     ------------------------------------------------------------ */
 
  /* MPI... */
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
 
  /* Measure CPU time... */
  TIMER("total", 1);
 
  /* Check arguments... */
  if (argc < 3)
    ERRMSG("Give parameters: <ctl> <filelist>");
 
  /* Read control parameters... */
  read_ctl(argc, argv, &ctl);
  read_ret(argc, argv, &ctl, &ret);
  debug = (int) scan_ctl(argc, argv, "DEBUG", -1, "1", NULL);
 
  /* Initialize look-up tables... */
  tbl_t *tbl = read_tbl(&ctl);
 
  /* Read retrieval grid... */
  nz = (int) scan_ctl(argc, argv, "NZ", -1, "", NULL);
  if (nz > NP)
    ERRMSG("Too many altitudes!");
  for (iz = 0; iz < nz; iz++)
    z[iz] = scan_ctl(argc, argv, "Z", iz, "", NULL);
 
  /* Read task range... */
  task0 = (int) scan_ctl(argc, argv, "TASK_MIN", -1, "0", NULL);
  task1 = (int) scan_ctl(argc, argv, "TASK_MAX", -1, "99999", NULL);
 
  /* Read track range... */
  track0 = (int) scan_ctl(argc, argv, "TRACK_MIN", -1, "0", NULL);
  track1 = (int) scan_ctl(argc, argv, "TRACK_MAX", -1, "99999", NULL);
 
  /* Read xtrack range... */
  xtrack0 = (int) scan_ctl(argc, argv, "XTRACK_MIN", -1, "0", NULL);
  xtrack1 = (int) scan_ctl(argc, argv, "XTRACK_MAX", -1, "59", NULL);
 
  /* Read height range... */
  np0 = (int) scan_ctl(argc, argv, "NP_MIN", -1, "0", NULL);
  np1 = (int) scan_ctl(argc, argv, "NP_MAX", -1, "100", NULL);
  np1 = GSL_MIN(np1, nz - 1);
 
  /* SZA threshold... */
  sza_thresh = scan_ctl(argc, argv, "SZA", -1, "96", NULL);
 
  /* ------------------------------------------------------------
     Distribute granules...
     ------------------------------------------------------------ */
 
  /* Open filelist... */
  printf("Read filelist: %s\n", argv[2]);
  if (!(in = fopen(argv[2], "r")))
    ERRMSG("Cannot open filelist!");
 
  /* Loop over netCDF files... */
  while (fscanf(in, "%s", filename) != EOF) {
 
    /* Distribute files with MPI... */
    if ((++ntask) % size != rank)
      continue;
 
    /* Check task range... */
    if (ntask < task0 || ntask > task1)
      continue;
 
    /* Write info... */
    printf("Retrieve file %s on rank %d of %d (with %d threads)...\n",
           filename, rank + 1, size, omp_get_max_threads());
 
    /* ------------------------------------------------------------
       Initialize retrieval...
       ------------------------------------------------------------ */
 
    /* Read netCDF file... */
    read_nc(filename, &ncd);
 
    /* Adjust number of tracks... */
    if (track1 >= ncd.ntrack)
      track1 = ncd.ntrack - 1;
 
    /* Identify radiance channels... */
    for (id = 0; id < ctl.nd; id++) {
      channel[id] = -999;
      for (i = 0; i < L1_NCHAN; i++)
        if (fabs(ctl.nu[id] - ncd.l1_nu[i]) < 0.1)
          channel[id] = i;
      if (channel[id] < 0)
        ERRMSG("Cannot identify radiance channel!");
    }
 
    /* Set climatological data for center of granule... */
    atm_clim.np = nz;
    for (iz = 0; iz < nz; iz++)
      atm_clim.z[iz] = z[iz];
    climatology(&ctl, &atm_clim);
 
    /* ------------------------------------------------------------
       Retrieval...
       ------------------------------------------------------------ */
 
    /* Loop over swaths... */
    for (track = track0; track <= track1; track++) {
 
      /* Loop over scan... */
      for (xtrack = xtrack0; xtrack <= xtrack1; xtrack++) {
 
        /* Init timer... */
        t0 = omp_get_wtime();
 
        /* Store observation data... */
        obs_meas.nr = 1;
        obs_meas.time[0] = ncd.l1_time[track][xtrack];
        obs_meas.obsz[0] = ncd.l1_sat_z[track];
        obs_meas.obslon[0] = ncd.l1_sat_lon[track];
        obs_meas.obslat[0] = ncd.l1_sat_lat[track];
        obs_meas.vplon[0] = ncd.l1_lon[track][xtrack];
        obs_meas.vplat[0] = ncd.l1_lat[track][xtrack];
        for (id = 0; id < ctl.nd; id++)
          obs_meas.rad[id][0] = ncd.l1_rad[track][xtrack][channel[id]];
 
        /* Flag out 4 micron channels for daytime measurements... */
        if (RAD2DEG(acos(cos_sza(obs_meas.time[0], obs_meas.obslon[0],
                                 obs_meas.obslat[0]))) < sza_thresh)
          for (id = 0; id < ctl.nd; id++)
            if (ctl.nu[id] >= 2000)
              obs_meas.rad[id][0] = GSL_NAN;
 
        /* Prepare atmospheric data... */
        copy_atm(&ctl, &atm_apr, &atm_clim, 0);
        for (ip = 0; ip < atm_apr.np; ip++) {
          atm_apr.time[ip] = obs_meas.time[0];
          atm_apr.lon[ip] = obs_meas.vplon[0];
          atm_apr.lat[ip] = obs_meas.vplat[0];
        }
 
        /* Merge Level-2 data... */
        init_l2(&ncd, track, xtrack, &ctl, &atm_apr);
 
        /* Retrieval... */
        optimal_estimation(&ret, &ctl, tbl, &obs_meas, &obs_i,
                           &atm_apr, &atm_i, &chisq);
 
        /* Buffer results... */
        buffer_nc(&atm_i, chisq, &ncd, track, xtrack, np0, np1);
 
        /* Write debug information... */
        if (debug >= 1)
          printf
            ("  task= %4d | track= %5d | xtrack= %3d | chi^2= %8.3f | time= %8.3f s\n",
             ntask, track, xtrack, chisq, omp_get_wtime() - t0);
        if (debug >= 2) {
          sprintf(filename2, "atm_apr_%d_%d_%d.tab", ntask, track, xtrack);
          write_atm(NULL, filename2, &ctl, &atm_apr);
          sprintf(filename2, "atm_i_%d_%d_%d.tab", ntask, track, xtrack);
          write_atm(NULL, filename2, &ctl, &atm_i);
          sprintf(filename2, "obs_meas_%d_%d_%d.tab", ntask, track, xtrack);
          write_obs(NULL, filename2, &ctl, &obs_meas);
          sprintf(filename2, "obs_i_%d_%d_%d.tab", ntask, track, xtrack);
          write_obs(NULL, filename2, &ctl, &obs_i);
        }
      }
    }
 
    /* ------------------------------------------------------------
       Finalize...
       ------------------------------------------------------------ */
 
    /* Write netCDF file... */
    write_nc(filename, &ncd);
 
    /* Write info... */
    printf("Retrieval finished on rank %d of %d!\n", rank, size);
  }
 
  /* Close file list... */
  fclose(in);
 
  /* Measure CPU time... */
  TIMER("total", 3);
 
  /* Report memory usage... */
  printf("MEMORY_ATM = %g MByte\n", 4. * sizeof(atm_t) / 1024. / 1024.);
  printf("MEMORY_CTL = %g MByte\n", 1. * sizeof(ctl_t) / 1024. / 1024.);
  printf("MEMORY_NCD = %g MByte\n", 1. * sizeof(ncd_t) / 1024. / 1024.);
  printf("MEMORY_OBS = %g MByte\n", 3. * sizeof(atm_t) / 1024. / 1024.);
  printf("MEMORY_RET = %g MByte\n", 1. * sizeof(ret_t) / 1024. / 1024.);
  printf("MEMORY_TBL = %g MByte\n", 1. * sizeof(tbl_t) / 1024. / 1024.);
 
  /* Report problem size... */
  printf("SIZE_TASKS = %d\n", size);
  printf("SIZE_THREADS = %d\n", omp_get_max_threads());
 
  /* MPI... */
  MPI_Finalize();
 
  return EXIT_SUCCESS;
}

Here is the call graph for this function: