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


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster