5                {
    6 
    8 
    9  char set[
LEN], pertname[
LEN];
 
   10 
   11  double nedt = 0, sza2 = 0;
   12 
   13  int orb = 0;
   14 
   15  FILE *out;
   16 
   17  
   18  if (argc < 4)
   19    ERRMSG(
"Give parameters: <ctl> <pert.nc> <map.tab>");
 
   20 
   21  
   22  scan_ctl(argc, argv, 
"PERTNAME", -1, 
"4mu", pertname);
 
   23  scan_ctl(argc, argv, 
"SET", -1, 
"full", set);
 
   24  const int orbit = (int) 
scan_ctl(argc, argv, 
"ORBIT", -1, 
"-999", NULL);
 
   25  int track0 = (int) 
scan_ctl(argc, argv, 
"TRACK0", -1, 
"0", NULL);
 
   26  int track1 = (int) 
scan_ctl(argc, argv, 
"TRACK1", -1, 
"100000", NULL);
 
   27  int xtrack0 = (int) 
scan_ctl(argc, argv, 
"XTRACK0", -1, 
"0", NULL);
 
   28  int xtrack1 = (int) 
scan_ctl(argc, argv, 
"XTRACK1", -1, 
"29", NULL);
 
   29  int ifov0 = (int) 
scan_ctl(argc, argv, 
"IFOV0", -1, 
"0", NULL);
 
   30  int ifov1 = (int) 
scan_ctl(argc, argv, 
"IFOV1", -1, 
"8", NULL);
 
   31  const double orblat = 
scan_ctl(argc, argv, 
"ORBLAT", -1, 
"0", NULL);
 
   32  const double t0 = 
scan_ctl(argc, argv, 
"T0", -1, 
"-1e100", NULL);
 
   33  const double t1 = 
scan_ctl(argc, argv, 
"T1", -1, 
"1e100", NULL);
 
   34  const double sza0 = 
scan_ctl(argc, argv, 
"SZA0", -1, 
"-1e100", NULL);
 
   35  const double sza1 = 
scan_ctl(argc, argv, 
"SZA1", -1, 
"1e100", NULL);
 
   36  const double dt230 = 
scan_ctl(argc, argv, 
"DT230", -1, 
"-999", NULL);
 
   37  const double nu = 
scan_ctl(argc, argv, 
"NU", -1, 
"-999", NULL);
 
   38  const int dc = (int) 
scan_ctl(argc, argv, 
"DC", -1, 
"0", NULL);
 
   39 
   40  
   42 
   43  
   45 
   46  
   47  track0 = GSL_MIN(GSL_MAX(track0, 0), pert->
ntrack - 1);
 
   48  track1 = GSL_MIN(GSL_MAX(track1, 0), pert->
ntrack - 1);
 
   49  xtrack0 = GSL_MIN(GSL_MAX(xtrack0, 0), pert->
nxtrack - 1);
 
   50  xtrack1 = GSL_MIN(GSL_MAX(xtrack1, 0), pert->
nxtrack - 1);
 
   51  ifov0 = GSL_MIN(GSL_MAX(ifov0, 0), pert->
nfov - 1);
 
   52  ifov1 = GSL_MIN(GSL_MAX(ifov1, 0), pert->
nfov - 1);
 
   53 
   54  
   55  LOG(1, 
"Write perturbation data: %s", argv[3]);
 
   56  if (!(out = fopen(argv[3], "w")))
   57    ERRMSG(
"Cannot create file!");
 
   58 
   59  
   60  fprintf(out,
   61          "# $1 = time (seconds since 01-JAN-2000, 00:00 UTC)\n"
   62          "# $2 = solar zenith angle [deg]\n"
   63          "# $3 = longitude [deg]\n"
   64          "# $4 = latitude [deg]\n"
   65          "# $5 = cloud channel brightness temperature [K]\n"
   66          "# $6 = %s brightness temperature [K]\n"
   67          "# $7 = %s brightness temperature perturbation [K]\n"
   68          "# $8 = %s brightness temperature variance [K^2]\n"
   69          "# $9 = along-track index\n"
   70          "# $10 = across-track index\n"
   71          "# $11 = field-of-view index\n", pertname, pertname, pertname);
   72 
   73  
   74  for (int itrack = track0; itrack <= track1; itrack++) {
   75 
   76    
   77    if (itrack > 0)
   78      if (pert->
lat[itrack - 1][pert->
nxtrack / 2][0] <= orblat
 
   79          && pert->
lat[itrack][pert->
nxtrack / 2][0] >= orblat)
 
   80        orb++;
   81 
   82    
   83    fprintf(out, "\n");
   84 
   85    
   86    if (itrack > 0 && pert->
time[itrack][pert->
nxtrack / 2][0]
 
   87        - pert->
time[itrack - 1][pert->
nxtrack / 2][0] >= 10)
 
   88      fprintf(out, "\n");
   89 
   90    
   91    for (int ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++)
   92      for (int ifov = ifov0; ifov <= ifov1; ifov++) {
   93 
   94        
   95        if (pert->
lon[itrack][ixtrack][ifov] < -180
 
   96            || pert->
lon[itrack][ixtrack][ifov] > 180
 
   97            || pert->
lat[itrack][ixtrack][ifov] < -90
 
   98            || pert->
lat[itrack][ixtrack][ifov] > 90)
 
   99          continue;
  100 
  101        
  102        int asc =
  103          (pert->
lat[itrack > 0 ? itrack : itrack + 1][pert->
nxtrack / 2]
 
  104           [ifov] > pert->
lat[itrack >
 
  105                              0 ? itrack -
  106                              1 : itrack][pert->
nxtrack / 2][ifov]);
 
  107 
  108        
  109        if (sza0 >= -1e10 && sza0 <= 1e10 && sza1 >= -1e10 && sza1 <= 1e10)
  110          sza2 = 
sza(pert->
time[itrack][ixtrack][ifov],
 
  111                     pert->
lon[itrack][ixtrack][ifov],
 
  112                     pert->
lat[itrack][ixtrack][ifov]);
 
  113 
  114        
  115        if (dt230 > 0 && nu > 0) {
  116          const double nesr = 
PLANCK(230.0 + dt230, nu) - 
PLANCK(230.0, nu);
 
  117          const double tbg =
  118            pert->
bt[itrack][ixtrack][ifov] - pert->
pt[itrack][ixtrack][ifov];
 
  120        }
  121 
  122        
  123        if (orbit < 0 || orb == orbit)
  124          if (set[0] == 'f' || (set[0] == 'a' && asc)
  125              || (set[0] == 'd' && !asc))
  126            if (pert->
time[itrack][ixtrack][ifov] >= t0
 
  127                && pert->
time[itrack][ixtrack][ifov] <= t1
 
  128                && sza2 >= sza0 && sza2 <= sza1)
  129              fprintf(out, "%.2f %g %g %g %g %g %g %g %d %d %d\n",
  130                      pert->
time[itrack][ixtrack][ifov],
 
  131                      sza(pert->
time[itrack][ixtrack][ifov],
 
  132                          pert->
lon[itrack][ixtrack][ifov],
 
  133                          pert->
lat[itrack][ixtrack][ifov]),
 
  134                      pert->
lon[itrack][ixtrack][ifov],
 
  135                      pert->
lat[itrack][ixtrack][ifov],
 
  136                      pert->
dc[itrack][ixtrack][ifov],
 
  137                      pert->
bt[itrack][ixtrack][ifov],
 
  138                      pert->
pt[itrack][ixtrack][ifov],
 
  139                      pert->
var[itrack][ixtrack][ifov] - gsl_pow_2(nedt),
 
  140                      itrack, ixtrack, ifov);
  141      }
  142  }
  143 
  144  
  145  fclose(out);
  146 
  147  
  148  free(pert);
  149 
  150  return EXIT_SUCCESS;
  151}
double sza(const double sec, const double lon, const double lat)
Calculate solar zenith angle.
 
double scan_ctl(int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
 
#define LEN
Maximum length of ASCII data lines.
 
#define BRIGHT(rad, nu)
Compute brightness temperature.
 
#define ERRMSG(...)
Print error message and quit program.
 
#define ALLOC(ptr, type, n)
Allocate memory.
 
#define PLANCK(T, nu)
Compute Planck function.
 
#define LOG(level,...)
Print log message.
 
void read_pert(char *filename, char *pertname, int dc, pert_t *pert)
Read radiance perturbation data.
 
double dc[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (cloud channel) [K].
 
double lat[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Latitude [deg].
 
int nfov
Number of field of views.
 
double pt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature perturbation (4 or 15 micron) [K].
 
int ntrack
Number of along-track values.
 
double var[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature variance (4 or 15 micron) [K].
 
double time[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Time (seconds since 2000-01-01T00:00Z).
 
double bt[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Brightness temperature (4.3 or 15 micron) [K].
 
int nxtrack
Number of across-track values.
 
double lon[PERT_NTRACK][PERT_NXTRACK][PERT_NFOV]
Longitude [deg].