OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
setanc.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 #include "anc_acq.h"
3 #include <ancproto.h>
4 
5 /* ============================================================================ */
6 /* no2conc() - retrieve no2 concentration from ancillary file */
7 /* */
8 /* Written By: B. Franz, NASA OBPG, June 2006. */
9 /* */
10 /* ============================================================================ */
11 
12 #define NXNO2 1440
13 #define NYNO2 720
14 #define NXANC 180
15 #define NYANC 90
16 
17 void no2_frac(float lon, float lat, float *no2_frac_200) {
18  static int firstCall = 1;
19  static int nx = NXANC;
20  static int ny = NYANC;
21  static float dx = 360.0 / NXANC;
22  static float dy = 180.0 / NYANC;
23 
24  typedef float map_frac_t[NXANC + 2];
25  static map_frac_t* map_frac;
26 
27  int i, j;
28  float xx, yy;
29  float t, u;
30 
31  float frac;
32 
33  *no2_frac_200 = 0.0;
34 
35  if (firstCall) {
36 
37  int32 sd_id;
38  int32 sds_id;
39  int32 status;
40  int32 sds_index;
41  int32 rank;
42  int32 nt;
43  int32 dims[H4_MAX_VAR_DIMS];
44  int32 nattrs;
45  int32 start[2];
46  int32 edges[2];
47  char name[H4_MAX_NC_NAME];
48  char sdsname[H4_MAX_NC_NAME];
49  float **map;
50  char *no2_frac_fil, no2_frac_file[300];
51 
52  firstCall = 0;
53 
54  // allocate data
55  map_frac = (map_frac_t*) allocateMemory((NYANC + 2) * sizeof (map_frac_t), "map_frac");
56 
57  if ((no2_frac_fil = getenv("OCDATAROOT")) == NULL) {
58  printf("-E- %s: Error looking up environmental variable OCDATAROOT\n",
59  __FILE__);
60  exit(1);
61  }
62  strcpy(no2_frac_file, no2_frac_fil);
63  strcat(no2_frac_file, "/common/trop_f_no2_200m.hdf");
64 
65  map = (float **) allocate2d_float(NYANC, NXANC);
66  if (map == NULL) {
67  printf("-E- %s: Error allocating NO2 frac space for %s.\n", __FILE__,
68  no2_frac_file);
69  exit(1);
70  }
71 
72  sd_id = SDstart(no2_frac_file, DFACC_RDONLY);
73  if (sd_id == -1) {
74  printf("-E- %s: Error opening NO2 frac file %s.\n", __FILE__,
75  no2_frac_file);
76  exit(1);
77  }
78 
79  printf("\nOpening NO2 frac file %s\n\n", no2_frac_file);
80 
81  strcpy(sdsname, "f_no2_200m");
82  sds_index = SDnametoindex(sd_id, sdsname);
83  if (sds_index == -1) {
84  printf("-E- %s: Error seeking %s SDS from %s.\n", __FILE__, sdsname,
85  no2_frac_file);
86  exit(1);
87  }
88  sds_id = SDselect(sd_id, sds_index);
89 
90  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
91  if (status != 0) {
92  printf("-E- %s: Error reading SDS info for %s from %s.\n", __FILE__,
93  sdsname, no2_frac_file);
94  exit(1);
95  }
96  if (dims[0] != ny || dims[1] != nx) {
97  printf("-E- %s: Dimension mis-match on %s array from %s.\n", __FILE__,
98  sdsname, no2_frac_file);
99  printf(" Expecting %d x %d\n", nx, ny);
100  printf(" Reading %d x %d\n", dims[1], dims[0]);
101  exit(1);
102  }
103 
104  start[0] = 0;
105  start[1] = 0;
106  edges[0] = ny;
107  edges[1] = nx;
108 
109  status = SDreaddata(sds_id, start, NULL, edges, (VOIDP) & map[0][0]);
110  if (status != 0) {
111  printf("-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
112  no2_frac_file);
113  exit(1);
114  }
115 
116  for (j = 0; j < ny; j++) {
117  for (i = 0; i < nx; i++) {
118  map_frac[j + 1][i + 1] = map[ny - j - 1][i];
119  }
120  }
121 
122  /* add boarders to simplify interpolation */
123 
124  for (j = 0; j < ny; j++) {
125  map_frac[j + 1][0] = map_frac[j + 1][nx];
126  map_frac[j + 1][nx + 1] = map_frac[j + 1][1];
127  }
128  for (i = 0; i < nx + 2; i++) {
129  map_frac[0][i] = map_frac[1][i];
130  map_frac[ny + 1][i] = map_frac[ny][i];
131  }
132 
133  SDendaccess(sds_id);
134  SDend(sd_id);
135  free2d_float(map);
136  }
137 
138  /* interpolate to pixel location */
139 
140  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
141  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
142 
143  xx = i * dx - 180.0 - dx / 2;
144  yy = j * dy - 90.0 - dy / 2;
145 
146  t = (lon - xx) / dx;
147  u = (lat - yy) / dy;
148 
149  frac = (1 - t) * (1 - u) * map_frac[j][i] + t * (1 - u) * map_frac[j][i + 1]
150  + t * u * map_frac[j + 1][i + 1] + (1 - t) * u * map_frac[j + 1][i];
151 
152  /* return components of stratospheric and tropospheric no2 */
153 
154  *no2_frac_200 = MAX(frac, 0.0);
155 
156  return;
157 }
158 
159 void no2conc(char *no2file, float lon, float lat, int32_t doy, float *no2_tropo,
160  float *no2_strat) {
161  static int firstCall = 1;
162  static int nx = NXNO2;
163  static int ny = NYNO2;
164  static float dx = 360.0 / NXNO2;
165  static float dy = 180.0 / NYNO2;
166 
167  typedef float map_t[NXNO2 + 2];
168  static map_t* map_total;
169  static map_t* map_tropo;
170 
171  int i, j;
172  float xx, yy;
173  float t, u;
174 
175  float total;
176  float tropo;
177 
178  *no2_tropo = 0.0;
179  *no2_strat = 0.0;
180 
181  if (firstCall) {
182 
183  int32 sd_id;
184  int32 sds_id;
185  int32 status;
186  int32 sds_index;
187  int32 rank;
188  int32 nt;
189  int32 dims[H4_MAX_VAR_DIMS];
190  int32 nattrs;
191  int32 start[2];
192  int32 edges[2];
193  char name[H4_MAX_NC_NAME];
194  char sdsname[H4_MAX_NC_NAME];
195  float **map;
196  char title[255];
197  int16 month;
198  char mstr[4] = "";
199 
200  firstCall = 0;
201 
202  // allocate data
203  map_total = (map_t*) allocateMemory((NYNO2 + 2) * sizeof (map_t), "map_total");
204  map_tropo = (map_t*) allocateMemory((NYNO2 + 2) * sizeof (map_t), "map_tropo");
205 
206  map = (float **) allocate2d_float(NYNO2, NXNO2);
207  if (map == NULL) {
208  printf("-E- %s: Error allocating space for %s.\n", __FILE__, no2file);
209  exit(1);
210  }
211 
212  sd_id = SDstart(no2file, DFACC_RDONLY);
213  if (sd_id == -1) {
214  printf("-E- %s: Error opening NO2 file %s.\n", __FILE__, no2file);
215  exit(1);
216  }
217 
218  printf("\nOpening NO2 file %s\n\n", no2file);
219 
220  if (SDreadattr(sd_id, SDfindattr(sd_id, "Title"), (VOIDP) title) == 0) {
221  if (strstr(title, "NO2 Climatology") != NULL) {
222  month = (int) doy / 31;
223  sprintf(mstr, "_%2.2i", month + 1);
224  }
225  }
226 
227  strcpy(sdsname, "tot_no2");
228  strcat(sdsname, mstr);
229  sds_index = SDnametoindex(sd_id, sdsname);
230  if (sds_index == -1) {
231  printf("-E- %s: Error seeking %s SDS from %s.\n", __FILE__, sdsname,
232  no2file);
233  exit(1);
234  }
235  sds_id = SDselect(sd_id, sds_index);
236 
237  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
238  if (status != 0) {
239  printf("-E- %s: Error reading SDS info for %s from %s.\n", __FILE__,
240  sdsname, no2file);
241  exit(1);
242  }
243  if (dims[0] != ny || dims[1] != nx) {
244  printf("-E- %s: Dimension mis-match on %s array from %s.\n", __FILE__,
245  sdsname, no2file);
246  printf(" Expecting %d x %d\n", nx, ny);
247  printf(" Reading %d x %d\n", dims[1], dims[0]);
248  exit(1);
249  }
250 
251  start[0] = 0;
252  start[1] = 0;
253  edges[0] = ny;
254  edges[1] = nx;
255 
256  status = SDreaddata(sds_id, start, NULL, edges, (VOIDP) & map[0][0]);
257  if (status != 0) {
258  printf("-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
259  no2file);
260  exit(1);
261  }
262 
263  for (j = 0; j < ny; j++) {
264  for (i = 0; i < nx; i++) {
265  map_total[j + 1][i + 1] = map[ny - j - 1][i];
266  }
267  }
268 
269  status = SDendaccess(sds_id);
270  if (status != 0) {
271  printf("-E- %s: Error closing SDS %s from %s.\n", __FILE__, sdsname,
272  no2file);
273  exit(1);
274  }
275 
276  strcpy(sdsname, "trop_no2");
277  strcat(sdsname, mstr);
278  sds_index = SDnametoindex(sd_id, sdsname);
279  if (sds_index == -1) {
280  printf("-E- %s: Error seeking %s SDS from %s.\n", __FILE__, sdsname,
281  no2file);
282  exit(1);
283  }
284  sds_id = SDselect(sd_id, sds_index);
285 
286  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
287  if (status != 0) {
288  printf("-E- %s: Error reading SDS info for %s from %s.\n", __FILE__,
289  sdsname, no2file);
290  exit(1);
291  }
292  if (dims[0] != ny || dims[1] != nx) {
293  printf("-E- %s: Dimension mis-match on %s array from %s.\n", __FILE__,
294  sdsname, no2file);
295  printf(" Expecting %d x %d\n", nx, ny);
296  printf(" Reading %d x %d\n", dims[1], dims[0]);
297  exit(1);
298  }
299 
300  start[0] = 0;
301  start[1] = 0;
302  edges[0] = ny;
303  edges[1] = nx;
304 
305  status = SDreaddata(sds_id, start, NULL, edges, (VOIDP) & map[0][0]);
306  if (status != 0) {
307  printf("-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
308  no2file);
309  exit(1);
310  }
311 
312  status = SDendaccess(sds_id);
313  if (status != 0) {
314  printf("-E- %s: Error closing SDS %s from %s.\n", __FILE__, sdsname,
315  no2file);
316  exit(1);
317  }
318 
319  for (j = 0; j < ny; j++) {
320  for (i = 0; i < nx; i++) {
321  map_tropo[j + 1][i + 1] = map[ny - j - 1][i];
322  }
323  }
324 
325  /* add boarders to simplify interpolation */
326 
327  for (j = 0; j < ny; j++) {
328  map_total[j + 1][0] = map_total[j + 1][nx];
329  map_total[j + 1][nx + 1] = map_total[j + 1][1];
330  map_tropo[j + 1][0] = map_tropo[j + 1][nx];
331  map_tropo[j + 1][nx + 1] = map_tropo[j + 1][1];
332  }
333  for (i = 0; i < nx + 2; i++) {
334  map_total[0][i] = map_total[1][i];
335  map_total[ny + 1][i] = map_total[ny][i];
336  map_tropo[0][i] = map_tropo[1][i];
337  map_tropo[ny + 1][i] = map_tropo[ny][i];
338  }
339 
340  free2d_float(map);
341  SDend(sd_id);
342  }
343 
344  /* interpolate to pixel location */
345 
346  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
347  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
348 
349  xx = i * dx - 180.0 - dx / 2;
350  yy = j * dy - 90.0 - dy / 2;
351 
352  t = (lon - xx) / dx;
353  u = (lat - yy) / dy;
354 
355  total = (1 - t) * (1 - u) * map_total[j][i]
356  + t * (1 - u) * map_total[j][i + 1] + t * u * map_total[j + 1][i + 1]
357  + (1 - t) * u * map_total[j + 1][i];
358 
359  tropo = (1 - t) * (1 - u) * map_tropo[j][i]
360  + t * (1 - u) * map_tropo[j][i + 1] + t * u * map_tropo[j + 1][i + 1]
361  + (1 - t) * u * map_tropo[j + 1][i];
362 
363  /* return components of stratospheric and tropospheric no2 */
364 
365  *no2_strat = MAX(total - tropo, 0.0);
366  *no2_tropo = MAX(tropo, 0.0);
367 
368  return;
369 }
370 
371 /* ============================================================================ */
372 /* ozone_climatology() - retrieve ozone concentration from daily climatology */
373 /* */
374 /* Written By: B. Franz, NASA OBPG, April 2009. */
375 /* W. Robinson, SAIC, 14 Feb 2014 generalize code for varying grid size */
376 /* */
377 
378 /* ============================================================================ */
379 float ozone_climatology(char *file, int day, float lon, float lat) {
380  static int firstCall = 1;
381  static int nx, ny, mapx;
382  static float dx, dy, *map;
383 
384  int i, j;
385  float xx, yy, *tmp;
386  float t, u, val_11, val_12, val_21, val_22;
387  float ozone;
388 
389  if (firstCall) {
390 
391  int32 sd_id;
392  int32 sds_id;
393  int32 status;
394  int32 sds_index;
395  int32 rank;
396  int32 nt;
397  int32 dims[H4_MAX_VAR_DIMS];
398  int32 nattrs;
399  int32 start[2];
400  int32 edges[2];
401  char name[H4_MAX_NC_NAME];
402  char sdsname[H4_MAX_NC_NAME];
403 
404  firstCall = 0;
405 
406  if (day < 1 || day > 366) {
407  printf("-E- %s: Bogus day number for ozone look-up: %d\n", __FILE__,
408  day);
409  exit(1);
410  }
411 
412  /*
413  if ((path = getenv("OCDATAROOT")) == NULL)
414  {
415  printf("-E- %s: Error looking up environmental variable OCDATAROOT\n",
416  __FILE__);
417  exit(1);
418  }
419  strcpy(file, path);
420  strcat(file, "/common/ozone_climatology.hdf");
421  */
422 
423  sd_id = SDstart(file, DFACC_RDONLY);
424  if (sd_id == -1) {
425  printf("-E- %s: Error openin file %s.\n", __FILE__, file);
426  exit(1);
427  }
428 
429  printf("\nOpening ozone file %s\n\n", file);
430 
431  sprintf(sdsname, "ozone_mean_%03d", day);
432  sds_index = SDnametoindex(sd_id, sdsname);
433  if (sds_index == -1) {
434  printf("-E- %s: Error seeking %s SDS from %s.\n", __FILE__, sdsname,
435  file);
436  exit(1);
437  }
438  sds_id = SDselect(sd_id, sds_index);
439 
440  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
441  if (status != 0) {
442  printf("-E- %s: Error reading SDS info for %s from %s.\n", __FILE__,
443  sdsname, file);
444  exit(1);
445  }
446  ny = dims[0];
447  nx = dims[1];
448  mapx = nx + 2;
449  if (nx < 0 || ny < 0) {
450  printf("-E- %s: grid has bad dimensions, sds %s from %s.\n", __FILE__,
451  sdsname, file);
452  printf(" Reading %d x %d\n", nx, ny);
453  exit(1);
454  }
455 
456  dx = 360. / nx;
457  dy = 180. / ny;
458  /*
459  * allocate the storage for initial read and final interplation grid
460  */
461  if (((tmp = (float *) malloc(nx * ny * sizeof (float))) == NULL) ||
462  ((map = (float *) malloc(mapx * (ny + 2) * sizeof (float)))
463  == NULL)) {
464  printf(
465  "-E- %s, %d: Unable to allocate space for climatology grid storage\n",
466  __FILE__, __LINE__);
467  exit(1);
468  }
469 
470  start[0] = 0;
471  start[1] = 0;
472  edges[0] = ny;
473  edges[1] = nx;
474 
475  status = SDreaddata(sds_id, start, NULL, edges, (VOIDP) tmp);
476  if (status != 0) {
477  printf("-E- %s: Error reading SDS %s from %s.\n", __FILE__, sdsname,
478  file);
479  exit(1);
480  }
481  /* upright the latitude during transfer ( ny - j - 1 part) */
482  for (j = 0; j < ny; j++) {
483  for (i = 0; i < nx; i++) {
484  *(map + i + 1 + mapx * (j + 1)) = *(tmp + i + nx * j);
485  *(map + i + 1 + mapx * (j + 1)) =
486  *(tmp + i + nx * (ny - j - 1));
487  }
488  }
489 
490  /* add boarders to simplify interpolation - on sides (lon),
491  duplicate the column from the other side, and at top, bottom (N, S),
492  duplicate nearest neighbor */
493 
494  for (j = 0; j < ny; j++) {
495  *(map + mapx * (j + 1)) = *(map + nx + mapx * (j + 1));
496  *(map + (nx + 1) + mapx * (j + 1)) =
497  *(map + 1 + mapx * (j + 1));
498  }
499  for (i = 0; i < mapx; i++) {
500  *(map + i) = *(map + mapx + i);
501  *(map + i + mapx * (ny + 1)) = *(map + i + mapx * ny);
502  }
503 
504  SDendaccess(sds_id);
505  SDend(sd_id);
506  free(tmp);
507  }
508 
509  /* interpolate to pixel location */
510 
511  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
512  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
513 
514  xx = i * dx - 180.0 - dx / 2;
515  yy = j * dy - 90.0 - dy / 2;
516 
517  t = (lon - xx) / dx;
518  u = (lat - yy) / dy;
519 
520  val_11 = *(map + i + mapx * j); /* [i,j] */
521  val_12 = *(map + i + 1 + mapx * j); /* [i+1,j] */
522  val_21 = *(map + i + mapx * (j + 1)); /* [i,j+1] */
523  val_22 = *(map + i + 1 + mapx * (j + 1)); /* [i+1,j+1] */
524 
525  if ((val_11 < 0) || (val_12 < 0) || (val_21 < 0) || (val_22 < 0)) {
526  ozone = BAD_FLT;
527  printf("I %s, %d: Attempt to use missing climatology points\n",
528  __FILE__, __LINE__);
529  } else
530  ozone = (1 - t) * (1 - u) * val_11 + t * (1 - u) * val_12
531  + t * u * val_22 + (1 - t) * u * val_21;
532  return (ozone);
533 }
534 
535 /* -----------------------------------------------------------------------
536  function setanc
537 
538  Returns windspeed, pressure, rel. humidity, water vapor, and ozone
539  concentration for a specified time and place. Uses climatology
540  or realtime ancillary data files in SeaWiFS hdf format. Ancillary
541  file names are retrieved from the ancfiles.dat file.
542 
543  Returns 1 on error.
544 
545  Inputs:
546 
547  Outputs:
548  zw(npix) Zonal Wind m/s
549  mw(npix) Meridional Wind m/s
550  ws(npix) Wind Speed m/s
551  pr(npix) Surface Pressure millibars
552  rh(npix) Relative Humidity %
553  pw(npix) Water Vapor g/cm^2
554  oz(npix) Ozone atm-cm
555  no2(npix) NO2 molecules cm^-2
556 
557 
558  Written By: BA Franz, GSC, 6/97
559  Conversion to C: G Fu, GSC, 3/99
560 
561  ----------------------------------------------------------------------- */
562 
563 int setanc(l1str *l1rec) {
564  static float r2d = 57.29577951;
565  static int firstCall = 1;
566  static int32_t anc_id[] ={
567  -1, -1
568  };
569  static short *ancqc = NULL;
570  static char *no2file = NULL;
571  static short npix; // need this for l3gen which changes the size of a line
572 
573  float u, v, u_u, v_u, ws_2;
574 
575  int32_t i, retval, status;
576  float *lat = (float *) l1rec->lat;
577  float *lon = (float *) l1rec->lon;
578  short year, jday;
579  double dsec;
580  unix2yds(l1rec->scantime, &year, &jday, &dsec);
581  int32_t msec = (int32_t) (dsec * 1.e3);
582 
583  short parmID;
584 
585  status = 0;
586 
587  if (firstCall) {
588  firstCall = 0;
589  npix = (short) l1rec->npix;
590  printf("\nOpening meteorological files.\n");
591  printf(" met1 = %s\n", input->met1);
592  printf(" met2 = %s\n", input->met2);
593  printf(" met3 = %s\n", input->met3);
594  printf(" ozone1 = %s\n", input->ozone1);
595  printf(" ozone2 = %s\n", input->ozone2);
596  printf(" ozone3 = %s\n", input->ozone3);
597  printf(" no2 = %s\n", input->no2file);
598  printf("\n");
599  if ((ancqc = calloc(npix, sizeof (short))) == NULL) {
600  fprintf(stderr, "-E- %s %d: Unable to allocate buffer space.\n", __FILE__,
601  __LINE__);
602  exit(1);
603  }
604  if (strcmp(input->no2file, "") != 0) {
605  no2file = input->no2file;
606  }
607  /*
608  * do the setup and identification of alternate anc files
609  * (currently, only the ECMWF source)
610  */
611  if (anc_acq_init(input, l1rec, anc_id) != 0)
612  return 1;
613  } else if (l1rec->npix > npix) {
614  npix = (short) l1rec->npix;
615  free(ancqc);
616  if ((ancqc = calloc(npix, sizeof (short))) == NULL) {
617  fprintf(stderr, "-E- %s %d: Unable to allocate buffer space.\n", __FILE__,
618  __LINE__);
619  exit(1);
620  }
621  }
622 
623 
624  /*
625  * get standard met if anc_id = 1 or get met from ECMWF if anc_id = 0
626  * and OLCI tie point data if anc_id = 2
627  */
628  if (anc_id[0] == 0) {
629  if (anc_acq_lin(0, l1rec) != 0)
630  return 1;
631  } else if (anc_id[0] == 2) {
632  if (anc_acq_lin_olci(0, input->met1, l1rec) != 0)
633  return 1;
634  } else if (anc_id[0] == 3) {
635  if (anc_acq_lin_met(l1rec) != 0)
636  return 1;
637  } else {
638  /* relative humidity */
639  /* ----------------- */
640 
641  parmID = 5;
642  if(input->relhumid != -2000) {
643  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec, input->met1,
644  input->met2, input->met3, input->anc_cor_file, parmID, l1rec->rh,
645  l1rec->rh_unc, l1rec->ancqc);
646  if (retval != 0) {
647  fprintf(stderr,
648  "-E- %s %d: Error loading relative humidity ancillary data. %s\n",
649  __FILE__, __LINE__, input->anc_cor_file);
650  status = 1;
651  }
652  }
653 
654  /* wind speed */
655  /* ---------- */
656 
657  parmID = 0;
658  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec, input->met1,
659  input->met2, input->met3, input->anc_cor_file, parmID, l1rec->zw,
660  l1rec->zw_unc, ancqc);
661  if (retval != 0) {
662  fprintf(stderr,
663  "-E- %s %d: Error loading Zonal wind speed ancillary data.\n", __FILE__,
664  __LINE__);
665  status = 1;
666  }
667 
668  parmID = 1;
669  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec, input->met1,
670  input->met2, input->met3, input->anc_cor_file, parmID, l1rec->mw,
671  l1rec->mw_unc, ancqc);
672  if (retval != 0) {
673  fprintf(stderr,
674  "-E- %s %d: Error loading Meridional wind speed ancillary data.\n",
675  __FILE__, __LINE__);
676  status = 1;
677  }
678  /* create wind speed, direction from u, v components */
679  for (i = 0; i < npix; i++) {
680  u = l1rec->zw[i];
681  v = l1rec->mw[i];
682  u_u = l1rec->zw_unc[i];
683  v_u = l1rec->mw_unc[i];
684  ws_2 = u * u + v * v;
685  /* */
686  if(input->windspeed != -2000)
687  l1rec->ws[i] = sqrt(ws_2);
688  if(input->windangle != -2000)
689  l1rec->wd[i] = atan2f(-l1rec->zw[i], -l1rec->mw[i]) * r2d;
690  /*
691  * make uncertainties in the speed, direction too
692  * when u+v real small use alternate formulation
693  */
694  u = fabs(u);
695  v = fabs(v);
696  if ((u + v) > 0.05 * (u_u + v_u)) {
697  l1rec->ws_unc[i] = sqrt((u * u * u_u * u_u + v * v * v_u * v_u) /
698  ws_2);
699  l1rec->wd_unc[i] = sqrt(v * v * u_u * u_u + u * u * v_u * v_u) / ws_2;
700  if (l1rec->wd_unc[i] > M_PI) l1rec->wd_unc[i] = M_PI;
701  } else {
702  l1rec->ws_unc[i] = sqrt(0.5 * (u_u * u_u + v_u * v_u));
703  l1rec->wd_unc[i] = M_PI;
704  }
705  l1rec->wd_unc[i] *= r2d;
706  l1rec->ancqc[i] |= ancqc[i];
707  }
708 
709  /* surface pressure */
710  /* ---------------- */
711 
712  parmID = 2;
713  if(input->pressure != -2000) {
714  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec, input->met1,
715  input->met2, input->met3, input->anc_cor_file, parmID, l1rec->pr,
716  l1rec->pr_unc, ancqc);
717  if (retval != 0) {
718  fprintf(stderr,
719  "-E- %s %d: Error loading surface pressure ancillary data.\n", __FILE__,
720  __LINE__);
721  status = 1;
722  }
723  }
724 
725  for (i = 0; i < npix; i++) {
726 
727  if (l1rec->pr[i] <= 0.0 || isnan(l1rec->pr[i]))
728  l1rec->pr[i] = 1013.25;
729  else if (l1rec->pr[i] < 900.0)
730  l1rec->pr[i] = 900.0;
731  else if (l1rec->pr[i] > 1100.0)
732  l1rec->pr[i] = 1100.0;
733 
734  l1rec->ancqc[i] |= ancqc[i];
735 
736  /* if processing land, adjust pressure for terrain height */
737  if (input->proc_land && l1rec->height[i] != 0.0) {
738  l1rec->pr[i] *= exp(-l1rec->height[i] / 8434);
739  }
740 
741  }
742 
743  /* precipitable water (water vapor) */
744  /* -------------------------------- */
745 
746  parmID = 3;
747  if(input->watervapor != -2000) {
748  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec, input->met1,
749  input->met2, input->met3, input->anc_cor_file, parmID, l1rec->wv,
750  l1rec->wv_unc, ancqc);
751  if (retval != 0) {
752  fprintf(stderr,
753  "-E- %s %d: Error loading precipitable water ancillary data.\n",
754  __FILE__, __LINE__);
755  status = 1;
756  }
757 
758  /* convert from kg/m^2 to g/cm^2 */
759  for (i = 0; i < npix; i++) {
760  l1rec->wv[i] = l1rec->wv[i] / 10.0;
761  l1rec->wv_unc[i] = l1rec->wv_unc[i] / 10.0;
762  l1rec->ancqc[i] |= ancqc[i];
763  }
764  }
765  } /* end of met portion ingest */
766 
767  /* here, get the anc_profile. Note that since only GMAO ancillary is
768  supported, no need to identify WHAT kind of profile it is in anc_acq_ck.
769  */
770  if (strlen(input->anc_profile1)){
771  if( anc_acq_lin_prof( l1rec ) != 0 )
772  return 1;
773  }
774  /* here, get the anc_aerosol. Note that since only GMAO ancillary is
775  supported, no need to identify WHAT kind of profile it is in anc_acq_ck.
776  */
777  if (strlen(input->anc_aerosol1)){
778  if( anc_acq_lin_aerosol( l1rec ) != 0 )
779  return 1;
780  }
781  // for the surface albedo
782  if( strlen(input->sfc_albedo) ){
783  if( acq_sfc_albedo( l1rec ) != 0 )
784  return 1;
785  } else
786  if( input->proc_cloud ) {
787  printf(
788 "%s, %d: Cloud processing, requires an input surface albedo (sfc_albedo=...)\n",
789  __FILE__, __LINE__ );
790  return 1;
791  }
792  /* ozone */
793  /* ----- */
794  if (anc_id[1] == 0) {
795  if (anc_acq_lin(1, l1rec) != 0) /* arg 1 = 1 to acess ozone data */
796  return 1;
797  } else if (anc_id[1] == 2) {
798  if (anc_acq_lin_olci(1, input->met1, l1rec) != 0)
799  return 1;
800  } else if (anc_id[1] == 3) {
801  if (anc_acq_lin_oz(l1rec) != 0)
802  return 1;
803  } else
804  if(input->ozone != -2000) {
805  if (strstr(input->ozone1, "ozone_climatology") != NULL) {
806  for (i = 0; i < npix; i++) {
807  l1rec->oz[i] = ozone_climatology(input->ozone1, jday, l1rec->lon[i], l1rec->lat[i]);
808  if (l1rec->oz[i] < 0.0)
809  ancqc[i] = 1;
810  else
811  ancqc[i] = 0;
812  }
813  } else {
814  parmID = 4;
815  retval = get_ancillary(lat, lon, npix, year, jday, jday, msec,
816  input->ozone1, input->ozone2, input->ozone3, input->anc_cor_file,
817  parmID, l1rec->oz, l1rec->oz_unc, l1rec->ancqc);
818  if (retval != 0) {
819  fprintf(stderr, "-E- %s %d: Error loading Ozone ancillary data.\n",
820  __FILE__, __LINE__);
821  status = 1;
822  }
823  }
824 
825  /* convert from Dobson units to atm-cm */
826  for (i = 0; i < npix; i++) {
827  l1rec->oz[i] = l1rec->oz[i] / 1000.0;
828  l1rec->oz_unc[i] = l1rec->oz_unc[i] / 1000.0;
829  l1rec->ancqc[i] |= ancqc[i];
830  }
831  }
832 
833  /*
834  * where l1rec->ancqc is set, turn on ATMFAIL in the flags
835  */
836  for (i = 0; i < npix; i++)
837  if (l1rec->ancqc[i] != 0)
838  l1rec->flags[i] |= ATMWARN;
839 
840  /* no2 and fraction */
841  /* ---------------- */
842 
843  if ((input->gas_opt & NO2_BIT) != 0)
844  for (i = 0; i < npix; i++) {
845  no2conc(no2file, l1rec->lon[i], l1rec->lat[i], jday,
846  &l1rec->no2_tropo[i], &l1rec->no2_strat[i]);
847  l1rec->no2_tropo[i] *= 1e15;
848  l1rec->no2_strat[i] *= 1e15;
849  // the no2_tropo_unc and no2_strat_unc will also need this whenever
850  // the unc is available
851  no2_frac(l1rec->lon[i], l1rec->lat[i], &l1rec->no2_frac[i]);
852  }
853 
854  if (status != 0) {
855  fprintf(stderr, "-E- %s %d: Error loading ancillary data.\n", __FILE__,
856  __LINE__);
857  return (1);
858  }
859 
860  return (0);
861 }
862 
int anc_acq_lin(int32_t anc_class, l1str *l1rec)
Definition: anc_acq.c:2308
int32_t anc_acq_lin_aerosol(l1str *l1rec)
Definition: anc_acq.c:584
int setanc(l1str *l1rec)
Definition: setanc.c:563
#define MAX(A, B)
Definition: swl0_utils.h:26
integer, parameter int16
Definition: cubeio.f90:3
#define MIN(x, y)
Definition: rice.h:169
data_t t[NROOTS+1]
Definition: decode_rs.h:77
int32_t anc_acq_lin_prof(l1str *l1rec)
Definition: anc_acq.c:738
int j
Definition: decode_rs.h:73
int32_t day
int status
Definition: l1_czcs_hdf.c:32
void no2conc(char *no2file, float lon, float lat, int32_t doy, float *no2_tropo, float *no2_strat)
Definition: setanc.c:159
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
#define NULL
Definition: decode_rs.h:63
#define NYANC
Definition: setanc.c:15
read l1rec
float ozone_climatology(char *file, int day, float lon, float lat)
Definition: setanc.c:379
float * lat
int32 * msec
Definition: l1_czcs_hdf.c:31
#define NXNO2
Definition: setanc.c:12
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
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
instr * input
int32_t anc_acq_lin_met(l1str *l1rec)
Definition: anc_acq.c:373
data_t tmp
Definition: decode_rs.h:74
#define NXANC
Definition: setanc.c:14
int get_ancillary(float *in_lat, float *in_lon, int16_t cnt, int16_t syear, int16_t sday, int16_t eday, int32_t time, char *filename1, char *filename2, char *filename3, char *anc_cor_file, int16_t parm_flag, float *interp, float *anc_unc, int16_t *qcflag)
Definition: getanc.c:297
MOD_PR02AQUA Production where MYD is the prefix denoting AQUA file output The total
void free2d_float(float **p)
Free a two-dimensional array created by allocate2d_float.
Definition: allocate2d.c:142
#define NO2_BIT
Definition: l12_parms.h:59
int anc_acq_lin_olci(int anc_class, char *file, l1str *l1rec)
Definition: anc_acq.c:2578
#define M_PI
Definition: pml_iop.h:15
void unix2yds(double usec, short *year, short *day, double *secs)
#define ATMWARN
Definition: l2_flags.h:33
#define NYNO2
Definition: setanc.c:13
#define BAD_FLT
Definition: jplaeriallib.h:19
data_t u
Definition: decode_rs.h:74
#define fabs(a)
Definition: misc.h:93
Extra metadata that will be written to the HDF4 file l2prod rank
int32_t anc_acq_lin_oz(l1str *l1rec)
Definition: anc_acq.c:989
float * lon
void no2_frac(float lon, float lat, float *no2_frac_200)
Definition: setanc.c:17
float ** allocate2d_float(size_t h, size_t w)
Allocate a two-dimensional array of type float of a given size.
Definition: allocate2d.c:125
int acq_sfc_albedo(l1str *l1rec)
int anc_acq_init(instr *input, l1str *l1rec, int32_t *anc_id)
Definition: anc_acq.c:74
int i
Definition: decode_rs.h:71
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int npix
Definition: get_cmp.c:27