OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
anc_acq.c
Go to the documentation of this file.
1  /* ========================================================================= */
2  /* module anc_acq.c - functions to read alternate ancillary data */
3  /* Intended to initially do the MET and OZ with ECMWF data */
4  /* */
5  /* Written By: W. Robinson, SAIC, Aug, 2013 */
6  /* */
7  /* ========================================================================= */
8  #include "l12_proto.h"
9  #include "met_cvt.h"
10  #include "anc_acq.h"
11  #include <gsl/gsl_errno.h>
12  #include <gsl/gsl_interp2d.h>
13  #include <gsl/gsl_spline2d.h>
14  /* status of the 3 ancillary files: MET1, 2, 3
15  or OZONE1, 2, 3 - also set it up so the 1,2,3 mean which files are usable */
16  #define ANC_STAT_1T 1 /* one time used */
17  #define ANC_STAT_2T_END 4 /* 2 different anc times in end of list
18  = files[0] = [1] != [2] */
19  #define ANC_STAT_2T_START 2 /*2 different anc times in start of list
20  = files[0] ! [1] = [2] */
21  #define ANC_STAT_3T 3 /* all 3 files with different times */
22  #define ANC_STAT_CLIM 0 /* a climatology in use */
23  #define NPRM 7
24  #define ANCBAD -999.
25  #define OZ_KG_M2_TO_DU 1. / 2.1415e-5
26  #define USE_PMSL 1 /* choice to use surface pressure (0) or use the MSL
27  (mean sea level) pressure (1) with appropriate
28  adjustment for land height */
29  #define ANC_SRC_TYP_ECMWF 0 /* types of anc file data */
30  #define ANC_SRC_TYP_STD_HDF 1 /* old HDF from NCEP (met) TOMS (oz) */
31  #define ANC_SRC_TYP_OLCI 2 /* OLCI tie point at data points */
32  #define ANC_SRC_TYP_GMAO 3 /* GMAP FP-IT or FP forecast */
33  #define ANC_SRC_TYP_BAD -1 /* a bad or undefined type */
34 
35  enum out_nam {
36  ZW, MW, PR, WV, RH, SFCP, SFCRH, SFCT, ICEFR
37  };
38 
39  enum out_prof {
40  TPROF, RHPROF, HPROF, QPROF, O3PROF
41  };
42 
43  enum out_aerosol {
44  BCEXTTAU,
45  BCSCATAU,
46  DUEXTTAU,
47  DUSCATAU,
48  SSEXTTAU,
49  SSSCATAU,
50  SUEXTTAU,
51  SUSCATAU,
52  OCEXTTAU,
53  OCSCATAU,
54  TOTEXTTAU,
55  TOTSCATAU
56  };
57 
58  struct met_sto_str_d {
59  float s_lon; /* start longitude for a grid */
60  float lon_step; /* longitude incriment */
61  int nlon; /* # longitude points */
62  int nlat; /* # latitude points */
63  float e_lon; /* end longitude for a grid */
64  float s_lat; /* start latitude for a grid */
65  float lat_step; /* latitude incriment */
66  float e_lat; /* end latitude for a grid */
67  double data_time[3]; /* time (in Julian days and fraction)
68  for the data 1, 2, 3 */
69  int anc_f_stat; /* status of the met data */
70  float *data[3]; /* storage for MET1, 2, 3 */
71  };
72 
73  typedef struct met_sto_str_d met_sto_str;
74  static met_sto_str met_sto[NPRM];
75 
76  static int proc_land;
77 
78  int anc_acq_init(instr *input, l1str *l1rec, int32_t *anc_id)
79  /*******************************************************************
80 
81  anc_acq_init
82 
83  purpose: Identify the incoming MET files and if netcdf, set up the
84  ancillary data so it can be accessed by anc_acq_line
85  For now, only ECMWF netcdf files can be processed which contain
86  only 1 time and (at least) the parameters listed in prm_nm
87 
88  Returns type: int - 0 - good -1 any trouble checking input anc files
89 
90  Parameters: (in calling order)
91  Type Name I/O Description
92  ---- ---- --- -----------
93  instr * input I program inputs
94  l1str * l1rec I level 1 record structure
95  int32_t * anc_id O size 2 ID of anc data being
96  handled: for [met, ozone]
97  0 - ECMWF files, 1 - NCEP/TOMS
98  files, -1 - bad
99 
100  Modification history:
101  Programmer Date Description of change
102  ---------- ---- ---------------------
103  W. Robinson 12-Aug-2013 Original development
104  W. Robinson, SAIC 22 Feb 2016 enhance to handle OLCI tie point
105  ancillary data
106 
107  *******************************************************************/ {
108  char *files[3]; /* the 3 MET file names */
109  /* ECMWF parameter names in the 2 arrays below */
110  /* name descrip ends up as value in met_sto
111  sp sfc pressure pressure at surface (not currently used)
112  tcwv precip water precip water
113  msl MSL pressure seal lvl pressure
114  u10 10 m u wind uwind meridional S->N (+)
115  v10 10 m v wind vwind zonal W->E (+)
116  t2m 2m temp - Relative humidity
117  d2m 2m dewpoint temp /
118  tco3 total Ozone ozone
119  */
120  char *prm_nm_met[] = {"sp", "tcwv", "msl", "u10", "v10", "t2m", "d2m"};
121  char *prm_nm_oz[] = {"tco3"};
122  int32_t n_prm_met = 7, n_prm_oz = 1, sto_ix;
123  static char file_olci_str[FILENAME_MAX];
124  char *file_olci = (char *) 0; /* for olci file name, 0 if not olci */
125  /*
126  * determine if rad data ia OLCI and get tie point met file name if so
127  */
128  if (l1rec->l1file->format == FT_OLCI) {
129  strcpy(file_olci_str, l1rec->l1file->name);
130  char* ptr = strrchr(file_olci_str, '/');
131  if(ptr) {
132  *ptr = 0;
133  strcat(file_olci_str, "/tie_meteo.nc");
134  } else {
135  strcpy(file_olci_str, "tie_meteo.nc");
136  }
137  file_olci = file_olci_str;
138  }
139  proc_land = input->proc_land; /* carry to the anc_acq_lin routine */
140  /*
141  * get a general classification of the ancillary data from the 1st file
142  * name for met and ozone
143  */
144  anc_id[0] = anc_acq_ck(input->met1, file_olci);
145  anc_id[1] = anc_acq_ck(input->ozone1, file_olci);
146  if ((anc_id[0] == -1) || (anc_id[1] == -1))
147  return -1;
148  else {
149  if (anc_id[0] == ANC_SRC_TYP_ECMWF) {
150  /* do ecmwf met data full check and init */
151  files[0] = input->met1;
152  files[1] = input->met2;
153  files[2] = input->met3;
154  sto_ix = 0; /* location of met in storage struct */
155 
156  anc_id[0] = anc_acq_ecmwf_init(files, prm_nm_met, n_prm_met, sto_ix);
157  }
158  /* do olci set-up in the read: anc_acq_lin_olci */
159  if (anc_id[1] == ANC_SRC_TYP_ECMWF) {
160  /* do ecmwf ozone data full check and init */
161  files[0] = input->ozone1;
162  files[1] = input->ozone2;
163  files[2] = input->ozone3;
164  sto_ix = 6; /* location of oz in storage struct */
165 
166  anc_id[1] = anc_acq_ecmwf_init(files, prm_nm_oz, n_prm_oz, sto_ix);
167  }
168  }
169  if ((anc_id[0] == -1) || (anc_id[1] == -1))
170  return -1;
171  else
172  return 0;
173  }
174 
175  int32_t anc_acq_ck(char *file, char *file_olci)
176  /*******************************************************************
177 
178  anc_acq_ck
179 
180  purpose: Identify the incoming MET or OZONE file #1 as either netcdf OLCI,
181  netcdf ECMWF, or not netcdf, which means it is a std older file format
182 
183  Returns type: int ancillary data file source type, see file top for defs:
184  form: ANC_SRC_TYP_<type>
185 
186  Parameters: (in calling order)
187  Type Name I/O Description
188  ---- ---- --- -----------
189  char * file I 1st anc file name
190  char * file_olci I if null, data type is not OLCI
191  else, the name of the olci met
192  tie point data file
193 
194  Modification history:
195  Programmer Date Description of change
196  ---------- ---- ---------------------
197  W. Robinson 22-Feb-2016 Original development
198 
199  *******************************************************************/ {
200  size_t attlen = 0;
201  int status, ncid, fstat, tie_pt;
202  char *str_attr, lcl_fil[FILENAME_MAX];
203  /*
204  * make sure 1st file has something
205  */
206  fstat = ANC_SRC_TYP_BAD;
207  if ((file == NULL) || (file[0] == 0)) {
208  fprintf(stderr, "-E- %s %d: First anc file name is null\n",
209  __FILE__, __LINE__);
210  return -1;
211  }
212  /*
213  * First, any keyword substitutions
214  * If OLCI data and file name is special OLCI designation: OLCI_TIE_METEO,
215  * set the file to the olci met tie point file
216  */
217  strcpy(lcl_fil, file);
218  if ((file_olci != (char *) 0) &&
219  (strcmp(upcase(lcl_fil), "OLCI_TIE_METEO") == 0)) {
220  fprintf(stderr,
221  "-I- %s %d: Ancillary file will be taken from OLCI metadata tie point meteo file: %s\n",
222  __FILE__, __LINE__, file_olci);
223  strcpy(file, file_olci);
224  }
225  /*
226  * next, see if file opens with netcdf
227  */
228  // first check if it is an HDF4 file
229  if(Hishdf(file)) {
230  return ANC_SRC_TYP_STD_HDF;
231  }
232 
233  if (nc_open(file, 0, &ncid) != NC_NOERR) {
234  //File is not netcdf, thus assume it's a "standard" (HDF4) file
235  return ANC_SRC_TYP_STD_HDF;
236  } else {
237  /*
238  * for the OLCI data only, see if the file is tie point data
239  */
240  tie_pt = 0;
241  if (file_olci != (char *) 0) /* data type is OLCI */ {
242  status = nc_inq_attlen(ncid, NC_GLOBAL, "title", &attlen);
243  if (status == NC_NOERR) {
244  str_attr = (char *) calloc(attlen + 1, sizeof (char));
245  status = nc_get_att_text(ncid, NC_GLOBAL, "title", str_attr);
246  if (status != NC_NOERR) {
247  fprintf(stderr, "-I- %s %d: nc_get_att_string error, file: %s\n",
248  __FILE__, __LINE__, file);
249  return ANC_SRC_TYP_BAD;
250  }
251 
252  if (strcmp(str_attr, "OLCI Level 1b Product, Tie-Point Meteo Data Set")
253  == 0) {
254  free(str_attr);
255  /* should free the space with nc_free_string() as said in nc_get_att_string call */
256  fstat = ANC_SRC_TYP_OLCI;
257  tie_pt = 1;
258  } else
259  fstat = ANC_SRC_TYP_BAD;
260  }
261  }
262  /*
263  * look to see if the file is a GMAO file
264  */
265  if (tie_pt == 0) {
266  status = nc_inq_attlen(ncid, NC_GLOBAL, "title", &attlen);
267  if (status == NC_NOERR) {
268  /* Check for any of the GMAO ids set up */
269  str_attr = (char *) calloc(attlen + 1, sizeof (char));
270  status = nc_get_att_text(ncid, NC_GLOBAL, "title", str_attr);
271  if (status != NC_NOERR) {
272  fprintf(stderr, "-I- %s %d: nc_get_att_string error, file: %s\n",
273  __FILE__, __LINE__, file);
274  return ANC_SRC_TYP_BAD;
275  }
276  if (strstr(str_attr, "GMAO") != 0) {
277  fstat = ANC_SRC_TYP_GMAO;
278  }
279  free(str_attr);
280  }
281  }
282  status = nc_close(ncid);
283  return fstat;
284  }
285  /*
286  * data type id not olci or GMAO, but file is netcdf, so only
287  * thing left is ecmwf
288  */
289  return ANC_SRC_TYP_ECMWF;
290 }
291 
292 int32_t anc_acq_f_stat(char **files, char prioritize_files, int32_t n_anc)
293 /*******************************************************************
294  anc_acq_f_stat
295 
296  purpose: find the status of the input set of ancillary files. Enlarged
297  from the treatment for the MET, OZONE 3-file set and added treatment
298  for 1 and 2 files
299 
300  Returns type: int32_t - the status of the files sent in: ANC_STAT_...
301  (see declarations at anc_acq start)
302  -1 for bad
303 
304  Parameters: (in calling order)
305  Type Name I/O Description
306  ---- ---- --- -----------
307  char ** files I list of files of the ancillary
308  data
309  char prioritize_files I 0 - do nothing, 1 - for 3 files
310  with only 2 that are different,
311  make sure it is always the 1st
312  and 2nd files
313  int32_t n_anc I Number of files to consider
314 
315  Modification history:
316  Programmer Date Description of change
317  ---------- ---- ---------------------
318  W. Robinson 2 May 2018 original development
319 
320  *******************************************************************/ {
321  int32_t f12_mtch, f23_mtch, anc_f_stat;
322 
323  anc_f_stat = -1;
324  if (n_anc == 3) {
325  if ((files[1] == NULL) || (files[2] == NULL) || (files[1][0] == 0)
326  || (files[2][0] == 0)) anc_f_stat = ANC_STAT_CLIM;
327  else {
328  f12_mtch = 0;
329  if (strcmp(files[0], files[1]) == 0) f12_mtch = 1;
330 
331  f23_mtch = 0;
332  if (strcmp(files[1], files[2]) == 0) f23_mtch = 1;
333 
334  if ((strcmp(files[0], files[2]) == 0) && (f12_mtch == 0)) {
335  printf("%s, %d E: ANC1 and 3 match while ANC2 different\n",
336  __FILE__, __LINE__);
337  return -1;
338  }
339  if ((f12_mtch == 1) && (f23_mtch == 1)) anc_f_stat = ANC_STAT_1T;
340  else if ((f12_mtch == 1) && (f23_mtch == 0)) anc_f_stat = ANC_STAT_2T_END;
341  else if ((f12_mtch == 0) && (f23_mtch == 1))
342  anc_f_stat = ANC_STAT_2T_START;
343  else
344  anc_f_stat = ANC_STAT_3T;
345  }
346  /* get the 2 different file names in 1st, 2nd slots */
347  if (prioritize_files && (anc_f_stat == ANC_STAT_2T_END)) {
348  anc_f_stat = ANC_STAT_2T_START;
349  files[1] = files[2];
350  }
351  } else if (n_anc == 2) {
352  if ((files[0] == files[1]) || files[1][0] == 0)
353  anc_f_stat = ANC_STAT_1T;
354  else
355  anc_f_stat = ANC_STAT_2T_START;
356  } else {
357  anc_f_stat = ANC_STAT_1T;
358  }
359  /* END of file set identification */
360  if (prioritize_files) {
361  switch (anc_f_stat) {
362  case ANC_STAT_1T:
363  return 1;
364  break;
365  case ANC_STAT_2T_START:
366  case ANC_STAT_2T_END:
367  return 2;
368  break;
369  case ANC_STAT_3T:
370  return 3;
371  break;
372  }
373  }
374  return anc_f_stat;
375 }
376 
377 int32_t anc_acq_lin_met(l1str *l1rec)
378 /*******************************************************************
379 
380  anc_acq_lin_met
381 
382  purpose: process the met data for the l1rec, for now, just do the GMAO
383 
384  Returns type: int32_t - status 0 is good
385 
386  Parameters: (in calling order)
387  Type Name I/O Description
388  ---- ---- --- -----------
389  l1str * l1rec I/O all L1 line info
390 
391  Modification history:
392  Programmer Date Description of change
393  ---------- ---- ---------------------
394  W. Robinson 7 May 2018 Original development
395 
396  *******************************************************************/ {
397  static int32_t firstcall = 0;
398  static int32_t ntim_int; /* # needed times for interp arrays */
399  static float r2d = 57.29577951;
400  static int32_t anc_id = -1;
401  static gen_int_str *met_int, *met_tim;
402 
403  char *files[3]; /* the 3 MET file names */
404  static char file_olci_str[FILENAME_MAX];
405  char *file_olci = (char *) 0; /* for olci file name, 0 if not olci */
406  int32_t n_met = 9; /* # met parms finally needed */
407  int32_t nlvl_int = 1; /* # levels in interpolation 1 in the 2-d case */
408  int32_t ilvl = 0, nlvl = 1; /* for met, only 1 level */
409  int32_t itim, ilon, npix, iprm, t_interp, data_ix[2];
410  float val, wt_t1, uwnd, vwnd, unc, u_unc, v_unc, ws_2;
411  double l_time, anc_times[3], lat, lon, last_time;
412 
413  /* initialize for the met processing */
414  if (firstcall == 0) {
415  firstcall = 1;
416  /*
417  * determine if rad data is OLCI and get tie point met file name if so
418  */
419  if (l1rec->l1file->format == FT_OLCI) {
420  strcpy(file_olci_str, l1rec->l1file->name);
421  char* ptr = strrchr(file_olci_str, '/');
422  if(ptr) {
423  *ptr = 0;
424  strcat(file_olci_str, "/tie_meteo.nc");
425  } else {
426  strcpy(file_olci_str, "tie_meteo.nc");
427  }
428  file_olci = file_olci_str;
429  }
430  /* get the files from l1rec->met1,2,3 */
431  files[0] = input->met1;
432  files[1] = input->met2;
433  files[2] = input->met3;
434 
435  /* although done previously (for now), do again to order the
436  important files */
437  if ((ntim_int = anc_acq_f_stat(files, 1, 3)) == -1)
438  return 1;
439 
440  /* set up the interpolation structure array */
441  /* note met_int( ntim_int, nlvl_int, n_met ) with met running fastest */
442  if ((met_int = (gen_int_str *) malloc(n_met * nlvl_int * ntim_int *
443  sizeof (gen_int_str))) == NULL) {
444  fprintf(stderr,
445  "-E- %s %d: Unable to allocate met interpolation array\n",
446  __FILE__, __LINE__);
447  return 1;
448  }
449  /* Check the type of the data, call anc_acq_ck
450  * (not needed as this is only GMAO
451  */
452  if ((anc_id = anc_acq_ck(files[0], file_olci)) == -1)
453  return 1;
454 
455  /* if GMAO type, set up the GMAO met data */
456  if (anc_id == ANC_SRC_TYP_GMAO) {
457  /* loop thru times / files from 1 - ntim_int */
458  last_time = 0;
459  for (itim = 0; itim < ntim_int; itim++) {
460  met_tim = met_int + itim * n_met;
461  /* call anc_acq_gmao_met_prep() *to get all parms for that file */
462  if (anc_acq_gmao_met_prep(files[itim], met_tim) != 0)
463  return 1;
464  if (last_time > met_tim[0].anc_time) {
465  fprintf(stderr,
466  "-E- %s %d: met file times are out of sequence\n",
467  __FILE__, __LINE__);
468  return 1;
469  }
470  last_time = met_tim[0].anc_time;
471  }
472  }/* END GMAO */
473  else {
474  /* report an error - something is wrong - no other types now */
475  fprintf(stderr,
476  "-E- %s %d: Error reading ancillary file: %s\n",
477  __FILE__, __LINE__, files[0]);
478  return -1;
479  }
480  } /* END INITIALIZE */
481 
482  /*
483  * Start getting the data for the line in l1rec
484  * get the time of the current line
485  */
486  l_time = l1rec->scantime;
487  npix = l1rec->npix;
488 
489  /* go through the parameters and interpolate in space and time */
490  for (iprm = 0; iprm < n_met; iprm++) {
491  /* find the times and weighting to use for this line */
492  for (itim = 0; itim < ntim_int; itim++)
493  anc_times[itim] = met_int[ iprm + n_met * itim ].anc_time;
494 
495  /* get the proper times to use and weights */
496  if (anc_acq_fnd_t_interp(l_time, anc_times, ntim_int, &t_interp,
497  data_ix, &wt_t1) != 0) return 1;
498 
499  /* interpolate for 1 or 2 times */
500  for (ilon = 0; ilon < npix; ilon++) {
501  lon = l1rec->lon[ilon];
502  lat = l1rec->lat[ilon];
503  // gsl gets cranky with bad inputs, so...
504  if (lon < -180.0 || lon > 180.0 ||
505  lat <-90.0 || lat > 90.0 ||
506  isnan(lat) || isnan(lon)) {
507  continue;
508  }
509  if (anc_acq_eval_pt(met_int, iprm, ilvl, lat, lon, t_interp,
510  data_ix, wt_t1, ntim_int, nlvl, n_met, &val, &unc) != 0) {
511  fprintf(stderr,
512  "-E- %s %d: Error interpolating to file: %s\n",
513  __FILE__, __LINE__, files[0]);
514  return -1;
515  }
516  /* fill in the proper ancillary data slot */
517  switch (iprm) {
518  case ZW:
519  l1rec->zw[ilon] = val;
520  l1rec->zw_unc[ilon] = unc;
521  break;
522  case MW: /* the mwind comes second, so then combinations cn be done */
523  l1rec->mw[ilon] = val;
524  l1rec->mw_unc[ilon] = unc;
525  /* compute the wind speed, angle with the 2 components */
526  uwnd = l1rec->zw[ilon];
527  vwnd = l1rec->mw[ilon];
528  u_unc = l1rec->zw_unc[ilon];
529  v_unc = l1rec->mw_unc[ilon];
530  if (input->windspeed != -2000) l1rec->ws[ilon] =
531  sqrt(pow(uwnd, 2.) + pow(vwnd, 2.));
532  if (input->windangle != -2000) l1rec->wd[ilon] = atan2f(-uwnd, -vwnd)
533  * r2d;
534  /* deal with the uncertainties in speed, direction */
535  uwnd = fabsf(uwnd);
536  vwnd = fabsf(vwnd);
537  ws_2 = uwnd * uwnd + vwnd * vwnd;
538  if ((uwnd + vwnd) > 0.05 * (u_unc + v_unc)) {
539  l1rec->ws_unc[ilon] = sqrt((uwnd * uwnd * u_unc * u_unc +
540  vwnd * vwnd * v_unc * v_unc) / ws_2);
541  l1rec->wd_unc[ilon] = sqrt(vwnd * vwnd * u_unc * u_unc +
542  uwnd * uwnd * v_unc * v_unc) / ws_2;
543  if (l1rec->wd_unc[ilon] > M_PI) l1rec->wd_unc[ilon] = M_PI;
544  } else {
545  l1rec->ws_unc[ilon] = sqrt(0.5 * (u_unc * u_unc + v_unc * v_unc));
546  l1rec->wd_unc[ilon] = M_PI;
547  }
548  l1rec->wd_unc[ilon] *= r2d;
549  break;
550  case PR:
551  if (input->pressure != -2000) {
552  if (proc_land && (l1rec->height[ilon] != 0))
553  val *= exp(-l1rec->height[ilon] / 8434.);
554  l1rec->pr[ilon] = val;
555  l1rec->pr_unc[ilon] = unc;
556  }
557  break;
558  case WV: /* to make the kg m^-2 to g cm^-2 */
559  if (input->watervapor != -2000)
560  l1rec->wv[ilon] = val / 10.;
561  l1rec->wv_unc[ilon] = unc / 10.;
562  break;
563  case RH:
564  if (input->relhumid != -2000)
565  l1rec->rh[ilon] = val;
566  l1rec->rh_unc[ilon] = unc;
567  break;
568  case SFCP:
569  l1rec->sfcp[ilon] = val;
570  break;
571  case SFCRH:
572  l1rec->sfcrh[ilon] = val;
573  break;
574  case SFCT:
575  l1rec->sfct[ilon] = val;
576  break;
577  case ICEFR:
578  l1rec->icefr[ilon] = val;
579  if (val > input->ice_threshold)
580  l1rec->ice[ilon] = ON;
581  break;
582  }
583  }
584  }
585  /* and end */
586  return 0;
587 }
588 int32_t anc_acq_lin_aerosol(l1str *l1rec) {
589  static int32_t firstcall = 0;
590  static int32_t ntim_int; /* # needed times for interp arrays */
591  static gen_int_str *aer_int, *aer_tim;
592 
593  char *files[3]; /* the 3 MET aerosol file names */
594  int32_t n_met = 12; /* # met parms finally needed */
595  int32_t itim, ilon, npix, iprm, t_interp, data_ix[2];
596  int32_t ilvl = 0, nlvl = 1; /* only 1 level */
597 
598  float val, wt_t1, unc;
599  double l_time, anc_times[3], lat, lon, last_time;
600  anc_aer_struc *anc_aerosol;
601 
602  npix = l1rec->npix;
603  /* initialize for the met processing */
604  if (firstcall == 0) {
605  firstcall = 1;
606  /* get the files from l1rec->anc_profile1,2,3 */
607  files[0] = input->anc_aerosol1;
608  files[1] = input->anc_aerosol2;
609  files[2] = input->anc_aerosol3;
610 
611  printf("\nOpening ancillary aerosol files.\n");
612  printf(" anc_aerosol1 = %s\n", input->anc_aerosol1);
613  printf(" anc_aerosol2 = %s\n", input->anc_aerosol2);
614  printf(" anc_aerosol3 = %s\n", input->anc_aerosol3);
615  printf("\n");
616 
617  /* although done previously (for now), do again to order the
618  important files */
619  if ((ntim_int = anc_acq_f_stat(files, 1, 3)) == -1)
620  return 1;
621 
622  /* for a climatology (ntim_int = 0) just say that climatology
623  profies not implemented yet */
624  if (ntim_int == ANC_STAT_CLIM) {
625  fprintf(stderr,
626  "-E- %s %d: Ancillary met profile climatology not implemented yet\n",
627  __FILE__, __LINE__);
628  return 1;
629  }
630 
631  /* set up the interpolation structure array */
632  if ((aer_int = (gen_int_str *) malloc(n_met * nlvl * ntim_int *
633  sizeof (gen_int_str))) == NULL) {
634  fprintf(stderr,
635  "-E- %s %d: Unable to allocate profile interpolation array\n",
636  __FILE__, __LINE__);
637  return 1;
638  }
639 
640  /* loop thru times / files from 1 - ntim_int */
641  last_time = 0;
642  for (itim = 0; itim < ntim_int; itim++) {
643  aer_tim = aer_int + itim * n_met;
644  /* call anc_acq_gmao_met_prep() *to get all parms for that file */
645  if (anc_acq_gmao_aer_prep(files[itim], aer_tim) != 0)
646  return 1;
647  if (last_time > aer_tim[0].anc_time) {
648  fprintf(stderr,
649  "-E- %s %d: met file times are out of sequence\n",
650  __FILE__, __LINE__);
651  return 1;
652  }
653  last_time = aer_tim[0].anc_time;
654  }
655  } /* END INITIALIZE */
656 
657  /* set up the storage for the profile data if needed */
658  anc_aerosol = l1rec->anc_aerosol;
659  if (anc_aerosol == NULL) {
660  if (init_anc_aerosol(l1rec) != 0)
661  return 1;
662  }
663  /*
664  * Start getting the data for the line in l1rec
665  * get the time of the current line
666  */
667  l_time = l1rec->scantime;
668  npix = l1rec->npix;
669 
670  /* go through the parameters and interpolate in space and time */
671  for (iprm = 0; iprm < n_met; iprm++) {
672  /* find the times and weighting to use for this line */
673  for (itim = 0; itim < ntim_int; itim++)
674  anc_times[itim] = aer_int[ iprm + n_met * itim ].anc_time;
675 
676  /* get the proper times to use and weights */
677  if (anc_acq_fnd_t_interp(l_time, anc_times, ntim_int, &t_interp,
678  data_ix, &wt_t1) != 0) return 1;
679 
680  /* interpolate for 1 or 2 times */
681  for (ilon = 0; ilon < npix; ilon++) {
682  lon = l1rec->lon[ilon];
683  lat = l1rec->lat[ilon];
684  // gsl gets cranky with bad inputs, so...
685  if (lon < -180.0 || lon > 180.0 ||
686  lat <-90.0 || lat > 90.0 ||
687  isnan(lat) || isnan(lon)){
688  continue;
689  }
690  if (anc_acq_eval_pt(aer_int, iprm, ilvl, lat, lon, t_interp,
691  data_ix, wt_t1, ntim_int, nlvl, n_met, &val, &unc) != 0) {
692  fprintf(stderr,
693  "-E- %s %d: Error interpolating to file: %s\n",
694  __FILE__, __LINE__, files[0]);
695  return -1;
696  }
697  /* fill in the proper ancillary data slot */
698  switch (iprm) {
699  case BCEXTTAU:
700  l1rec->anc_aerosol->black_carbon_ext[ilon] = val;
701  break;
702  case BCSCATAU:
703  l1rec->anc_aerosol->black_carbon_scat[ilon] = val;
704  break;
705  case DUEXTTAU:
706  l1rec->anc_aerosol->dust_ext[ilon] = val;
707  break;
708  case DUSCATAU:
709  l1rec->anc_aerosol->dust_scat[ilon] = val;
710  break;
711  case SSEXTTAU:
712  l1rec->anc_aerosol->sea_salt_ext[ilon] = val;
713  break;
714  case SSSCATAU:
715  l1rec->anc_aerosol->sea_salt_scat[ilon] = val;
716  break;
717  case SUEXTTAU:
718  l1rec->anc_aerosol->sulphur_ext[ilon] = val;
719  break;
720  case SUSCATAU:
721  l1rec->anc_aerosol->sulphur_scat[ilon] = val;
722  break;
723  case OCEXTTAU:
724  l1rec->anc_aerosol->organic_carbon_ext[ilon] = val;
725  break;
726  case OCSCATAU:
727  l1rec->anc_aerosol->organic_carbon_scat[ilon] = val;
728  break;
729  case TOTEXTTAU:
730  l1rec->anc_aerosol->total_aerosol_ext[ilon] = val;
731  break;
732  case TOTSCATAU:
733  l1rec->anc_aerosol->total_aerosol_scat[ilon] = val;
734  break;
735  }
736  }
737  }
738  /* and end */
739  return 0;
740 }
741 
742 int32_t anc_acq_lin_prof(l1str *l1rec)
743 /*******************************************************************
744 
745  anc_acq_lin_prof
746 
747  purpose: process the met profile data for the l1rec, from the GMAO
748  FP-IT data
749 
750  Returns type: int32_t - status 0 is good
751 
752  Parameters: (in calling order)
753  Type Name I/O Description
754  ---- ---- --- -----------
755  l1str * l1rec I/O all L1 line info
756 
757  Modification history:
758  Programmer Date Description of change
759  ---------- ---- ---------------------
760  W. Robinson 7 May 2018 Original development
761 
762  *******************************************************************/ {
763  static int32_t firstcall = 0;
764  static int32_t ntim_int; /* # needed times for interp arrays */
765  static gen_int_str *prof_int, *prof_tim;
766 
767  char *files[3]; /* the 3 MET profile file names */
768  int32_t n_met = 5; /* # met parms finally needed */
769  int32_t nlvl = l1rec->l1file->nlvl; /* # levels in interpolation 42 for
770  GMAO, set in filehandle_init.c */
771  int32_t itim, ilon, ilvl, npix, iprm, loc, t_interp, data_ix[2];
772  float val, wt_t1, unc;
773  double l_time, anc_times[3], lat, lon, last_time;
774  anc_struc *anc_add;
775 
776  npix = l1rec->npix;
777  /* initialize for the met processing */
778  if (firstcall == 0) {
779  firstcall = 1;
780  /* get the files from l1rec->anc_profile1,2,3 */
781  files[0] = input->anc_profile1;
782  files[1] = input->anc_profile2;
783  files[2] = input->anc_profile3;
784 
785  printf("\nOpening meteorological profile files.\n");
786  printf(" anc_profile1 = %s\n", input->anc_profile1);
787  printf(" anc_profile2 = %s\n", input->anc_profile2);
788  printf(" anc_profile3 = %s\n", input->anc_profile3);
789  printf("\n");
790 
791  /* although done previously (for now), do again to order the
792  important files */
793  if ((ntim_int = anc_acq_f_stat(files, 1, 3)) == -1)
794  return 1;
795 
796  /* for a climatology (ntim_int = 0) just say that climatology
797  profies not implemented yet */
798  if (ntim_int == ANC_STAT_CLIM) {
799  fprintf(stderr,
800  "-E- %s %d: Ancillary met profile climatology not implemented yet\n",
801  __FILE__, __LINE__);
802  return 1;
803  }
804 
805  /* set up the interpolation structure array */
806  /* note prof_int( ntim_int, nlvl, n_met ) with met running fastest */
807  if ((prof_int = (gen_int_str *) malloc(n_met * nlvl * ntim_int *
808  sizeof (gen_int_str))) == NULL) {
809  fprintf(stderr,
810  "-E- %s %d: Unable to allocate profile interpolation array\n",
811  __FILE__, __LINE__);
812  return 1;
813  }
814 
815  /* get GMAO profile set-up */
816  /* loop thru times / files from 1 - ntim_int */
817  last_time = 0;
818  for (itim = 0; itim < ntim_int; itim++) {
819  prof_tim = prof_int + n_met * nlvl * itim;
820  /* call anc_acq_gmao_prof_prep() *to get all parms for that file */
821  if (anc_acq_gmao_prof_prep(files[itim], prof_tim, nlvl) != 0)
822  return 1;
823  if (last_time > prof_tim[0].anc_time) {
824  fprintf(stderr,
825  "-E- %s %d: met file times are out of sequence\n",
826  __FILE__, __LINE__);
827  return 1;
828  }
829  last_time = prof_tim[0].anc_time;
830  }
831  } /* END INITIALIZE */
832 
833  /* set up the storage for the profile data if needed */
834  anc_add = l1rec->anc_add;
835  if (anc_add == NULL) {
836  if (init_anc_add(l1rec) != 0)
837  return 1;
838  }
839 
840  /*
841  * Start getting the data for the line in l1rec
842  * get the time of the current line
843  */
844  l_time = l1rec->scantime;
845 
846  /* go through the parameters, levels and interpolate in space and time */
847  for (iprm = 0; iprm < n_met; iprm++) {
848  for (ilvl = 0; ilvl < nlvl; ilvl++) {
849  /* find the times and weighting to use for this line */
850  for (itim = 0; itim < ntim_int; itim++) {
851  loc = iprm + n_met * (ilvl + nlvl * itim);
852  anc_times[itim] = prof_int[ loc ].anc_time;
853  }
854 
855  /* get the proper times to use and weights */
856  if (anc_acq_fnd_t_interp(l_time, anc_times, ntim_int, &t_interp,
857  data_ix, &wt_t1) != 0) return 1;
858 
859  /* interpolate for 1 or 2 times */
860  for (ilon = 0; ilon < npix; ilon++) {
861  lon = l1rec->lon[ilon];
862  lat = l1rec->lat[ilon];
863  // gsl gets cranky with bad inputs, so...
864  if (lon < -180.0 || lon > 180.0 ||
865  lat <-90.0 || lat > 90.0 ||
866  isnan(lat) || isnan(lon)){
867  continue;
868  }
869  if (anc_acq_eval_pt(prof_int, iprm, ilvl, lat, lon, t_interp,
870  data_ix, wt_t1, ntim_int, nlvl, n_met, &val, &unc) != 0) {
871  fprintf(stderr,
872  "-E- %s %d: Error interpolating file: %s\n",
873  __FILE__, __LINE__, files[0]);
874  return -1;
875  }
876  /* fill in the proper ancillary data slot */
877  loc = ilvl + nlvl * ilon;
878  switch (iprm) {
879  case TPROF:
880  l1rec->anc_add->prof_temp[loc] = val;
881  break;
882  case RHPROF:
883  l1rec->anc_add->prof_rh[loc] = val;
884  break;
885  case HPROF:
886  l1rec->anc_add->prof_height[loc] = val;
887  break;
888  case QPROF:
889  l1rec->anc_add->prof_q[loc] = val;
890  break;
891  case O3PROF:
892  l1rec->anc_add->prof_o3[loc] = val;
893  break;
894  }
895  }
896  }
897  }
898  /* and end */
899  return 0;
900 }
901 
902 int32_t init_anc_aerosol(l1str *l1rec) {
903 
904  if ((l1rec->anc_aerosol =
905  (anc_aer_struc *) malloc(sizeof ( anc_aer_struc))) == NULL) {
906  fprintf(stderr,
907  "-E- %s %d: Unable to allocate additional ancillary data storage 1\n",
908  __FILE__, __LINE__);
909  return 1;
910  }
911  if (((l1rec->anc_aerosol->black_carbon_ext =
912  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
913  ((l1rec->anc_aerosol->black_carbon_scat =
914  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
915  ((l1rec->anc_aerosol->dust_ext =
916  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
917  ((l1rec->anc_aerosol->dust_scat =
918  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
919  ((l1rec->anc_aerosol->sea_salt_ext =
920  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
921  ((l1rec->anc_aerosol->sea_salt_scat =
922  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
923  ((l1rec->anc_aerosol->sulphur_ext =
924  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
925  ((l1rec->anc_aerosol->sulphur_scat =
926  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
927  ((l1rec->anc_aerosol->organic_carbon_ext =
928  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
929  ((l1rec->anc_aerosol->organic_carbon_scat =
930  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
931  ((l1rec->anc_aerosol->total_aerosol_ext =
932  (float *) malloc(l1rec->npix * sizeof (float))) == NULL) ||
933  ((l1rec->anc_aerosol->total_aerosol_scat =
934  (float *) malloc(l1rec->npix * sizeof (float))) == NULL)) {
935  fprintf(stderr,
936  "-E- %s %d: Unable to allocate additional ancillary data storage 2\n",
937  __FILE__, __LINE__);
938  return 1;
939  }
940  return 0;
941 }
942 
943 int init_anc_add(l1str *l1rec)
944 /*-----------------------------------------------------------------------
945  init_anc_add
946 
947  purpose: general storage allocation for added ancillary parameters
948 
949  Returns 0 if all checks are OK
950  Parameters: (in calling order)
951  Type Name I/O Description
952  ---- ---- --- -----------
953  l1str * l1rec I record containing the band-
954  dependent geometry fields and
955  the sizes in pixels, bands
956 
957  Modification history:
958  Programmer Date Description of change
959  ---------- ---- ---------------------
960  Wayne Robinson 22 June 2018 Original development
961 
962  -----------------------------------------------------------------------*/ {
963  int32_t npix, nlvl = 42;
964 
965  npix = l1rec->npix;
966 
967  if ((l1rec->anc_add =
968  (anc_struc *) malloc(sizeof ( anc_struc))) == NULL) {
969  fprintf(stderr,
970  "-E- %s %d: Unable to allocate additional ancillary data storage 1\n",
971  __FILE__, __LINE__);
972  return 1;
973  }
974  if (((l1rec->anc_add->prof_temp =
975  (float *) malloc(nlvl * npix * sizeof (float))) == NULL) ||
976  ((l1rec->anc_add->prof_rh =
977  (float *) malloc(nlvl * npix * sizeof (float))) == NULL) ||
978  ((l1rec->anc_add->prof_height =
979  (float *) malloc(nlvl * npix * sizeof (float))) == NULL) ||
980  ((l1rec->anc_add->prof_q =
981  (float *) malloc(nlvl * npix * sizeof (float))) == NULL) ||
982  ((l1rec->anc_add->prof_o3 =
983  (float *) malloc(nlvl * npix * sizeof (float))) == NULL)) {
984  fprintf(stderr,
985  "-E- %s %d: Unable to allocate additional ancillary data storage 2\n",
986  __FILE__, __LINE__);
987  return 1;
988  }
989  l1rec->anc_add->nlvl = nlvl;
990  return 0;
991 }
992 
993 int32_t anc_acq_lin_oz(l1str *l1rec)
994 /*******************************************************************
995 
996  anc_acq_lin_met
997 
998  purpose: process the ozone data for the l1rec, for now, just do the GMAO
999 
1000  Returns type: int32_t - status 0 is good
1001 
1002  Parameters: (in calling order)
1003  Type Name I/O Description
1004  ---- ---- --- -----------
1005  l1str * l1rec I/O all L1 line info
1006 
1007  Modification history:
1008  Programmer Date Description of change
1009  ---------- ---- ---------------------
1010  W. Robinson 12 June 2018 Original development
1011 
1012  *******************************************************************/ {
1013  static int32_t firstcall = 0;
1014  static int32_t ntim_int; /* # needed times for interp arrays */
1015  static int32_t anc_id = -1;
1016  static gen_int_str *oz_int, *oz_tim;
1017 
1018  char *files[3]; /* the 3 OZONE file names */
1019  static char file_olci_str[FILENAME_MAX];
1020  char *file_olci = (char *) 0; /* for olci file name, 0 if not olci */
1021  int32_t n_oz = 1; /* # met parms finally needed */
1022  int32_t nlvl_int = 1; /* # levels in interpolation 1 in the 2-d case */
1023  int32_t ilvl = 0, nlvl = 1; /* for oz, only 1 level */
1024  int32_t itim, ilon, npix, iprm, t_interp, data_ix[2];
1025  float val, wt_t1, unc;
1026  double l_time, anc_times[3], lat, lon, last_time;
1027 
1028  /* initialize for the ozone processing */
1029  if (firstcall == 0) {
1030  firstcall = 1;
1031  /*
1032  * determine if rad data is OLCI and get tie point met file name if so
1033  */
1034  if (l1rec->l1file->format == FT_OLCI) {
1035  strcpy(file_olci_str, l1rec->l1file->name);
1036  char* ptr = strrchr(file_olci_str, '/');
1037  if(ptr) {
1038  *ptr = 0;
1039  strcat(file_olci_str, "/tie_meteo.nc");
1040  } else {
1041  strcpy(file_olci_str, "tie_meteo.nc");
1042  }
1043  file_olci = file_olci_str;
1044  }
1045  /* get the files from l1rec->met1,2,3 */
1046  files[0] = input->ozone1;
1047  files[1] = input->ozone2;
1048  files[2] = input->ozone3;
1049 
1050  /* although done previously (for now), do again to order the
1051  important files */
1052  if ((ntim_int = anc_acq_f_stat(files, 1, 3)) == -1)
1053  return 1;
1054 
1055  /* set up the interpolation structure array */
1056  /* note oz_int( ntim_int, nlvl_int, n_oz ) with oz running fastest */
1057  if ((oz_int = (gen_int_str *) malloc(n_oz * nlvl_int * ntim_int *
1058  sizeof (gen_int_str))) == NULL) {
1059  fprintf(stderr,
1060  "-E- %s %d: Unable to allocate met interpolation array\n",
1061  __FILE__, __LINE__);
1062  return 1;
1063  }
1064  /* Check the type of the data, call anc_acq_ck
1065  * (not needed as this is only GMAO
1066  */
1067  if ((anc_id = anc_acq_ck(files[0], file_olci)) == -1)
1068  return 1;
1069 
1070  /* if GMAO type, set up the GMAO met data */
1071  if (anc_id == ANC_SRC_TYP_GMAO) {
1072  /* loop thru times / files from 1 - ntim_int */
1073  last_time = 0;
1074  for (itim = 0; itim < ntim_int; itim++) {
1075  oz_tim = oz_int + itim * n_oz;
1076  /* call anc_acq_gmao_oz_prep() *to get all parms for that file */
1077  if (anc_acq_gmao_oz_prep(files[itim], oz_tim) != 0)
1078  return 1;
1079  if (last_time > oz_tim[0].anc_time) {
1080  fprintf(stderr,
1081  "-E- %s %d: ozone file times are out of sequence\n",
1082  __FILE__, __LINE__);
1083  return 1;
1084  }
1085  last_time = oz_tim[0].anc_time;
1086  }
1087  }/* END GMAO */
1088  else {
1089  /* report an error - something is wrong - no other types now */
1090  fprintf(stderr,
1091  "-E- %s %d: Error reading file: %s\n",
1092  __FILE__, __LINE__, files[0]);
1093  return -1;
1094  }
1095  } /* END INITIALIZE */
1096 
1097  /*
1098  * Start getting the data for the line in l1rec
1099  * get the time of the current line
1100  */
1101  l_time = l1rec->scantime;
1102  npix = l1rec->npix;
1103 
1104  /* go through the parameters and interpolate in space and time */
1105  for (iprm = 0; iprm < n_oz; iprm++) {
1106  /* find the times and weighting to use for this line */
1107  for (itim = 0; itim < ntim_int; itim++)
1108  anc_times[itim] = oz_int[ iprm + n_oz * itim ].anc_time;
1109 
1110  /* get the proper times to use and weights */
1111  if (anc_acq_fnd_t_interp(l_time, anc_times, ntim_int, &t_interp,
1112  data_ix, &wt_t1) != 0) return 1;
1113 
1114  /* interpolate for 1 or 2 times */
1115  if (input->ozone != -2000) {
1116  for (ilon = 0; ilon < npix; ilon++) {
1117  lon = l1rec->lon[ilon];
1118  lat = l1rec->lat[ilon];
1119  // gsl gets cranky with bad inputs, so...
1120  if (lon < -180.0 || lon > 180.0 ||
1121  lat <-90.0 || lat > 90.0 ||
1122  isnan(lat) || isnan(lon)){
1123  continue;
1124  }
1125  if (anc_acq_eval_pt(oz_int, iprm, ilvl, lat, lon, t_interp,
1126  data_ix, wt_t1, ntim_int, nlvl, n_oz, &val, &unc) != 0) {
1127  fprintf(stderr,
1128  "-E- %s %d: Error interpolating file: %s\n",
1129  __FILE__, __LINE__, files[0]);
1130  return -1;
1131  }
1132  /* fill in the proper ozone data slot converted from DU to
1133  cm thickness*/
1134  l1rec->oz[ilon] = val / 1000.;
1135  l1rec->oz_unc[ilon] = unc / 1000.;
1136  }
1137  }
1138  }
1139  /* and end */
1140  return 0;
1141 }
1142 
1143 int32_t anc_acq_gmao_met_prep(char *file, gen_int_str *met_int)
1144 /*******************************************************************
1145 
1146  anc_acq_gmao_met_prep
1147 
1148  purpose: set up the interpolation objects and qc arrays for all the
1149  GMAO met parms for a file
1150 
1151  Returns type: int32_t - status 0 is good
1152 
1153  Parameters: (in calling order)
1154  Type Name I/O Description
1155  ---- ---- --- -----------
1156  char * file I GMAO single-level file
1157  gen_int_str * met_int I/O array of interpolation
1158  structures for the met data
1159 
1160  Modification history:
1161  Programmer Date Description of change
1162  ---------- ---- ---------------------
1163  W. Robinson 7 May 2018 Original development
1164  W. Robinson 16 Jan 2019 adapt to read the T10M product for
1165  sfc T
1166 
1167  *******************************************************************/ {
1168  /* list of GMAO groups, parm names and special processing value */
1169  /* note that for the second RH, use the the 1st RH */
1170  char *ob_gmao_grp[] = {"met", "met", "met", "met", "met", "met", "met", "met", "ocn_ice"};
1171  char *ob_gmao_prm_nm[] = {"U10M", "V10M", "SLP", "TQV", "PS", "PS", "PS", "T10M", "FRSEAICE"};
1172  /* also will need for rh: met, T10M and met, QV10M */
1173  /* also will need for icefr: lnd_ice, FRSNO (future) */
1174  int32_t n_raw_gmao = 9;
1175  int32_t iprm, nlat, nlon, nlvl, iv, nv;
1176  float *data2, *data3, *data_rh, p_lcl;
1177  unsigned char *qual, *qual2;
1178  double time, *lat_coord, *lon_coord, *ddata;
1179  double *lat_coord2, *lon_coord2;
1180  static float *data = 0;
1181 
1182  /* out prm GMAO grp(s) GMAO name(s) do
1183  zwind (10 m) zw met U10M as-is
1184  mwind (10 m) mw met V10M as-is
1185  pressure (at MSL) pr met SLP cvt Pa ->HPa
1186  water_vapor (PW) wv met TQV as-is
1187  humidity (1000 mb) rh met PS convert to RH with:
1188  met T10M
1189  met QV10M
1190  sfc pressure sfcp met PS cvt Pa ->HPa
1191  RH sfc sfcrh take above rh for now
1192  T surface sfct met T10M as-is
1193  ice fraction icefr ocn_ice FRSEAICE as-is
1194  lnd_ice FRSNO Future combination with ocean
1195  */
1196 
1197  /* loop for output parms and get the data array(s) and do any special
1198  processing to make final param in final units */
1199  for (iprm = 0; iprm < n_raw_gmao; iprm++) {
1200  /* get the GMAO array */
1201  if (anc_acq_read_gmao(file, ob_gmao_grp[iprm], ob_gmao_prm_nm[iprm],
1202  &data, &qual, &time, &nlon, &nlat, &nlvl, &lon_coord,
1203  &lat_coord) != 0)
1204  return 1;
1205 
1206  nv = nlon * nlat;
1207  if ((ddata = (double *) malloc(nv * sizeof ( double))) == NULL) {
1208  fprintf(stderr, "-E- %s %d: Unable to allocate array ddata\n",
1209  __FILE__, __LINE__);
1210  return 1;
1211  }
1212  /* convert parm to needed type and combine parms to get field needed if
1213  necessary */
1214  switch (iprm) {
1215  case ZW:
1216  case MW:
1217  case WV:
1218  break;
1219  case PR:
1220  for (iv = 0; iv < nv; iv++) {
1221  if (*(qual + iv) == 0) {
1222  p_lcl = *(data + iv);
1223  p_lcl = met_cvt_p_cvt(p_lcl, MET_UNITS__P_PA, MET_UNITS__P_HPA);
1224  *(data + iv) = p_lcl;
1225  }
1226  }
1227  break;
1228  case RH:
1229  if (anc_acq_read_gmao(file, "met", "T10M", &data2, &qual2, &time,
1230  &nlon, &nlat, &nlvl, &lon_coord2, &lat_coord2) != 0)
1231  return 1;
1232  free(lat_coord2);
1233  free(lon_coord2);
1234  for (iv = 0; iv < nv; iv++)
1235  *(qual + iv) = *(qual + iv) | *(qual2 + iv);
1236  free(qual2);
1237 
1238  if (anc_acq_read_gmao(file, "met", "QV10M", &data3, &qual2, &time,
1239  &nlon, &nlat, &nlvl, &lon_coord2, &lat_coord2) != 0)
1240  return 1;
1241  free(lat_coord2);
1242  free(lon_coord2);
1243  for (iv = 0; iv < nv; iv++)
1244  *(qual + iv) = *(qual + iv) | *(qual2 + iv);
1245  free(qual2);
1246 
1247  if ((data_rh = (float *) malloc(nlon * nlat * sizeof ( float)))
1248  == NULL) {
1249  fprintf(stderr, "-E- %s, %d: malloc failure\n", __FILE__, __LINE__);
1250  return 1;
1251  }
1252  met_cvt_q_to_rh(nv, data, MET_UNITS__P_PA, data2, MET_UNITS__T_K,
1253  data3, MET_UNITS__Q_KG_KG, data_rh);
1254 
1255  free(data2);
1256  free(data3);
1257  for (iv = 0; iv < nv; iv++)
1258  *(data + iv) = *(data_rh + iv);
1259  free(data_rh);
1260  break;
1261  case SFCP:
1262  for (iv = 0; iv < nv; iv++) {
1263  if (*(qual + iv) == 0) {
1264  p_lcl = *(data + iv);
1265  p_lcl = met_cvt_p_cvt(p_lcl, MET_UNITS__P_PA, MET_UNITS__P_HPA);
1266  *(data + iv) = p_lcl;
1267  }
1268  }
1269  break;
1270  case SFCRH:
1271  if (anc_acq_read_gmao(file, "met", "T10M", &data2, &qual2, &time,
1272  &nlon, &nlat, &nlvl, &lon_coord2, &lat_coord2) != 0)
1273  return 1;
1274  free(lat_coord2);
1275  free(lon_coord2);
1276  for (iv = 0; iv < nv; iv++)
1277  *(qual + iv) = *(qual + iv) | *(qual2 + iv);
1278  free(qual2);
1279 
1280  if (anc_acq_read_gmao(file, "met", "QV10M", &data3, &qual2, &time,
1281  &nlon, &nlat, &nlvl, &lon_coord2, &lat_coord2) != 0)
1282  return 1;
1283  free(lat_coord2);
1284  free(lon_coord2);
1285  for (iv = 0; iv < nv; iv++)
1286  *(qual + iv) = *(qual + iv) | *(qual2 + iv);
1287  free(qual2);
1288 
1289  if ((data_rh = (float *) malloc(nlon * nlat * sizeof ( float)))
1290  == NULL) {
1291  fprintf(stderr, "-E- %s, %d: malloc failure\n", __FILE__, __LINE__);
1292  return 1;
1293  }
1294  met_cvt_q_to_rh(nv, data, MET_UNITS__P_PA, data2, MET_UNITS__T_K,
1295  data3, MET_UNITS__Q_KG_KG, data_rh);
1296 
1297  free(data2);
1298  free(data3);
1299  for (iv = 0; iv < nv; iv++)
1300  *(data + iv) = *(data_rh + iv);
1301  free(data_rh);
1302  break;
1303  case SFCT:
1304  for (iv = 0; iv < nv; iv++) {
1305  if (*(qual + iv) == 0) {
1306  p_lcl = *(data + iv);
1307  p_lcl = met_cvt_t_cvt(p_lcl, MET_UNITS__T_K, MET_UNITS__T_C);
1308  *(data + iv) = p_lcl;
1309  }
1310  }
1311  break;
1312  case ICEFR:
1313  /* we may want to add the lnd_ice in with the ocn_ice */
1314 
1315  if (anc_acq_read_gmao(file, "lnd_ice", "FRSNO", &data2, &qual2,
1316  &time, &nlon, &nlat, &nlvl, &lon_coord2, &lat_coord2) != 0)
1317  return 1;
1318  free(lat_coord2);
1319  free(lon_coord2);
1320  for (iv = 0; iv < nv; iv++)
1321  *(qual + iv) = *(qual + iv) | *(qual2 + iv);
1322 
1323  /* over land (good qual), choose the snow fraction values */
1324  for (iv = 0; iv < nv; iv++) {
1325  if (*(qual2 + iv) == 0)
1326  *(data + iv) = *(data2 + iv);
1327  }
1328  free(qual2);
1329  free(data2);
1330  break;
1331  default:
1332  fprintf(stderr, "-E- %s %d: Unknown output identifier: %d\n",
1333  __FILE__, __LINE__, iprm);
1334  }
1335  /* make the grid into an interpolation object
1336  have grid: data lat_coord, lon_coord as scales nlon, nlat as lengths
1337  -> met_int[iprm].int_id */
1338 
1339  met_int[iprm].accel_lat = gsl_interp_accel_alloc();
1340  met_int[iprm].accel_lon = gsl_interp_accel_alloc();
1341  met_int[iprm].int_id = gsl_spline2d_alloc(gsl_interp2d_bilinear, nlon, nlat);
1342  met_int[iprm].lat_coord = lat_coord;
1343  met_int[iprm].lon_coord = lon_coord;
1344 
1345  for (iv = 0; iv < nv; iv++)
1346  ddata[iv] = data[iv];
1347  if (gsl_spline2d_init(met_int[iprm].int_id, met_int[iprm].lon_coord,
1348  met_int[iprm].lat_coord, ddata, nlon, nlat) != 0) {
1349  fprintf(stderr, "-E- %s %d: GSL 2-D initialization failed, file: %s\n",
1350  __FILE__, __LINE__, file);
1351  return 1;
1352  }
1353 
1354  /* delete the parameter data array here
1355  and assign the qc, lon and lat coords to the int struct element
1356  */
1357  met_int[iprm].anc_time = time;
1358  met_int[iprm].qual = qual;
1359  met_int[iprm].nlat = nlat;
1360  met_int[iprm].nlon = nlon;
1361  free(data);
1362  free(ddata);
1363  }
1364  return 0;
1365 }
1366 
1367 int32_t anc_acq_gmao_prof_prep(char *file, gen_int_str *prof_int,
1368  int32_t nlvl_expect)
1369 /*******************************************************************
1370 
1371  anc_acq_gmao_prof_prep
1372 
1373  purpose: set up the interpolation objects and qc arrays for all the
1374  GMAO profile parms for a file (at 1 time)
1375 
1376  Returns type: int32_t - status 0 is good
1377 
1378  Parameters: (in calling order)
1379  Type Name I/O Description
1380  ---- ---- --- -----------
1381  char * file I GMAO single-level file
1382  gen_int_str * prof_int I/O array of interpolation
1383  structures for the met
1384  profile data
1385  int32_t nlvl_expect I expected # profile levels
1386 
1387  Modification history:
1388  Programmer Date Description of change
1389  ---------- ---- ---------------------
1390  W. Robinson 20 June 2018 Original development
1391 
1392  *******************************************************************/ {
1393  /* list of GMAO groups, parm names and special processing value */
1394  /* note that for the second RH, use the the 1st RH */
1395  char *ob_gmao_prm_nm[] = {"T", "RH", "H", "QV", "O3"};
1396  int32_t n_raw_gmao = 5;
1397  int32_t iprm, nlat, nlon, nlvl, ilvl, loc, iv, nv, ntot;
1398  float p_lcl;
1399  unsigned char *qual;
1400  double time, *lat_coord, *lon_coord, *ddata;
1401  static float *data = 0;
1402 
1403  /* out prm GMAO grp(s) GMAO name(s) do
1404  temp profile prof_temp NONE T convert K -> C
1405  RH profile prof_rh NONE RH convert fraction to %
1406  height profile prof_h NONE H as-is
1407  specific humidity profile prof_q NONE QV as-is
1408  ozone profile prof_o3 NONE O3 as-is (unless asked)
1409  */
1410 
1411  /* loop for output parms and get the data array(s) and do any special
1412  processing to make final param in final units */
1413  for (iprm = 0; iprm < n_raw_gmao; iprm++) {
1414  /* get the GMAO array */
1415  if (anc_acq_read_gmao(file, "", ob_gmao_prm_nm[iprm], &data, &qual,
1416  &time, &nlon, &nlat, &nlvl, &lon_coord, &lat_coord) != 0)
1417  return 1;
1418 
1419  /* verify standard # levels */
1420  if (nlvl != nlvl_expect) {
1421  fprintf(stderr,
1422  "-E- %s %d: unexpected # profile levels: %d were read from file: %s\n",
1423  __FILE__, __LINE__, nlvl, file);
1424  return 1;
1425  }
1426  /* convert temperature parm to degrees C */
1427  if (iprm == TPROF) {
1428  ntot = nlon * nlat * nlvl;
1429  for (iv = 0; iv < ntot; iv++) {
1430  if (*(qual + iv) == 0) {
1431  p_lcl = *(data + iv);
1432  p_lcl = met_cvt_t_cvt(p_lcl, MET_UNITS__T_K, MET_UNITS__T_C);
1433  *(data + iv) = p_lcl;
1434  }
1435  }
1436  }/* make the RH in % */
1437  else if (iprm == RHPROF) {
1438  ntot = nlon * nlat * nlvl;
1439  for (iv = 0; iv < ntot; iv++) {
1440  if (*(qual + iv) == 0) {
1441  *(data + iv) *= 100.;
1442  }
1443  }
1444  }
1445  /* loop through each level to make separate interpolation objects */
1446  nv = nlon * nlat;
1447  if ((ddata = (double *) malloc(nv * sizeof ( double))) == NULL) {
1448  fprintf(stderr, "-E- %s %d: Unable to allocate array ddata\n",
1449  __FILE__, __LINE__);
1450  return 1;
1451  }
1452  for (ilvl = 0; ilvl < nlvl; ilvl++) {
1453  loc = iprm + n_raw_gmao * ilvl;
1454  /* make the grid into an interpolation object
1455  have grid: data lat_coord, lon_coord as scales nlon, nlat as lengths
1456  -> prof_int[loc].int_id */
1457 
1458  prof_int[loc].accel_lat = gsl_interp_accel_alloc();
1459  prof_int[loc].accel_lon = gsl_interp_accel_alloc();
1460  prof_int[loc].int_id = gsl_spline2d_alloc(gsl_interp2d_bilinear,
1461  nlon, nlat);
1462  prof_int[loc].lat_coord = lat_coord;
1463  prof_int[loc].lon_coord = lon_coord;
1464 
1465  for (iv = 0; iv < nv; iv++)
1466  ddata[iv] = data[ iv + nv * ilvl ];
1467  if (gsl_spline2d_init(prof_int[loc].int_id, prof_int[loc].lon_coord,
1468  prof_int[loc].lat_coord, ddata, nlon, nlat) != 0) {
1469  fprintf(stderr, "-E- %s %d: GSL 2-D initialization failed, file: %s\n",
1470  __FILE__, __LINE__, file);
1471  return 1;
1472  }
1473  prof_int[loc].anc_time = time;
1474  prof_int[loc].qual = qual + nv * ilvl;
1475  prof_int[loc].nlat = nlat;
1476  prof_int[loc].nlon = nlon;
1477  }
1478 
1479  /* delete the parameter data array here
1480  and assign the qc, lon and lat coords to the int struct element
1481  */
1482  free(data);
1483  free(ddata);
1484  }
1485  return 0;
1486 }
1487 
1488 int32_t anc_acq_gmao_oz_prep(char *file, gen_int_str *oz_int)
1489 /*******************************************************************
1490 
1491  anc_acq_gmao_oz_prep
1492 
1493  purpose: set up the interpolation objects and qc arrays for all the
1494  GMAO ozone parms for a file
1495 
1496  Returns type: int32_t - status 0 is good
1497 
1498  Parameters: (in calling order)
1499  Type Name I/O Description
1500  ---- ---- --- -----------
1501  char * file I GMAO single-level file
1502  gen_int_str * oz_int I/O array of interpolation
1503  structures for the ozone data
1504 
1505  Modification history:
1506  Programmer Date Description of change
1507  ---------- ---- ---------------------
1508  W. Robinson 12 Jun 2018 Original development
1509 
1510  *******************************************************************/ {
1511  /* Only the GMAO TO3 is needed from the 'met' file */
1512  int32_t nlat, nlon, nlvl, iv, nv;
1513  int32_t iprm = 0;
1514  unsigned char *qual;
1515  double time, *lat_coord, *lon_coord, *ddata;
1516  static float *data = 0;
1517 
1518  /* read the ozone */
1519  if (anc_acq_read_gmao(file, "met", "TO3", &data, &qual, &time, &nlon,
1520  &nlat, &nlvl, &lon_coord, &lat_coord) != 0)
1521  return 1;
1522 
1523  nv = nlon * nlat;
1524  if ((ddata = (double *) malloc(nv * sizeof ( double))) == NULL) {
1525  fprintf(stderr, "-E- %s %d: Unable to allocate array ddata\n",
1526  __FILE__, __LINE__);
1527  return 1;
1528  }
1529  /* make the grid into an interpolation object
1530  have grid: data lat_coord, lon_coord as scales nlon, nlat as lengths
1531  -> oz_int[iprm].int_id */
1532 
1533  oz_int[iprm].accel_lat = gsl_interp_accel_alloc();
1534  oz_int[iprm].accel_lon = gsl_interp_accel_alloc();
1535  oz_int[iprm].int_id = gsl_spline2d_alloc(gsl_interp2d_bilinear, nlon, nlat);
1536  oz_int[iprm].lat_coord = lat_coord;
1537  oz_int[iprm].lon_coord = lon_coord;
1538 
1539  for (iv = 0; iv < nv; iv++)
1540  ddata[iv] = data[iv];
1541  if (gsl_spline2d_init(oz_int[iprm].int_id, oz_int[iprm].lon_coord,
1542  oz_int[iprm].lat_coord, ddata, nlon, nlat) != 0) {
1543  fprintf(stderr, "-E- %s %d: GSL 2-D initialization failed, file: %s\n",
1544  __FILE__, __LINE__, file);
1545  return 1;
1546  }
1547  /* delete the parameter data array here
1548  and assign the qc, lon and lat coords to the int struct element
1549  */
1550  oz_int[iprm].anc_time = time;
1551  oz_int[iprm].qual = qual;
1552  oz_int[iprm].nlat = nlat;
1553  oz_int[iprm].nlon = nlon;
1554  free(data);
1555  free(ddata);
1556  /* and end */
1557  return 0;
1558 }
1559 
1560 int32_t anc_acq_gmao_aer_prep(char *file, gen_int_str *aer_int) {
1561  /* Only the GMAO TO3 is needed from the 'met' file */
1562  int32_t nlat, nlon, nlvl, iv, nv;
1563  int32_t iprm = 0;
1564  unsigned char *qual;
1565  double time, *lat_coord, *lon_coord, *ddata;
1566  static float *data = 0;
1567 
1568  /* list of GMAO groups, parm names */
1569  char *ob_gmao_prm_nm[] = {"BCEXTTAU", "BCSCATAU", "DUEXTTAU", "DUSCATAU",
1570  "SSEXTTAU", "SSSCATAU", "SUEXTTAU", "SUSCATAU", "OCEXTTAU", "OCSCATAU",
1571  "TOTEXTTAU", "TOTSCATAU"};
1572 
1573  int32_t n_raw_gmao = 12;
1574  /* read the aerosol parameters */
1575  for (iprm = 0; iprm < n_raw_gmao; iprm++) {
1576  if (anc_acq_read_gmao(file, "", ob_gmao_prm_nm[iprm], &data, &qual,
1577  &time, &nlon, &nlat, &nlvl, &lon_coord, &lat_coord) != 0)
1578  return 1;
1579 
1580  nv = nlon * nlat;
1581  if ((ddata = (double *) malloc(nv * sizeof ( double))) == NULL) {
1582  fprintf(stderr, "-E- %s %d: Unable to allocate array ddata\n",
1583  __FILE__, __LINE__);
1584  return 1;
1585  }
1586  /* make the grid into an interpolation object
1587  have grid: data lat_coord, lon_coord as scales nlon, nlat as lengths
1588  -> aer_int[iprm].int_id */
1589 
1590  aer_int[iprm].accel_lat = gsl_interp_accel_alloc();
1591  aer_int[iprm].accel_lon = gsl_interp_accel_alloc();
1592  aer_int[iprm].int_id = gsl_spline2d_alloc(gsl_interp2d_bilinear, nlon, nlat);
1593  aer_int[iprm].lat_coord = lat_coord;
1594  aer_int[iprm].lon_coord = lon_coord;
1595 
1596  for (iv = 0; iv < nv; iv++)
1597  ddata[iv] = data[iv];
1598  if (gsl_spline2d_init(aer_int[iprm].int_id, aer_int[iprm].lon_coord,
1599  aer_int[iprm].lat_coord, ddata, nlon, nlat) != 0) {
1600  fprintf(stderr, "-E- %s %d: GSL 2-D initialization failed, file: %s\n",
1601  __FILE__, __LINE__, file);
1602  return 1;
1603  }
1604  /* delete the parameter data array here
1605  and assign the qc, lon and lat coords to the int struct element
1606  */
1607  aer_int[iprm].anc_time = time;
1608  aer_int[iprm].qual = qual;
1609  aer_int[iprm].nlat = nlat;
1610  aer_int[iprm].nlon = nlon;
1611  free(data);
1612  free(ddata);
1613  /* and end */
1614  }
1615  return 0;
1616 }
1617 
1618 int32_t anc_acq_eval_pt(gen_int_str *met_int, int32_t iprm, int32_t ilvl,
1619  float lat, float lon, int32_t t_interp, int32_t *data_ix, float wt_t1,
1620  int32_t ntim_int, int32_t nlvl, int32_t nprm, float *final_val, float *unc)
1621 /*******************************************************************
1622 
1623  anc_acq_eval_pt
1624 
1625  purpose: find a parameter for a certain lat, lon by interpolating in
1626  space and time
1627 
1628  Returns type: int32_t - status 0 is good
1629 
1630  Parameters: (in calling order)
1631  Type Name I/O Description
1632  ---- ---- --- -----------
1633  gen_int_str * met_int I array of all interpolation
1634  info for parameters, levels,
1635  times
1636  int32_t iprm I parameter index
1637  int32_t ilvl I level index
1638  float lat I latitude
1639  float lon I longitude
1640  int32_t t_interp I time interpolation indicator:
1641  1 - interpolate, 0 - none
1642  int32_t data_ix I indicies of first and possibly
1643  2nd times to interpolate
1644  float wt_t1 I weight to place on 1st time
1645  int32_t iprm I parameter to work on
1646  int32_t ilvl I level to work on
1647  int32_t ntim_int I total # times available
1648  int32_t nlvl I # levels available
1649  int32_t nprm I # parameters available
1650  float * final_val O final interpolated value
1651  float * unc O Uncertainty in value
1652  (absolute diff in values at
1653  2 times used in interpolation)
1654 
1655  Modification history:
1656  Programmer Date Description of change
1657  ---------- ---- ---------------------
1658  W. Robinson 7 May 2018 Original development
1659 
1660  *******************************************************************/ {
1661  float wt_t2;
1662  gsl_interp_accel *accel_lat, *accel_lon;
1663  int32_t met_ptr, tim_ix, itim, ntim;
1664  double val;
1665  int32_t ilat, ilon, nlon, nlat;
1666 
1667  wt_t2 = 1. - wt_t1;
1668  /* the # times is just the interpolation indicator + 1 */
1669  ntim = t_interp + 1;
1670  *unc = 0.;
1671  for (itim = 0; itim < ntim; itim++) {
1672  tim_ix = data_ix[itim];
1673  met_ptr = iprm + nprm * (ilvl + nlvl * tim_ix);
1674  accel_lon = met_int[met_ptr].accel_lon;
1675  accel_lat = met_int[met_ptr].accel_lat;
1676  nlat = met_int[met_ptr].nlat;
1677  nlon = met_int[met_ptr].nlon;
1678 
1679  if (gsl_spline2d_eval_e(met_int[met_ptr].int_id, lon, lat, accel_lon,
1680  accel_lat, &val) != 0) {
1681  fprintf(stderr, "-E- %s %d: gsl_spline2d_eval_e error: %d, %d, %d\n",
1682  __FILE__, __LINE__, itim, iprm, ilvl);
1683  return -1;
1684  }
1685  /* make sure the value is OK and make bad if not */
1686  ilat = gsl_interp_accel_find(accel_lat,
1687  met_int[met_ptr].lat_coord, nlat, lat);
1688  ilon = gsl_interp_accel_find(accel_lon,
1689  met_int[met_ptr].lon_coord, nlon, lon);
1690 
1691  if ((met_int[met_ptr].qual[ilon + nlon * ilat ] == 1) ||
1692  (met_int[met_ptr].qual[ilon + 1 + nlon * ilat ] == 1) ||
1693  (met_int[met_ptr].qual[ilon + nlon * (ilat + 1) ] == 1) ||
1694  (met_int[met_ptr].qual[ilon + 1 + nlon * (ilat + 1) ] == 1))
1695  val = BAD_FLT;
1696  /* for the case where time interpolation is done, get the other
1697  time's interpolated value */
1698  if (itim == 1) {
1699  if (*final_val == BAD_FLT) {
1700  *final_val = (val == BAD_FLT) ? BAD_FLT : val;
1701  } else if (val == BAD_FLT) {
1702  *final_val = (*final_val == BAD_FLT) ? BAD_FLT : *final_val;
1703  } else {
1704  *unc = fabsf(*final_val - (float)val);
1705  *final_val = *final_val * wt_t1 + val * wt_t2;
1706  }
1707  } else {
1708  *final_val = val;
1709  }
1710  }
1711  if (*final_val == BAD_FLT) *unc = BAD_FLT;
1712  return 0;
1713 }
1714 
1715 int32_t anc_acq_fnd_t_interp(double s_time, double *anc_time,
1716  int32_t anc_f_stat, int32_t *t_interp, int32_t *data_ix, float *wt)
1717 /*******************************************************************
1718 
1719  anc_acq_fnd_t_interp
1720 
1721  purpose: determine the anc times to use and their weight from the anc
1722  times and the scan time
1723 
1724  Returns type: int32_t - status 0 is good
1725 
1726  Parameters: (in calling order)
1727  Type Name I/O Description
1728  ---- ---- --- -----------
1729  double s_time I scan time
1730  double * anc_time I times for ancillary data
1731  int32_t anc_f_stat I ancillary time status
1732  (ANC_STAT...) with 1-3 being
1733  the number of active anc times
1734  and 0, >3 other possible
1735  configurations
1736  int32_t * t_interp O interpolation flag: 0 - none
1737  1 - use the 2 files pointed to
1738  int32_t * data_ix O indecies of 1st, possibly 2nd
1739  anc time to use
1740  float * wt O for t_interp = 1, the weights
1741  of time 1 to use for
1742  interpolation. weight for 2nd
1743  time = 1. - wt
1744  NOTE that there could be a anc_t_warn set (1) if we still want to use the
1745  early definitions for bad anc. I'll comment those cases but not impliment
1746 
1747  Modification history:
1748  Programmer Date Description of change
1749  ---------- ---- ---------------------
1750  W. Robinson 7 May 2018 Original development
1751 
1752  *******************************************************************/ {
1753  int32_t data1_ix=0, data2_ix=0;
1754 
1755 
1756  if (anc_f_stat == ANC_STAT_1T) {
1757  /* one time only available */
1758  *t_interp = 0;
1759  data1_ix = 0;
1760  } else if (anc_f_stat == ANC_STAT_2T_START) {
1761  /* 2 different times in the 1, 2 positions */
1762  if (s_time < anc_time[0]) {
1763  /* scan is earlier than 1st time */
1764  *t_interp = 0;
1765  data1_ix = 0;
1766  /* anc_t_warn = 1 */
1767  } else if (s_time > anc_time[1]) {
1768  /* use anc 2 only */
1769  *t_interp = 0;
1770  data1_ix = 1;
1771  } else {
1772  /* in-between anc 1, 2 use data 0, 1 and time interpolate */
1773  *t_interp = 1;
1774  data1_ix = 0;
1775  data2_ix = 1;
1776  *wt = (anc_time[1] - s_time) / (anc_time[1] - anc_time[0]);
1777  }
1778  } else if (anc_f_stat == ANC_STAT_2T_END) {
1779  if (s_time < anc_time[0]) {
1780  /* outside on the low end, use 1st t only */
1781  *t_interp = 0;
1782  data1_ix = 0;
1783  } else if (s_time > anc_time[2]) {
1784  /* beyond the high end, use 2nd time only */
1785  *t_interp = 0;
1786  data1_ix = 2;
1787  /* anc_t_warn = 1 */
1788  } else {
1789  /* between the MET 1 and 3 */
1790  *t_interp = 1;
1791  data1_ix = 0;
1792  data2_ix = 2;
1793  *wt = (anc_time[2] - s_time) / (anc_time[2] - anc_time[0]);
1794  }
1795  } else if (anc_f_stat == ANC_STAT_3T) {
1796  if (s_time < anc_time[0]) {
1797  *t_interp = 0;
1798  data1_ix = 0;
1799  /* anc_t_warn = 1 */
1800  } else if (s_time > anc_time[2]) {
1801  *t_interp = 0;
1802  data1_ix = 2;
1803  /* anc_t_warn = 1 */
1804  } else if (s_time < anc_time[1]) {
1805  /* between times 0 and 1 */
1806  *t_interp = 1;
1807  data1_ix = 0;
1808  data2_ix = 1;
1809  *wt = (anc_time[1] - s_time) / (anc_time[1] - anc_time[0]);
1810  } else {
1811  /* what's left: between data 1 and 2 */
1812  *t_interp = 1;
1813  data1_ix = 1;
1814  data2_ix = 2;
1815  *wt = (anc_time[2] - s_time) / (anc_time[2] - anc_time[1]);
1816  }
1817  } else {
1818  /* this should not happen at this time - a status that is either
1819  climatology or not defined */
1820  printf("%s, %d: Undefined anc_f_stat - should not happen\n",
1821  __FILE__, __LINE__);
1822  return -1;
1823  }
1824  *data_ix = data1_ix;
1825  *(data_ix + 1) = data2_ix;
1826 
1827  return 0;
1828 }
1829 
1830 int32_t anc_acq_read_gmao(char *file, char *group, char *ds_name,
1831  float **data, unsigned char **qa, double *time, int32_t *nlon, int32_t *nlat,
1832  int32_t *nlvl, double **lon_coord, double **lat_coord)
1833 /*******************************************************************
1834  anc_acq_read_gmao
1835 
1836  purpose: read in the (OBPG modified or not) GMAO parameter from the
1837  specified group. Also, add a column so longitude runs -180 -> 180
1838 
1839  Returns type: int32_t - status 0 is good
1840 
1841  Parameters: (in calling order)
1842  Type Name I/O Description
1843  ---- ---- --- -----------
1844  char * file I GMAO single-level file
1845  char * group I group to set, if "", the
1846  dataset is at the top
1847  char * ds_name I name of dataset to read
1848  float ** data O array of parameter data,
1849  BAD_FLT for bad values
1850  unsigned char ** qa O quality of data 0 good, 1 bad
1851  double * time O time for this data
1852  int32_t * nlon O # longitudes
1853  int32_t * nlat O # latitudes
1854  int32_t * nlvl O # levels
1855  double ** lon_coord O array of longitude values
1856  double ** lat_coord O array of latitude values
1857 
1858  Modification history:
1859  Programmer Date Description of change
1860  ---------- ---- ---------------------
1861  W. Robinson 7 May 2018 Original development
1862 
1863  Note that the level coordinates are not returned (though they could easily be)
1864  - the GMAO data has a fairly fixed definition of what the pressure levels is
1865 
1866  *******************************************************************/ {
1867  static char lcl_fil[FILENAME_MAX] = "";
1868  static int ncid;
1869  int var_id, var2_id, var3_id, dim_id;
1870  int32_t loc;
1871  int32_t nv, ilvl, nlon_pre, il, ip;
1872  size_t tlvl;
1873  float *data_tmp, fillv, missv;
1874 
1875  /* if input file is different, close old file (if open) and open new file */
1876  if (strcmp(file, lcl_fil) != 0) {
1877  if (lcl_fil[0] != 0) {
1878  if (nc_close(ncid) != NC_NOERR) {
1879  fprintf(stderr, "-E- %s %d: nc_close error, file: %s\n",
1880  __FILE__, __LINE__, file);
1881  return 1;
1882  }
1883  }
1884  if (nc_open(file, 0, &ncid) != NC_NOERR) {
1885  fprintf(stderr,
1886  "-E- %s %d: file: %s is not netcdf, not acceptable GMAO file\n",
1887  __FILE__, __LINE__, file);
1888  return -1;
1889  }
1890  strcpy(lcl_fil, file);
1891  }
1892  /* get the main dimensions */
1893  if (((*nlat = ncio_dim_siz(ncid, "lat")) == -1) ||
1894  ((nlon_pre = ncio_dim_siz(ncid, "lon")) == -1)) {
1895  fprintf(stderr,
1896  "-E- %s %d: file: %s error reading lat, lon dimension size\n",
1897  __FILE__, __LINE__, file);
1898  return 1;
1899  }
1900  /* separate level treatment */
1901  if (nc_inq_dimid(ncid, "lev", &dim_id) != NC_NOERR) {
1902  *nlvl = 1;
1903  } else {
1904  if (nc_inq_dimlen(ncid, dim_id, &tlvl) != NC_NOERR) {
1905  fprintf(stderr,
1906  "-E- %s %d: file: %s error reading level dimension size\n",
1907  __FILE__, __LINE__, file);
1908  return 1;
1909  }
1910  *nlvl = tlvl;
1911  }
1912  *nlon = nlon_pre + 1;
1913  nv = *nlat * *nlon * *nlvl;
1914 
1915  /* get the time for this data */
1916  char isodate[32];
1917  memset( isodate, '\0', 32 );
1918  nc_get_att_text(ncid, NC_GLOBAL, "time_coverage_start", isodate);
1919  *time = isodate2unix((const char*) isodate);
1920 
1921  /* set up data and lon, lat arrays using the sizes */
1922  if (((*data = (float *) malloc(nv * sizeof ( float))) == NULL) ||
1923  ((*lon_coord = (double *) malloc(*nlon * sizeof ( double))) == NULL) ||
1924  ((*lat_coord = (double *) malloc(*nlat * sizeof ( double))) == NULL) ||
1925  ((*qa = (unsigned char *) calloc(nv, sizeof ( char))) == NULL) ||
1926  ((data_tmp = (float *) malloc(nlon_pre * *nlat * *nlvl * sizeof ( float))) == NULL)) {
1927  fprintf(stderr,
1928  "-E- %s %d: file: %s error allocating gmao parameter allocation\n",
1929  __FILE__, __LINE__, file);
1930  return 1;
1931  }
1932 
1933  /* read the data and the lon, lat values */
1934  if ((nc_inq_varid(ncid, ds_name, &var_id) != NC_NOERR) ||
1935  (nc_inq_varid(ncid, "lat", &var2_id) != NC_NOERR) ||
1936  (nc_inq_varid(ncid, "lon", &var3_id) != NC_NOERR)) {
1937  fprintf(stderr,
1938  "-E- %s %d: file: %s error setting an id for product: %s\n",
1939  __FILE__, __LINE__, file, ds_name);
1940  return 1;
1941  }
1942 
1943  if (((nc_get_var_double(ncid, var2_id, *lat_coord)) != NC_NOERR) ||
1944  ((nc_get_var_double(ncid, var3_id, *lon_coord)) != NC_NOERR)) {
1945  fprintf(stderr,
1946  "-E- %s %d: file: %s error reading the scales and parameter\n",
1947  __FILE__, __LINE__, file);
1948  return 1;
1949  }
1950  /*
1951  * Get the data and its bad values
1952  */
1953  if ((nc_get_att_float(ncid, var_id, "_FillValue", &fillv) != NC_NOERR) ||
1954  (nc_get_att_float(ncid, var_id, "missing_value", &missv) != NC_NOERR)) {
1955  printf("%s, %d: nc_get_att_float error for parm %s\n",
1956  __FILE__, __LINE__, ds_name);
1957  return -1;
1958  }
1959  /*
1960  * read the dataset and replace either fill or missing with bad value
1961  */
1962  if (nc_get_var_float(ncid, var_id, data_tmp) != NC_NOERR) {
1963  printf("%s, %d: nc_get_var_float error for parm %s\n",
1964  __FILE__, __LINE__, ds_name);
1965  return -1;
1966  }
1967  for (ilvl = 0; ilvl < *nlvl; ilvl++) {
1968  for (il = 0; il < *nlat; il++) {
1969  for (ip = 0; ip < *nlon; ip++) {
1970  loc = ip + *nlon * (il + *nlat * ilvl);
1971  if (ip < nlon_pre) {
1972  *(*data + loc) =
1973  *(data_tmp + ip + nlon_pre * (il + *nlat * ilvl));
1974  } else {
1975  *(*data + loc) = *(data_tmp + nlon_pre * (il + *nlat * ilvl));
1976  }
1977  if ((*(*data + loc) == missv) ||
1978  (*(*data + loc) == fillv)) {
1979  *(*data + loc) = BAD_FLT;
1980  *(*qa + loc) = 1;
1981  }
1982  }
1983  }
1984  }
1985 
1986  free(data_tmp);
1987  (*lon_coord)[nlon_pre] = 180.;
1988  /*
1989  * All is set up, return the data
1990  */
1991  return 0;
1992 }
1993 
1994 /*
1995 int32_t anc_acq_fin()
1996  *******************************************************************
1997  anc_acq_fin
1998 
1999  purpose: finish the ancillary processing
2000 
2001  Returns type: int32_t - return status 0 is OK
2002 
2003  Parameters: (in calling order)
2004  Type Name I/O Description
2005  ---- ---- --- -----------
2006  char ** files I list of files of the ancillary
2007 
2008  Modification history:
2009  Programmer Date Description of change
2010  ---------- ---- ---------------------
2011  W. Robinson 17 May 2018 original development
2012 
2013  *******************************************************************/
2014 
2015 /*
2016 Mainly, free the int objects, lat, lon arrays in the interpolation struct and
2017 also the interpolation structure
2018  */
2019 
2020 int32_t anc_acq_ecmwf_init(char **files, char **prm_nm, int n_prm,
2021  int32_t sto_ix)
2022 /*******************************************************************
2023 
2024  anc_acq_ecmwf_init
2025 
2026  purpose: Identify the incoming MET or OZONE file set and if netcdf,
2027  set up the ancillary data so it can be accessed by anc_acq_line.
2028  Otherwise, do nothing and have getanc.c process the NCEP/TOMS data.
2029  For now, only ECMWF netcdf files can be processed which contain
2030  only 1 time and (at least) the parameters listed in prm_nm
2031 
2032  Returns type: int ancillary data identification: 0 - ECMWF data,
2033  1 non-ECMWF data, -1 any trouble checking input anc files
2034 
2035  Parameters: (in calling order)
2036  Type Name I/O Description
2037  ---- ---- --- -----------
2038  char ** files I 3 file name set, either MET
2039  or OZONE source
2040  char ** prm_nm I array of parameter names to read
2041  from the ECMWF
2042  int n_prm I # parameters in prm_nm. Note
2043  that the met parms (n_prm >1)
2044  have 7 fields to read but the
2045  t(2 m) and td(2 m) combine to
2046  make RH so the storage for met
2047  is 1 less (nprm_sto)
2048  int32_t sto_ix I position in storage array to
2049  place result
2050  int32_t anc_typ O type of anc data 0 - ECMWF,
2051  1 - non ECMWF, -1 - problem
2052 
2053  Modification history:
2054  Programmer Date Description of change
2055  ---------- ---- ---------------------
2056  W. Robinson 24-Sep-2013 Original development
2057 
2058  *******************************************************************/ {
2059  int ids, iprm, npix, nlin, ilin, ipix;
2060  int t_days, t_hrs, lon_gt_180, ird;
2061  int dstpix, npix0, nlin0, status;
2062  int npix_ext, nlin_ext, ntim, f12_mtch, f23_mtch, anc_f_stat;
2063  int ncid, iprm_sto, nprm_sto;
2064  float s_lon, lon_step, e_lon, s_lat, lat_step, e_lat, time;
2065  float prm_bad[] = {ANCBAD, ANCBAD, ANCBAD, ANCBAD, ANCBAD, ANCBAD,
2066  ANCBAD, ANCBAD};
2067  float *base_data, *lat, *lon, *comp1, *comp2;
2068  double data_time;
2069  int64_t jd1900;
2070 
2071  nprm_sto = (n_prm > 1) ? n_prm - 1 : n_prm;
2072  /*
2073  * identify the type of files set from their existance and name matching
2074  */
2075  anc_f_stat = anc_acq_f_stat(files, 0, 3);
2076 
2077  f12_mtch = 0;
2078  if ((anc_f_stat == ANC_STAT_2T_END) || (anc_f_stat == ANC_STAT_1T))
2079  f12_mtch = 1;
2080  f23_mtch = 0;
2081  if ((anc_f_stat == ANC_STAT_2T_START) || (anc_f_stat == ANC_STAT_1T))
2082  f23_mtch = 1;
2083 
2084  if (anc_f_stat == ANC_STAT_CLIM) {
2085  /*
2086  printf( "%s, %d I: Assuming standard (not ECMWF) climatology\n",
2087  __FILE__, __LINE__ );
2088  */
2089  return 1;
2090  }
2091  jd1900 = jd4713bc_get_jd(1900, 1, 1);
2092  for (ids = 0; ids < 3; ids++) {
2093  /*
2094  * ECMWF is only anc data in netcdf format for now, so if it opens
2095  * it should be ECMWF format
2096  */
2097  if ((ids > 0) && (anc_f_stat == ANC_STAT_1T)) break;
2098  if ((ids == 1) && (f12_mtch == 1)) continue;
2099  if ((ids == 2) && (f23_mtch == 1)) break;
2100 
2101  if(Hishdf(files[ids]))
2102  status = NC2_ERR;
2103  else
2104  status = nc_open(files[ids], 0, &ncid);
2105 
2106  if (status != NC_NOERR) {
2107  if (ids == 0) {
2108  /*
2109  printf(
2110  "%s, %d E: file: %s ECMWF is not openable\n",
2111  __FILE__, __LINE__, files[ids] );
2112  */
2113  return -1;
2114  } else {
2115  printf("%s, %d: nc_open failed on file: %s\n", __FILE__, __LINE__,
2116  files[ids]);
2117  printf(" Mismatch or bad ECMWF file\n");
2118  return -1;
2119  }
2120  }
2121 
2122  /*
2123  * get the basic information for MET1
2124  */
2125  if (((npix0 = ncio_dim_siz(ncid, "longitude")) == -1) ||
2126  ((nlin0 = ncio_dim_siz(ncid, "latitude")) == -1) ||
2127  ((ntim = ncio_dim_siz(ncid, "time")) == -1)) {
2128  printf(
2129  "%s, %d: ncio_dim_siz error reading longitude, latitude or time datasets\n",
2130  __FILE__, __LINE__);
2131  return -1;
2132  }
2133  if (ids == 0) {
2134  npix = npix0;
2135  nlin = nlin0;
2136  } else {
2137  if ((npix != npix0) || (nlin != nlin0)) {
2138  printf("%s, %d: mismatch in size of MET array[%d]\n", __FILE__,
2139  __LINE__, ids);
2140  return -1;
2141  }
2142  }
2143  /*
2144  * for now, if more than 1 time, we can't proces it
2145  */
2146  if (ntim > 1) {
2147  printf("%s, %d: Number of times > 1, can't deal with at this time\n",
2148  __FILE__, __LINE__);
2149  return -1;
2150  }
2151  npix_ext = npix + 2;
2152  nlin_ext = nlin + 2;
2153  /*
2154  * allocate storage in the structures for the data
2155  */
2156  for (iprm = 0; iprm < nprm_sto; iprm++) {
2157  if ((met_sto[iprm + sto_ix].data[ids] = (float *)
2158  malloc(npix_ext * nlin_ext * sizeof (float))) == NULL) {
2159  printf("%s, %d: malloc failed for data[%d] in met_sto %d\n", __FILE__,
2160  __LINE__, ids, iprm);
2161  return -1;
2162  }
2163  }
2164  /*
2165  * for 1st dataset, make array to read data into initially and
2166  * arrays for lat, lon
2167  */
2168  if (ids == 0) {
2169  if (((base_data = (float *)
2170  malloc(npix * nlin * sizeof (float))) == NULL) ||
2171  ((comp1 = (float *)
2172  malloc(npix * nlin * sizeof (float))) == NULL) ||
2173  ((comp2 = (float *)
2174  malloc(npix * nlin * sizeof (float))) == NULL)) {
2175  printf("%s, %d: malloc failed for base_data or comp arrays\n",
2176  __FILE__, __LINE__);
2177  return -1;
2178  }
2179 
2180  if ((lat = (float *) malloc(nlin * sizeof (float))) == NULL) {
2181  printf("%s, %d: malloc failed for latitude\n", __FILE__,
2182  __LINE__);
2183  return -1;
2184  }
2185  if ((lon = (float *) malloc(npix * sizeof (float))) == NULL) {
2186  printf("%s, %d: malloc failed for longitude\n", __FILE__,
2187  __LINE__);
2188  return -1;
2189  }
2190  if (ncio_grab_f_ds(ncid, "latitude", lat) != 0) {
2191  printf("%s, %d: ncio_grab_f_ds failed on latitude\n",
2192  __FILE__, __LINE__);
2193  return -1;
2194  }
2195  if (ncio_grab_f_ds(ncid, "longitude", lon) != 0) {
2196  printf("%s, %d: ncio_grab_f_ds failed on longitude\n",
2197  __FILE__, __LINE__);
2198  return -1;
2199  }
2200  /*
2201  * from the latitude, longitude arrays, determine the nav properties
2202  * ECMWF longitudes go 0 -> 360 and we do -180 -> 180
2203  */
2204  s_lat = lat[0];
2205  e_lat = lat[nlin - 1];
2206  lat_step = lat[1] - lat[0];
2207 
2208  lon_gt_180 = -1;
2209  for (ipix = 0; ipix < npix; ipix++) {
2210  if (lon[ipix] > 180.) {
2211  s_lon = lon[ipix] - 360.;
2212  lon_gt_180 = ipix;
2213  e_lon = lon[ipix - 1]; /* if lon_gt_180 = 0, need to upgrade this */
2214  break;
2215  }
2216  }
2217  lon_step = lon[1] - lon[0];
2218  }
2219  /*
2220  * get the time from the dataset and convert to julian days
2221  */
2222  if (ncio_grab_f_ds(ncid, "time", &time) != 0) {
2223  printf("%s, %d: error reading the time in\n", __FILE__, __LINE__);
2224  return -1;
2225  }
2226 
2227  t_days = time / 24;
2228  t_hrs = (int) time % 24;
2229  data_time = (double) (jd1900 + t_days) + (double) t_hrs / 24.;
2230  /*
2231  * for all params, read the data, put lon in range -180 - 180
2232  * and add extra layer to make interpolation easier
2233  * The RH is made from T and Td in else below
2234  */
2235  ird = 0;
2236  for (iprm = 0; iprm < nprm_sto; iprm++) {
2237  iprm_sto = iprm + sto_ix;
2238  if (ird != 5) {
2239  if (ncio_grab_stdsclf_ds(ncid, prm_nm[ ird ],
2240  prm_bad[ird], base_data) != 0) {
2241  printf("%s, %d: ncio_grab_stdsclf_ds failed on %s\n",
2242  __FILE__, __LINE__, prm_nm[ird]);
2243  return -1;
2244  }
2245  ird++;
2246  } else {
2247  if ((ncio_grab_stdsclf_ds(ncid, prm_nm[ird],
2248  prm_bad[ird + sto_ix], comp1) != 0) ||
2249  (ncio_grab_stdsclf_ds(ncid, prm_nm[ird + 1],
2250  prm_bad[ird + sto_ix + 1], comp2) != 0)) {
2251  printf("%s, %d: ncio_grab_stdsclf_ds failed on %s or %s\n",
2252  __FILE__, __LINE__, prm_nm[ird], prm_nm[ird + 1]);
2253  return -1;
2254  }
2255  ird = ird + 2;
2256  /*
2257  * for the td and t at 2 m, we need to make a RH
2258  */
2259  if (met_cvt_ttd_to_rh(npix * nlin, comp1, MET_UNITS__T_K, comp2,
2260  MET_UNITS__T_K, base_data) != 0) {
2261  printf("met_cvt_ttd_to_rh had an error\n");
2262  printf("%s, %d: met_cvt_ttd_to_rh failure\n",
2263  __FILE__, __LINE__);
2264  return -1;
2265  }
2266  }
2267  /* rotate to -180 -> 180 */
2268  for (ilin = 0; ilin < nlin; ilin++) {
2269  for (ipix = 0; ipix < npix; ipix++) {
2270  dstpix = ipix - lon_gt_180; /* put in with lon -180 -> 180 */
2271  if (dstpix < 0) dstpix += npix;
2272  *(met_sto[ iprm_sto ].data[ids]
2273  + dstpix + 1 + (ilin + 1) * npix_ext) =
2274  *(base_data + ipix + npix * ilin);
2275  }
2276  }
2277  /* now, the extra boarder: lat, then lon */
2278  /* for lat, repeat the nearest value */
2279  for (ipix = 0; ipix < npix; ipix++) {
2280  *(met_sto[iprm_sto].data[ids] + ipix + 1) =
2281  *(met_sto[iprm_sto].data[ids] + ipix + 1 + npix_ext);
2282  *(met_sto[iprm_sto].data[ids] + ipix + 1 + (nlin + 1) * npix_ext) =
2283  *(met_sto[iprm_sto].data[ids] + ipix + 1 + nlin * npix_ext);
2284  }
2285  /* for lon, use the opposite side value */
2286  for (ilin = 0; ilin < nlin_ext; ilin++) {
2287  *(met_sto[iprm_sto].data[ids] + ilin * npix_ext) =
2288  *(met_sto[iprm_sto].data[ids] + npix + ilin * npix_ext);
2289  *(met_sto[iprm_sto].data[ids] + npix + 1 + ilin * npix_ext) =
2290  *(met_sto[iprm_sto].data[ids] + 1 + ilin * npix_ext);
2291  }
2292  /* put in the controls found above */
2293  met_sto[iprm_sto].s_lon = s_lon;
2294  met_sto[iprm_sto].lon_step = lon_step;
2295  met_sto[iprm_sto].nlon = npix;
2296  met_sto[iprm_sto].e_lon = e_lon;
2297  met_sto[iprm_sto].s_lat = s_lat;
2298  met_sto[iprm_sto].lat_step = lat_step;
2299  met_sto[iprm_sto].nlat = nlin;
2300  met_sto[iprm_sto].e_lat = e_lat;
2301  met_sto[iprm_sto].data_time[ids] = data_time;
2302  met_sto[iprm_sto].anc_f_stat = anc_f_stat;
2303  }
2304  /*
2305  * close the dataset
2306  */
2307  nc_close(ncid);
2308  }
2309  return 0;
2310 }
2311 
2312 int anc_acq_lin(int32_t anc_class, l1str *l1rec)
2313 /*******************************************************************
2314 
2315  anc_acq_lin
2316 
2317  purpose: get proper ancillary parameters for a particular line of
2318  points
2319 
2320  Returns type: int - 0 if good, else -1
2321 
2322  Parameters: (in calling order)
2323  Type Name I/O Description
2324  ---- ---- --- -----------
2325  int32_t anc_class I anc data class to access the
2326  correct stored grids: 0 for
2327  MET grids and 1 for ozone grid
2328  l1str * l1rec I/O structure with information
2329  for the line
2330  Modification history:
2331  Programmer Date Description of change
2332  ---------- ---- ---------------------
2333  W. Robinson 16 Aug 2013 original development
2334 
2335  *******************************************************************/ {
2336  double l_time;
2337  float data_val, uwnd, vwnd, data_val1, data_val2, dx, dy, lon_frac, lat_frac;
2338  float trg_lon, trg_lat, wt_t1, wt_t2, s_lon, s_lat;
2339  float *data;
2340  int iprm, xbox_st, ybox_st, nx, ny, t_interp, data1_ix, data2_ix, anc_f_stat;
2341  int npix, ipix, sto_st, sto_en;
2342  static float r2d = 57.29577951;
2343  /*
2344  * find places in the met_sto structure to look in
2345  */
2346  if (anc_class == 0) {
2347  sto_st = 0;
2348  sto_en = 5;
2349  } else {
2350  sto_st = 6;
2351  sto_en = 6;
2352  }
2353  /*
2354  * get the time of the current line
2355  */
2356  int16_t year, month, day;
2357  double sec;
2358  unix2ymds(l1rec->scantime, &year, &month, &day, &sec);
2359  l_time = (double) jd4713bc_get_jd((int32_t) year, (int32_t) month, (int32_t) day);
2360  l_time += sec / 86400.0;
2361 
2362  npix = l1rec->npix;
2363  /*
2364  * for this line, decide which of the 3 anc files will be needed based on
2365  * the line's time and the ancillary times
2366  */
2367  /* ***** In this set up, all grids are the same, so only one
2368  determination is needed
2369  */
2370  anc_f_stat = met_sto[sto_st].anc_f_stat;
2371  if (anc_f_stat == ANC_STAT_1T) {
2372  /* use data[0] only */
2373  t_interp = 0;
2374  data1_ix = 0;
2375  /* further along, when interpolating, use met_sto[0].data[data1_ix]
2376  to access the data */
2377  } else if (anc_f_stat == ANC_STAT_2T_START) {
2378  /* 2 different times in the 1, 2 positions */
2379  if (l_time < met_sto[sto_st].data_time[0]) {
2380  printf("%s, %d: data time is before the ancillary data start time\n",
2381  __FILE__, __LINE__);
2382  return -1;
2383  } else if (l_time > met_sto[sto_st].data_time[1]) {
2384  /* use MET2 only */
2385  t_interp = 0;
2386  data1_ix = 1;
2387  } else {
2388  /* in-between MET 1, 2 use data 0, 1 and time interpolate */
2389  t_interp = 1;
2390  data1_ix = 0;
2391  data2_ix = 1;
2392  wt_t1 = (met_sto[sto_st].data_time[1] - l_time) /
2393  (met_sto[sto_st].data_time[1] - met_sto[sto_st].data_time[0]);
2394  wt_t2 = (l_time - met_sto[sto_st].data_time[0]) /
2395  (met_sto[sto_st].data_time[1] - met_sto[sto_st].data_time[0]);
2396  }
2397  } else if (anc_f_stat == ANC_STAT_2T_END) {
2398  if (l_time < met_sto[sto_st].data_time[0]) {
2399  /* outside on the low end, use data[0] */
2400  t_interp = 0;
2401  data1_ix = 0;
2402  } else if (l_time > met_sto[sto_st].data_time[2]) {
2403  /* beyond the high end, Can't use end time alone */
2404  printf("%s, %d: data time is after the ancillary data end time\n",
2405  __FILE__, __LINE__);
2406  return -1;
2407  } else {
2408  /* between the MET 1 and 3 */
2409  t_interp = 1;
2410  data1_ix = 0;
2411  data2_ix = 2;
2412  wt_t1 = (met_sto[sto_st].data_time[2] - l_time) /
2413  (met_sto[sto_st].data_time[2] - met_sto[sto_st].data_time[0]);
2414  wt_t2 = (l_time - met_sto[sto_st].data_time[0]) /
2415  (met_sto[sto_st].data_time[2] - met_sto[sto_st].data_time[0]);
2416  }
2417  } else if (anc_f_stat == ANC_STAT_3T) {
2418  if (l_time < met_sto[sto_st].data_time[0]) {
2419  printf("%s, %d: data time is before the ancillary data start time\n",
2420  __FILE__, __LINE__);
2421  return -1;
2422  } else if (l_time > met_sto[sto_st].data_time[2]) {
2423  printf("%s, %d: data time is after the ancillary data end time\n",
2424  __FILE__, __LINE__);
2425  return -1;
2426  } else if (l_time < met_sto[sto_st].data_time[1]) {
2427  /* between data 0 and 1 */
2428  t_interp = 1;
2429  data1_ix = 0;
2430  data2_ix = 1;
2431  wt_t1 = (met_sto[sto_st].data_time[1] - l_time) /
2432  (met_sto[sto_st].data_time[1] - met_sto[sto_st].data_time[0]);
2433  wt_t2 = (l_time - met_sto[sto_st].data_time[0]) /
2434  (met_sto[sto_st].data_time[1] - met_sto[sto_st].data_time[0]);
2435  } else {
2436  /* what's left: between data 1 and 2 */
2437  t_interp = 1;
2438  data1_ix = 1;
2439  data2_ix = 2;
2440  wt_t1 = (met_sto[sto_st].data_time[2] - l_time) /
2441  (met_sto[sto_st].data_time[2] - met_sto[sto_st].data_time[1]);
2442  wt_t2 = (l_time - met_sto[sto_st].data_time[1]) /
2443  (met_sto[sto_st].data_time[2] - met_sto[sto_st].data_time[1]);
2444  }
2445  } else {
2446  /* this should not happen at this time - a status that is either
2447  climatology or not defined */
2448  printf("%s, %d: Undefined anc_f_stat - should not happen\n",
2449  __FILE__, __LINE__);
2450  return -1;
2451  }
2452  /*
2453  * this found if time interpolation is needed and the and grids to use.
2454  * next, for each pixel, find the bounding grid box and weights
2455  * for bi-linear interpolation to be applied to each parameter
2456  *
2457  * AGAIN note that since all parameters are from the same size grid,
2458  * we can compute this information once.
2459  */
2460  dx = met_sto[sto_st].lon_step;
2461  dy = met_sto[sto_st].lat_step;
2462  s_lon = met_sto[sto_st].s_lon;
2463  s_lat = met_sto[sto_st].s_lat;
2464  nx = met_sto[sto_st].nlon;
2465  ny = met_sto[sto_st].nlat;
2466  for (ipix = 0; ipix < npix; ipix++) {
2467  trg_lat = l1rec->lat[ipix];
2468  trg_lon = l1rec->lon[ipix];
2469 
2470  /*
2471  xbox_st =
2472  MAX( MIN( (INT)( ( trg_lon - s_lon + dx / 2. ) / dx ), nx + 1 ), 0 );
2473  ybox_st =
2474  MAX( MIN( (INT)( ( trg_lat - s_lat + dy / 2. ) / dy ), ny + 1 ), 0 );
2475  x_dist = xbox_st * dx + s_lon - dx / 2;
2476  y_dist = ybox_st * dy + s_lat - dy / 2;
2477 
2478  I think below is correct for data at the grid points
2479  */
2480  xbox_st = MAX(MIN((int) ((trg_lon - s_lon + dx) / dx), nx + 1), 0);
2481  ybox_st = MAX(MIN((int) ((trg_lat - s_lat + dy) / dy), ny + 1), 0);
2482 
2483  lon_frac = (trg_lon - s_lon) / dx - (float) (xbox_st - 1);
2484  lat_frac = (trg_lat - s_lat) / dy - (float) (ybox_st - 1);
2485 
2486  for (iprm = sto_st; iprm < (sto_en + 1); iprm++) {
2487  data = met_sto[iprm].data[data1_ix];
2488  data_val1 = bilin_interp(data, xbox_st, (nx + 2), ybox_st, lon_frac,
2489  lat_frac);
2490 
2491  if (t_interp == 1) {
2492  data = met_sto[iprm].data[data2_ix];
2493  data_val2 = bilin_interp(data, xbox_st, (nx + 2), ybox_st, lon_frac,
2494  lat_frac);
2495 
2496  /*
2497  * do time interpolation
2498  */
2499  if (data_val1 < ANCBAD + 1) {
2500  if (data_val2 < ANCBAD + 1)
2501  data_val = ANCBAD;
2502  else
2503  data_val = data_val2;
2504  } else
2505  data_val = wt_t1 * data_val1 + wt_t2 * data_val2;
2506  } else
2507  data_val = data_val1;
2508  /*
2509  * place this interpolated value in proper l1rec slot
2510  */
2511  switch (iprm) {
2512  case 0: /* sfc press */
2513  /* Currently, no use for this, but it may be better than what
2514  is done with MSL pressure to take it to a height above sea
2515  level. USE_PMSL of 0 will use this. Note that pressure on Mt
2516  Everest is nominally 337 mb, so enlarge range accordingly */
2517  if (USE_PMSL == 1)
2518  break;
2519  else {
2520  if (input->pressure != -2000) {
2521  data_val = (data_val < ANCBAD + 1) ? ANCBAD : data_val / 100.;
2522  if (data_val < 0) l1rec->pr[ipix] = 1013.25;
2523  else if (data_val < 250.) l1rec->pr[ipix] = 250.;
2524  else if (data_val > 1100.) l1rec->pr[ipix] = 1100.;
2525  else l1rec->pr[ipix] = data_val;
2526  }
2527  }
2528  break;
2529  case 1: /* precip water */
2530  /* need to make from kg m^-2 into g cm^-2 */
2531  if (input->watervapor != -2000)
2532  l1rec->wv[ipix] = (data_val < ANCBAD + 1) ? 0. : data_val / 10.;
2533  break;
2534  case 2: /* sea level pressure */
2535  /* need to make from pascals to hectopascals (millibars) */
2536  if (USE_PMSL == 0)
2537  break;
2538  else {
2539  if (input->pressure != -2000) {
2540  data_val = (data_val < ANCBAD + 1) ? ANCBAD : data_val / 100.;
2541  if (data_val < 0) l1rec->pr[ipix] = 1013.25;
2542  else if (data_val < 900.) l1rec->pr[ipix] = 900.;
2543  else if (data_val > 1100.) l1rec->pr[ipix] = 1100.;
2544  else l1rec->pr[ipix] = data_val;
2545  }
2546 
2547  /* if processing land, adjust pressure for terrain height */
2548  if (proc_land && l1rec->height[ipix] != 0.0)
2549  l1rec->pr[ipix] *= exp(-l1rec->height[ipix] / 8434);
2550  }
2551  break;
2552  case 3: /* u wind, zonal W-E */
2553  uwnd = (data_val < ANCBAD + 1) ? 0. : data_val;
2554  l1rec->zw[ipix] = uwnd;
2555  break;
2556  case 4: /* v wind, meridional S-N */
2557  vwnd = (data_val < ANCBAD + 1) ? 0. : data_val;
2558  l1rec->mw[ipix] = vwnd;
2559  if (input->windspeed != -2000)
2560  l1rec->ws[ipix] = sqrt(pow(uwnd, 2.) + pow(vwnd, 2.));
2561  if (input->windangle != -2000)
2562  l1rec->wd[ipix] = atan2f(-uwnd, vwnd) * r2d;
2563  break;
2564  case 5: /* rel humidity % */
2565  if (input->relhumid != -2000)
2566  l1rec->rh[ipix] = (data_val < ANCBAD + 1) ? 0. : data_val;
2567  break;
2568  case 6: /* ozone */
2569  if (input->ozone != -2000) {
2570  /* convert from kg m^-2 to DU (units of 10 um STP OZ thickness)
2571  to cm of thickness */
2572  l1rec->oz[ipix] = (data_val < ANCBAD + 1) ?
2573  0. : data_val * OZ_KG_M2_TO_DU / 1000.;
2574  }
2575  break;
2576  }
2577  }
2578  }
2579  return 0;
2580 }
2581 
2582 int anc_acq_lin_olci(int anc_class, char *file, l1str *l1rec)
2583 /*******************************************************************
2584 
2585  anc_acq_lin_olci
2586 
2587  purpose: Read a line of ancillary data from the OLCI tie point meteo
2588  file and fill the met or ozone fields of teh l1 structure
2589 
2590  Returns type: status - 0 if all is good
2591 
2592  Parameters: (in calling order)
2593  Type Name I/O Description
2594  ---- ---- --- -----------
2595  int anc_class I class of ancillary: 0 - met,
2596  1 - ozone
2597  char * file I ancillary file name
2598  l1str * l1rec I/O structure to fill with data
2599 
2600  Modification history:
2601  Programmer Date Description of change
2602  ---------- ---- ---------------------
2603  W. Robinson, SAIC 24 Feb 2016 original development
2604 
2605  *******************************************************************/ {
2606  static int ncid[2], varid[5]; /* there are 2 file ids for the
2607  met and ozone tie point files, and 4
2608  met dataset or variable ids:
2609  0 - wind - horizontal_wind
2610  1 - rh - humidity
2611  2 - msl pressure - sea_level_pressure
2612  3 - precip water - total_column_water_vapor
2613  4 - ozone - total_ozone
2614  */
2615  static float fill[5], valid_min[5], valid_max[5];
2616  static int32_t firstcalls[2] = {1, 1}, pix_smp[2];
2617  static int32_t npix, *qual, *qual_met, *qual_oz;
2618  static int32_t tie_epix, spix, epix, dpix;
2619  static size_t tie_nlin[2], tie_npix[2];
2620  static float *tie_data, *tie_met, *tie_oz;
2621  static float r2d = 57.29577951;
2622  int32_t anc_field_per_parm[5] = {1, 2, 1, 1, 1};
2623  int32_t anc_class_n_ds[2] = {4, 1}, class_off, n_ds_prm, ptr_prm;
2624  int32_t n_field_per_parm, ifld, ipx, px_tie1, px_tie2;
2625  int nid, lin_smp, ids, stat, ipx_dat;
2626  size_t start[3], count[3];
2627  char *anc_prm_nm[] = {"humidity", "horizontal_wind", "sea_level_pressure",
2628  "total_columnar_water_vapour", "total_ozone"};
2629  float *ar, fr_dist, w1, w2;
2630 
2631  int32_t iscan = l1rec->iscan;
2632  /*
2633  * do the initialization on the first calls for met and ozone files
2634  */
2635  if (anc_class == 0)
2636  class_off = 0;
2637  else
2638  class_off = 4;
2639 
2640  n_ds_prm = anc_class_n_ds[anc_class];
2641 
2642  if (firstcalls[anc_class]) {
2643  npix = l1rec->npix;
2644  spix = l1_input->spixl - 1;
2645  epix = l1_input->epixl - 1;
2646  dpix = l1_input->dpixl;
2647  /*
2648  * open the file and check that the sampling is 64 in pixls, 1 in lines
2649  */
2650  if (nc_open(file, NC_NOWRITE, (ncid + anc_class)) != NC_NOERR) {
2651  fprintf(stderr,
2652  "-E- %s %d: Unable to open OLCI tie point anc file: %s\n",
2653  __FILE__, __LINE__, file);
2654  return -1;
2655  }
2656  if (nc_get_att_int(ncid[anc_class], NC_GLOBAL, "ac_subsampling_factor",
2657  &pix_smp[anc_class]) != NC_NOERR) {
2658  fprintf(stderr,
2659  "-E- %s %d: Unable to read column sampling attrib from OLCI tie point anc file: %s\n",
2660  __FILE__, __LINE__, file);
2661  return -1;
2662  }
2663  if (nc_get_att_int(ncid[anc_class], NC_GLOBAL, "al_subsampling_factor",
2664  &lin_smp) != NC_NOERR) {
2665  fprintf(stderr,
2666  "-E- %s %d: Unable to read row sampling attrib from OLCI tie point anc file: %s\n",
2667  __FILE__, __LINE__, file);
2668  return -1;
2669  }
2670  /*
2671  * get the tie dataset dim sizes
2672  */
2673  if (nc_inq_dimid(ncid[anc_class], "tie_columns", &nid) != NC_NOERR) {
2674  fprintf(stderr,
2675  "-E- %s %d: Unable to get OLCI tie point meteo # pixels\n",
2676  __FILE__, __LINE__);
2677  return -1;
2678  }
2679  if (nc_inq_dimlen(ncid[anc_class], nid, &tie_npix[anc_class])
2680  != NC_NOERR) {
2681  fprintf(stderr,
2682  "-E- %s %d: Unable to get OLCI tie point meteo # pixels\n",
2683  __FILE__, __LINE__);
2684  return -1;
2685  }
2686 
2687  if (nc_inq_dimid(ncid[anc_class], "tie_rows", &nid) != NC_NOERR) {
2688  fprintf(stderr,
2689  "-E- %s %d: Unable to get OLCI tie point meteo # lines\n",
2690  __FILE__, __LINE__);
2691  return -1;
2692  }
2693  if (nc_inq_dimlen(ncid[anc_class], nid, &tie_nlin[anc_class]) != NC_NOERR) {
2694  fprintf(stderr,
2695  "-E- %s %d: Unable to get OLCI tie point meteo # lines\n",
2696  __FILE__, __LINE__);
2697  return -1;
2698  }
2699  /*
2700  * If not the expected values, leave
2701  */
2702  if (lin_smp != 1) {
2703  fprintf(stderr,
2704  "-E- %s %d: OLCI and tie point line sampling: %d, not = 1\n",
2705  __FILE__, __LINE__, lin_smp);
2706  return -1;
2707  }
2708  tie_epix = pix_smp[anc_class] * (tie_npix[anc_class] - 1) + 1;
2709  if (tie_epix < epix) {
2710  fprintf(stderr,
2711  "-E- %s %d: tie point range out to pixel %d is < data range of %d\n",
2712  __FILE__, __LINE__, tie_epix, epix);
2713  fprintf(stderr,
2714  "tie point sampling: %d and # pixels: %d\n",
2715  pix_smp[anc_class], (int) tie_npix[anc_class]);
2716  return -1;
2717  }
2718  /*
2719  * warn if the # lines in tie and dataset don't match
2720  * UNFORTUNATELY, the # scans is not set in that data area at this time
2721  olci_dat = (olci_t *) ( l1rec->l1file->private_data );
2722  nlin = olci_dat->nscan;
2723  if( tie_nlin[anc_class] != nlin )
2724  fprintf(stderr,
2725  "-W- %s %d: OLCI tie point and radiance data have unequal # lines\n",
2726  __FILE__, __LINE__ );
2727  */
2728  /*
2729  * set up storage for the tie point data and quality line
2730  */
2731  if ((tie_data = (float *) malloc(tie_npix[anc_class] * sizeof (float)))
2732  == NULL) {
2733  fprintf(stderr,
2734  "-E- %s %d: Unable to allocate tie point data array\n",
2735  __FILE__, __LINE__);
2736  return -1;
2737  }
2738  /* assign address for tie data storage */
2739  if (anc_class == 0) tie_met = tie_data;
2740  else tie_oz = tie_data;
2741 
2742  if ((qual = (int32_t *) malloc(tie_npix[anc_class] * sizeof (int32_t)))
2743  == NULL) {
2744  fprintf(stderr,
2745  "-E- %s %d: Unable to allocate tie point quality array\n",
2746  __FILE__, __LINE__);
2747  return -1;
2748  }
2749  if (anc_class == 0) qual_met = qual;
2750  else qual_oz = qual;
2751  /*
2752  * loop through the parameters and set the access id, get fill value,
2753  * valid min and max
2754  */
2755  for (ids = 0; ids < n_ds_prm; ids++) {
2756  ptr_prm = ids + class_off;
2757  if (nc_inq_varid(ncid[anc_class], anc_prm_nm[ptr_prm], &varid[ptr_prm])
2758  != NC_NOERR) {
2759  fprintf(stderr,
2760  "-E- %s %d: Can't find id for OLCI anc tie point dataset: %s\n",
2761  __FILE__, __LINE__, anc_prm_nm[ptr_prm]);
2762  return -1;
2763  }
2764  if (nc_get_att_float(ncid[anc_class], varid[ptr_prm], "_FillValue",
2765  &fill[ptr_prm]) != NC_NOERR) {
2766  fprintf(stderr,
2767  "-E- %s %d: Can't get _FillValue for OLCI anc tie point dataset: %s\n",
2768  __FILE__, __LINE__, anc_prm_nm[ptr_prm]);
2769  return -1;
2770  }
2771 
2772  if (nc_get_att_float(ncid[anc_class], varid[ptr_prm], "valid_min",
2773  &valid_min[ptr_prm]) != NC_NOERR) {
2774  fprintf(stderr,
2775  "-E- %s %d: Can't get valid_min for OLCI anc tie point dataset: %s\n",
2776  __FILE__, __LINE__, anc_prm_nm[ptr_prm]);
2777  return -1;
2778  }
2779  if (nc_get_att_float(ncid[anc_class], varid[ptr_prm], "valid_max",
2780  &valid_max[ptr_prm]) != NC_NOERR) {
2781  fprintf(stderr,
2782  "-E- %s %d: Can't get valid_max for OLCI anc tie point dataset: %s\n",
2783  __FILE__, __LINE__, anc_prm_nm[ptr_prm]);
2784  return -1;
2785  }
2786  }
2787  firstcalls[anc_class] = 0;
2788  } /* end of init portion */
2789  /*
2790  * read the line of each parameter
2791  */
2792  start[1] = 0;
2793  count[0] = 1;
2794 
2795  for (ids = 0; ids < n_ds_prm; ids++) {
2796  ptr_prm = ids + class_off;
2797  /* note that the description of the horizontal wind in the tie point
2798  file does not indicate the u, v components. from all descriptions
2799  of wind components that ecmwf has (the source), they mention u,
2800  then v. We can only assume that the 1st index is u (zonal)
2801  followed by v (meridional) */
2802  switch (ptr_prm) {
2803  case 0: ar = l1rec->rh;
2804  break;
2805  case 1: ar = l1rec->zw;
2806  break;
2807  case 2: ar = l1rec->pr;
2808  break;
2809  case 3: ar = l1rec->wv;
2810  break;
2811  case 4: ar = l1rec->oz;
2812  break;
2813  }
2814  /*
2815  * set the start and count to read the current scan
2816  */
2817  start[0] = iscan;
2818  count[1] = tie_npix[anc_class];
2819  count[2] = 0;
2820  start[2] = 0;
2821 
2822  tie_data = (anc_class == 0) ? tie_met : tie_oz;
2823  qual = (anc_class == 0) ? qual_met : qual_oz;
2824 
2825  n_field_per_parm = anc_field_per_parm[ptr_prm];
2826  for (ifld = 0; ifld < n_field_per_parm; ifld++) {
2827  if (ids == 1) {
2828  count[2] = 1;
2829  if (ifld == 1) {
2830  start[2] = 1; /* for 2nd wind field */
2831  ar = l1rec->mw;
2832  }
2833  }
2834  /* make sure we have not exceeded the tie point # lines */
2835  if (iscan < tie_nlin[anc_class]) {
2836  /* read the line */
2837  if ((stat = nc_get_vara_float(ncid[anc_class], varid[ptr_prm],
2838  start, count, tie_data)) != NC_NOERR) {
2839  fprintf(stderr,
2840  "-E- %s %d: Can't read OLCI anc tie point line, dataset: %s\n",
2841  __FILE__, __LINE__, anc_prm_nm[ptr_prm]);
2842  return -1;
2843  }
2844  /* set a quality for all the tie points and convert to expected units */
2845  for (ipx = 0; ipx < tie_npix[anc_class]; ipx++) {
2846  if ((*(tie_data + ipx) == fill[ptr_prm]) ||
2847  (*(tie_data + ipx) < valid_min[ptr_prm]) ||
2848  (*(tie_data + ipx) > valid_max[ptr_prm]))
2849  *(qual + ipx) = 1;
2850  else {
2851  *(qual + ipx) = 0;
2852  if (ptr_prm == 3) *(tie_data + ipx) *= 0.1; /* for pw conversion
2853  kg m-2 -> g cm-2 */
2854  if (ptr_prm == 4) *(tie_data + ipx) *= 46.698; /* for oz
2855  conversion kg m-2 -> cm std atmos */
2856  }
2857  }
2858  /* interpolate */
2859  for (ipx = 0; ipx < npix; ipx++) {
2860  /* need to find pixel in un-subsetted line */
2861  ipx_dat = ipx * dpix + spix;
2862  px_tie1 = ipx_dat / pix_smp[anc_class];
2863  px_tie2 = px_tie1 + 1;
2864  /* for very last pixel, we need to make this adjustment*/
2865  if (px_tie2 > (tie_npix[anc_class] - 1)) {
2866  px_tie1 -= 1;
2867  px_tie2 -= 1;
2868  }
2869  if ((*(qual + px_tie1) == 1) || (*(qual + px_tie2) == 1)) {
2870  /* fill with filler */
2871  *(ar + ipx) = anc_miss_fill(ptr_prm);
2872  l1rec->flags[ipx] |= ATMWARN;
2873  } else {
2874  /* interpolate */
2875  fr_dist = (float) (ipx_dat - px_tie1 * pix_smp[anc_class]) /
2876  (float) pix_smp[anc_class];
2877  *(ar + ipx) = tie_data[px_tie1] * (1. - fr_dist) +
2878  tie_data[px_tie2] * fr_dist;
2879  }
2880  }
2881  } else {
2882  /* place a fill value and set flag to atmwarn for whole line */
2883  for (ipx = 0; ipx < npix; ipx++) {
2884  *(ar + ipx) = anc_miss_fill(ptr_prm);
2885  l1rec->flags[ipx] |= ATMWARN;
2886  }
2887  }
2888  }
2889  /*
2890  * set up wind speed, direction
2891  */
2892  if (ids == 1) {
2893  for (ipx = 0; ipx < npix; ipx++) {
2894  w1 = l1rec->zw[ipx];
2895  w2 = l1rec->mw[ipx];
2896  if (input->windspeed != -2000)
2897  l1rec->ws[ipx] = sqrt(w1 * w1 + w2 * w2);
2898  if (input->windangle != -2000)
2899  l1rec->wd[ipx] = atan2f(-w1, w2) * r2d;
2900  }
2901  }
2902  /*
2903  * adjust the pressure over land
2904  */
2905  if (ids == 2) {
2906  for (ipx = 0; ipx < npix; ipx++) {
2907  if (l1rec->height[ipx] != 0.0)
2908  l1rec->pr[ipx] *= exp(-l1rec->height[ipx] / 8434);
2909  }
2910  }
2911  }
2912  return 0;
2913 }
2914 
2915 float anc_miss_fill(int32_t prod_ix)
2916 /*******************************************************************
2917 
2918  anc_miss_fill
2919 
2920  purpose: return a ancillary fill value for a product type
2921  using the same fillers as getanc.c/interpolate()
2922 
2923  Returns type: float of the fill value
2924 
2925  Parameters: (in calling order)
2926  Type Name I/O Description
2927  ---- ---- --- -----------
2928  int23_t prod_ix I product number being done
2929 
2930  Modification history:
2931  Programmer Date Description of change
2932  ---------- ---- ---------------------
2933  W. Robinson 16 Aug 2013 original development
2934 
2935  *******************************************************************/ {
2936  float fill;
2937 
2938  switch (prod_ix) {
2939  case 0: fill = 80.;
2940  break; /* humidity % */
2941  case 1: fill = 6.;
2942  break; /* wind m s-1 */
2943  case 2: fill = 1013.;
2944  break; /* msl pressure hPa */
2945  case 3: fill = 50.;
2946  break; /* pw in kg m-2 */
2947  case 4: fill = 0.36;
2948  break; /* ozone in cm at std atmos */
2949  }
2950  return fill;
2951 }
2952 
2953 float bilin_interp(float *data, int xbox_st, int nx, int ybox_st,
2954  float xfrac, float yfrac)
2955 /*******************************************************************
2956 
2957  bilin_interp
2958 
2959  purpose: quick bi-linear interpolation.
2960 
2961  Returns type: float of interpolated result
2962 
2963  Parameters: (in calling order)
2964  Type Name I/O Description
2965  ---- ---- --- -----------
2966  float * data I 2-d grid of data
2967  int xbox_st I x (longitude) index of
2968  grid box to interpolate
2969  int nx I # pixels in x
2970  int ybox_st I x (latitude) index of
2971  grid box to interpolate
2972  float xfrac I fractional grid box distance
2973  from xbox_st to the point
2974  float yfrac I fractional grid box distance
2975  from ybox_st to the point
2976 
2977  Modification history:
2978  Programmer Date Description of change
2979  ---------- ---- ---------------------
2980  W. Robinson 16 Aug 2013 original development
2981 
2982  *******************************************************************/ {
2983  float data_val;
2984 
2985  if ((*(data + xbox_st + nx * ybox_st) < (ANCBAD + 1)) ||
2986  (*(data + xbox_st + nx * (ybox_st + 1)) < (ANCBAD + 1)) ||
2987  (*(data + (xbox_st + 1) + nx * ybox_st) < (ANCBAD + 1)) ||
2988  (*(data + (xbox_st + 1) + nx * (ybox_st + 1)) < (ANCBAD + 1)))
2989  data_val = ANCBAD;
2990  else
2991  data_val =
2992  (1 - xfrac) * (1 - yfrac) *
2993  *(data + xbox_st + nx * ybox_st) +
2994  (1 - xfrac) * yfrac *
2995  *(data + xbox_st + nx * (ybox_st + 1)) +
2996  xfrac * (1 - yfrac) *
2997  *(data + (xbox_st + 1) + nx * ybox_st) +
2998  xfrac * yfrac *
2999  *(data + (xbox_st + 1) + nx * (ybox_st + 1));
3000  return data_val;
3001 }
3002 
3003 int64_t jd4713bc_get_jd(int32_t year, int32_t month, int32_t day)
3004 /*******************************************************************
3005 
3006  jd4713bc_get_jd
3007 
3008  purpose: get the julian day (from 4713 BC) for the year,
3009  month and day of month
3010  taken from the idl jd routine
3011 
3012  Returns type: int64_t - the julian date
3013 
3014  Parameters: (in calling order)
3015  Type Name I/O Description
3016  ---- ---- --- -----------
3017  int32_t year I standard 4 digit year
3018  int32_t month I month of the year
3019  int32_t day I day of month
3020 
3021  Modification history:
3022  Programmer Date Description of change
3023  ---------- ---- ---------------------
3024  W. Robinson 5 Aug 2013 original development
3025 
3026  *******************************************************************/ {
3027  int64_t lyear, lmonth, lday, jday;
3028 
3029  lyear = (int64_t) year;
3030  lmonth = (int64_t) month;
3031  lday = (int64_t) day;
3032  jday = (367 * lyear - 7 * (lyear + (lmonth + 9) / 12) / 4
3033  + 275 * lmonth / 9 + lday + 1721014);
3034  /*
3035  * this additional step is only needed if you expect to work on dates
3036  * outside March 1, 1900 to February 28, 2100
3037  */
3038  jday = jday + 15 - 3 * ((lyear + (lmonth - 9) / 7) / 100 + 1) / 4;
3039  return jday;
3040 }
3041 
3042 int jd4713bc_get_date(int64_t jd, int32_t *year, int32_t *month, int32_t *day)
3043 /*******************************************************************
3044 
3045  jd4713bc_get_date
3046 
3047  purpose: get the year, month, day from julian date
3048  (from 4713 BC)
3049  taken from the idl jddate routine
3050 
3051  Returns type: int - no set value now
3052 
3053  Parameters: (in calling order)
3054  Type Name I/O Description
3055  ---- ---- --- -----------
3056  int64_t jd I julian date
3057  int32_t * year O standard 4 digit year
3058  int32_t * month O month
3059  int32_t * day O day of month
3060 
3061  Modification history:
3062  Programmer Date Description of change
3063  ---------- ---- ---------------------
3064  W. Robinson 5 Aug 2013 original development
3065 
3066  *******************************************************************/ {
3067  int64_t v1, v2, v3, v4;
3068  v1 = jd + 68569;
3069  v2 = 4 * v1 / 146097;
3070  v1 = v1 - (146097 * v2 + 3) / 4;
3071  v3 = 4000 * (v1 + 1) / 1461001;
3072  v1 = v1 - 1461 * v3 / 4 + 31;
3073  v4 = 80 * v1 / 2447;
3074  *day = v1 - 2447 * v4 / 80;
3075  v1 = v4 / 11;
3076  *month = v4 + 2 - 12 * v1;
3077  *year = 100 * (v2 - 49) + v3 + v1;
3078  return 0;
3079 }
anc_acq_lin
int anc_acq_lin(int32_t anc_class, l1str *l1rec)
Definition: anc_acq.c:2308
anc_acq_lin_aerosol
int32_t anc_acq_lin_aerosol(l1str *l1rec)
Definition: anc_acq.c:584
out_aerosol
out_aerosol
Definition: anc_acq.c:39
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
ZW
@ ZW
Definition: anc_acq.c:32
HPROF
@ HPROF
Definition: anc_acq.c:36
MIN
#define MIN(x, y)
Definition: rice.h:169
anc_acq_ck
int32_t anc_acq_ck(char *file, char *file_olci)
Definition: anc_acq.c:171
anc_acq_f_stat
int32_t anc_acq_f_stat(char **files, char prioritize_files, int32_t n_anc)
Definition: anc_acq.c:288
BCSCATAU
@ BCSCATAU
Definition: anc_acq.c:41
ANC_STAT_1T
#define ANC_STAT_1T
Definition: anc_acq.c:16
ANC_STAT_CLIM
#define ANC_STAT_CLIM
Definition: anc_acq.c:20
anc_acq_lin_prof
int32_t anc_acq_lin_prof(l1str *l1rec)
Definition: anc_acq.c:738
TOTSCATAU
@ TOTSCATAU
Definition: anc_acq.c:51
met_sto_str_d::e_lon
float e_lon
Definition: anc_acq.c:59
day
int32_t day
Definition: atrem_corl1v2.h:308
MET_UNITS__T_K
#define MET_UNITS__T_K
Definition: met_cvt.h:26
status
int status
Definition: l1_czcs_hdf.c:32
met_sto_str_d
Definition: anc_acq.c:54
FT_OLCI
@ FT_OLCI
Definition: filetype.h:39
SFCT
@ SFCT
Definition: anc_acq.c:32
USE_PMSL
#define USE_PMSL
Definition: anc_acq.c:24
anc_acq_gmao_aer_prep
int32_t anc_acq_gmao_aer_prep(char *file, gen_int_str *aer_int)
Definition: anc_acq.c:1556
anc_acq.h
ANC_SRC_TYP_ECMWF
#define ANC_SRC_TYP_ECMWF
Definition: anc_acq.c:25
TOTEXTTAU
@ TOTEXTTAU
Definition: anc_acq.c:50
ncio_dim_siz
int ncio_dim_siz(int, char *)
Definition: ncio.c:18
out_prof
out_prof
Definition: anc_acq.c:35
qual
int16_t * qual
Definition: l2bin.cpp:86
NULL
#define NULL
Definition: decode_rs.h:63
TPROF
@ TPROF
Definition: anc_acq.c:36
met_sto_str_d::data
float * data[3]
Definition: anc_acq.c:66
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
out_nam
out_nam
Definition: anc_acq.c:31
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
ANC_SRC_TYP_BAD
#define ANC_SRC_TYP_BAD
Definition: anc_acq.c:29
hico.proc_hico.valid_min
valid_min
Definition: proc_hico.py:246
lat
float * lat
Definition: l1_hmodis_hdf.c:756
ANC_SRC_TYP_STD_HDF
#define ANC_SRC_TYP_STD_HDF
Definition: anc_acq.c:26
ON
#define ON
Definition: l1.h:43
ANC_STAT_2T_END
#define ANC_STAT_2T_END
Definition: anc_acq.c:17
prod_ix
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through prod_ix
Definition: HOWTO_Add_a_product.txt:32
jday
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
hico.proc_hico.valid_max
valid_max
Definition: proc_hico.py:247
WV
@ WV
Definition: anc_acq.c:32
met_sto_str_d::lat_step
float lat_step
Definition: anc_acq.c:61
DUEXTTAU
@ DUEXTTAU
Definition: anc_acq.c:42
nlin
int nlin
Definition: get_cmp.c:28
met_sto_str_d::e_lat
float e_lat
Definition: anc_acq.c:62
spacetrack_tles.year
year
Definition: spacetrack_tles.py:170
file
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
DUSCATAU
@ DUSCATAU
Definition: anc_acq.c:43
OCSCATAU
@ OCSCATAU
Definition: anc_acq.c:49
input
instr * input
Definition: main_l2binmatch.cpp:27
anc_acq_gmao_met_prep
int32_t anc_acq_gmao_met_prep(char *file, gen_int_str *met_int)
Definition: anc_acq.c:1139
RHPROF
@ RHPROF
Definition: anc_acq.c:36
l12_proto.h
bilin_interp
float bilin_interp(float *data, int xbox_st, int nx, int ybox_st, float xfrac, float yfrac)
Definition: anc_acq.c:2949
anc_acq_lin_met
int32_t anc_acq_lin_met(l1str *l1rec)
Definition: anc_acq.c:373
ANCBAD
#define ANCBAD
Definition: anc_acq.c:22
met_cvt.h
ICEFR
@ ICEFR
Definition: anc_acq.c:32
NPRM
#define NPRM
Definition: anc_acq.c:21
MET_UNITS__P_PA
#define MET_UNITS__P_PA
Definition: met_cvt.h:24
ANC_STAT_2T_START
#define ANC_STAT_2T_START
Definition: anc_acq.c:18
anc_acq_ecmwf_init
int32_t anc_acq_ecmwf_init(char **files, char **prm_nm, int n_prm, int32_t sto_ix)
Definition: anc_acq.c:2016
SFCRH
@ SFCRH
Definition: anc_acq.c:32
met_sto_str_d::lon_step
float lon_step
Definition: anc_acq.c:56
ncio_grab_stdsclf_ds
int ncio_grab_stdsclf_ds(int, char *, float, float *)
Definition: ncio.c:101
met_sto_str_d::anc_f_stat
int anc_f_stat
Definition: anc_acq.c:65
met_cvt_t_cvt
double met_cvt_t_cvt(double val, int in_type, int out_type)
Definition: met_cvt.c:284
init_anc_add
int init_anc_add(l1str *l1rec)
Definition: anc_acq.c:939
SUSCATAU
@ SUSCATAU
Definition: anc_acq.c:47
PR
@ PR
Definition: anc_acq.c:32
met_sto_str_d::s_lat
float s_lat
Definition: anc_acq.c:60
jd4713bc_get_date
int jd4713bc_get_date(int64_t jd, int32_t *year, int32_t *month, int32_t *day)
Definition: anc_acq.c:3038
jd
Definition: jd.py:1
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
anc_acq_lin_olci
int anc_acq_lin_olci(int anc_class, char *file, l1str *l1rec)
Definition: anc_acq.c:2578
MW
@ MW
Definition: anc_acq.c:32
M_PI
#define M_PI
Definition: pml_iop.h:15
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
met_sto_str_d::nlon
int nlon
Definition: anc_acq.c:57
O3PROF
@ O3PROF
Definition: anc_acq.c:36
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
anc_acq_gmao_prof_prep
int32_t anc_acq_gmao_prof_prep(char *file, gen_int_str *prof_int, int32_t nlvl_expect)
Definition: anc_acq.c:1363
met_cvt_q_to_rh
int met_cvt_q_to_rh(int nval, float *pres, int p_type, float *temp, int t_type, float *q, int q_type, float *rh)
Definition: met_cvt.c:15
OZ_KG_M2_TO_DU
#define OZ_KG_M2_TO_DU
Definition: anc_acq.c:23
SSSCATAU
@ SSSCATAU
Definition: anc_acq.c:45
met_sto_str_d::s_lon
float s_lon
Definition: anc_acq.c:55
MET_UNITS__T_C
#define MET_UNITS__T_C
Definition: met_cvt.h:27
loc
data_t loc[NROOTS]
Definition: decode_rs.h:78
ncio_grab_f_ds
int ncio_grab_f_ds(int, char *, float *)
Definition: ncio.c:57
OCEXTTAU
@ OCEXTTAU
Definition: anc_acq.c:48
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
anc_acq_fnd_t_interp
int32_t anc_acq_fnd_t_interp(double s_time, double *anc_time, int32_t anc_f_stat, int32_t *t_interp, int32_t *data_ix, float *wt)
Definition: anc_acq.c:1711
met_cvt_ttd_to_rh
int met_cvt_ttd_to_rh(int nval, float *temp, int t_type, float *dwp_temp, int dwp_type, float *rh)
Definition: met_cvt.c:182
upcase
char * upcase(char *instr)
Definition: upcase.c:10
v2
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame the required RAM for each execution is MB on the DEC ALPHA and MB on the SGI Octane v2
Definition: HISTORY.txt:728
dpix
int32 dpix
Definition: l1_czcs_hdf.c:22
anc_acq_read_gmao
int32_t anc_acq_read_gmao(char *file, char *group, char *ds_name, float **data, unsigned char **qa, double *time, int32_t *nlon, int32_t *nlat, int32_t *nlvl, double **lon_coord, double **lat_coord)
Definition: anc_acq.c:1826
time
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start time
Definition: HISTORY.txt:248
spix
int32 spix
Definition: l1_czcs_hdf.c:21
iscan
int32_t iscan
Definition: l1_hmodis_hdf.c:745
init_anc_aerosol
int32_t init_anc_aerosol(l1str *l1rec)
Definition: anc_acq.c:898
SFCP
@ SFCP
Definition: anc_acq.c:32
met_sto_str_d::nlat
int nlat
Definition: anc_acq.c:58
files
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT files
Definition: HISTORY.txt:442
SSEXTTAU
@ SSEXTTAU
Definition: anc_acq.c:44
anc_acq_lin_oz
int32_t anc_acq_lin_oz(l1str *l1rec)
Definition: anc_acq.c:989
lon
float * lon
Definition: l1_hmodis_hdf.c:755
epix
int32 epix
Definition: l1_czcs_hdf.c:23
MET_UNITS__P_HPA
#define MET_UNITS__P_HPA
Definition: met_cvt.h:25
anc_acq_gmao_oz_prep
int32_t anc_acq_gmao_oz_prep(char *file, gen_int_str *oz_int)
Definition: anc_acq.c:1484
QPROF
@ QPROF
Definition: anc_acq.c:36
met_cvt_p_cvt
double met_cvt_p_cvt(double val, int in_type, int out_type)
Definition: met_cvt.c:236
RH
@ RH
Definition: anc_acq.c:32
SUEXTTAU
@ SUEXTTAU
Definition: anc_acq.c:46
anc_acq_init
int anc_acq_init(instr *input, l1str *l1rec, int32_t *anc_id)
Definition: anc_acq.c:74
MET_UNITS__Q_KG_KG
#define MET_UNITS__Q_KG_KG
Definition: met_cvt.h:28
val
msiBandIdx val
Definition: l1c_msi.cpp:34
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
met_sto_str_d::data_time
double data_time[3]
Definition: anc_acq.c:63
npix
int npix
Definition: get_cmp.c:27
ANC_SRC_TYP_OLCI
#define ANC_SRC_TYP_OLCI
Definition: anc_acq.c:27
anc_miss_fill
float anc_miss_fill(int32_t prod_ix)
Definition: anc_acq.c:2911
ANC_STAT_3T
#define ANC_STAT_3T
Definition: anc_acq.c:19
jd4713bc_get_jd
int64_t jd4713bc_get_jd(int32_t year, int32_t month, int32_t day)
Definition: anc_acq.c:2999
anc_acq_eval_pt
int32_t anc_acq_eval_pt(gen_int_str *met_int, int32_t iprm, int32_t ilvl, float lat, float lon, int32_t t_interp, int32_t *data_ix, float wt_t1, int32_t ntim_int, int32_t nlvl, int32_t nprm, float *final_val, float *unc)
Definition: anc_acq.c:1614
isodate2unix
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61
ANC_SRC_TYP_GMAO
#define ANC_SRC_TYP_GMAO
Definition: anc_acq.c:28
count
int count
Definition: decode_rs.h:79
BCEXTTAU
@ BCEXTTAU
Definition: anc_acq.c:40


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster