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


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