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
sstref.c
Go to the documentation of this file.
1 #include <sys/types.h>
2 #include <unistd.h>
3 
4 /* ============================================================================ */
5 /* module sstref.c - retrieve sst reference temperature from ancillary file */
6 /* */
7 /* Written By: B. Franz, NASA/SIMBIOS, August 2003. */
8 /* */
9 /* ============================================================================ */
10 
11 /*
12  * hardcoded the minatsrcnt value to 1 so I could eliminate the input parameter.
13  * Only used for get_atsr anyway, and thus only valid for 2006-2009...
14  */
15 #include "l12_proto.h"
16 #include <netcdf.h>
17 #include "l12_parms.h"
18 #include "l1_aci.h"
19 #define PATHCLIM 1
20 #define OISSTBIN 2
21 #define OISSTV2D 3
22 #define AMSRE3DAY 4
23 #define AMSREDAY 5
24 #define ATSR 6
25 #define NTEV2 7
26 #define AMSRE3DN 8
27 #define ATSRDAY 9
28 #define WINDSAT3DAY 10
29 #define WINDSATDAY 11
30 #define WINDSAT3DN 12
31 #define CMC 13
32 #define L3M 14
33 
34 #define CtoK 273.15 /* conversion between degree C and Kelvin*/
35 static float sstbad = BAD_FLT;
36 static int32_t format = -1;
37 
38 static int MiddleOfMonth[2][12] = {
39  {15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349},
40  {15, 45, 75, 106, 136, 167, 197, 228, 259, 289, 320, 350}};
41 
42 
43 /* ----------------------------------------------------------------------------------- */
44 /* get_oisstv2d() - read and interpolate Reynolds 0.25-deg daily V2 netcdf OI files */
45 /* */
46 /* B. Franz, SAIC, July 2008. */
47 /* ----------------------------------------------------------------------------------- */
48 
49 #define OI4NX 1440
50 #define OI4NY 720
51 
52 float get_oisstv2d(char *sstfile, float lon, float lat) {
53  static int firstCall = 1;
54  static int nx = OI4NX;
55  static int ny = OI4NY;
56  static float dx = 360.0 / OI4NX;
57  static float dy = 180.0 / OI4NY;
58 
59  typedef float sstref_t[OI4NX + 2];
60  static sstref_t* sstref;
61  //static float sstref[OI4NY+2][OI4NX+2];
62 
63  float sst = sstbad;
64  int i, j, ii;
65  int ntmp;
66  float xx, yy;
67  float t, u;
68  float reftmp[2][2];
69  float ftmp;
70 
71  if (firstCall) {
72  sstref = (sstref_t*) malloc((OI4NY + 2) * sizeof (sstref_t));
73 
74  typedef int16 ssttmp_t[OI4NX];
75  ssttmp_t *ssttmp;
76  ssttmp = (ssttmp_t*) malloc(OI4NY * sizeof (ssttmp_t));
77 
78  char name [H4_MAX_NC_NAME] = "";
79  char sdsname[H4_MAX_NC_NAME] = "";
80  int ncid, ndims, nvars, ngatts, unlimdimid;
81  int32 sd_id;
82  int32 sds_id;
83  int32 rank;
84  int32 nt;
85  int32 dims[H4_MAX_VAR_DIMS];
86  int32 nattrs;
87  int32 start[4] = {0, 0, 0, 0};
88  int32 status;
89  float slope;
90  float offset;
91 
92  firstCall = 0;
93 
94  printf("Loading Daily V2 0.25-deg OI Reynolds SST reference from %s\n", sstfile);
95  printf("\n");
96 
97  /* Open the file */
98  strcpy(sdsname, "sst");
99 
100  /* Try HDF first... */
101 
102  // SDstart seg faults on the mac, so need to protect the
103  // call with Hishdf()
104  if (Hishdf(sstfile))
105  sd_id = SDstart(sstfile, DFACC_RDONLY);
106  else
107  sd_id = FAIL;
108 
109  if (sd_id != FAIL) {
110 
111  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
112  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
113  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) & ssttmp[0][0]);
114  if (status != 0) {
115  fprintf(stderr, "-E- %s Line %d: Error reading SDS %s from %s.\n",
116  __FILE__, __LINE__, sdsname, sstfile);
117  exit(1);
118  }
119  if (getHDFattr(sd_id, "scale_factor", sdsname, (VOIDP) & slope) != 0) {
120  fprintf(stderr, "-E- %s line %d: error reading scale factor.\n",
121  __FILE__, __LINE__);
122  exit(1);
123  }
124  if (getHDFattr(sd_id, "add_offset", sdsname, (VOIDP) & offset) != 0) {
125  fprintf(stderr, "-E- %s line %d: error reading scale offset.\n",
126  __FILE__, __LINE__);
127  exit(1);
128  }
129  status = SDendaccess(sds_id);
130  status = SDend(sd_id);
131  } else {
132 
133  /* try netCDF */
134  if (nc_open(sstfile, NC_NOWRITE, &ncid) == NC_NOERR) {
135 
136  status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
137  status = nc_inq_varid(ncid, sdsname, &sds_id);
138 
139  /* Read the data. */
140  if (nc_get_var(ncid, sds_id, &ssttmp[0][0]) != NC_NOERR) {
141  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
142  __FILE__, __LINE__, sdsname, sstfile);
143  exit(1);
144  }
145 
146  if (nc_get_att_float(ncid, sds_id, "scale_factor", &slope) != NC_NOERR) {
147  fprintf(stderr, "-E- %s line %d: error reading scale factor.\n",
148  __FILE__, __LINE__);
149  exit(1);
150  }
151 
152  if (nc_get_att_float(ncid, sds_id, "add_offset", &offset) != NC_NOERR) {
153  fprintf(stderr, "-E- %s line %d: error reading scale offset.\n",
154  __FILE__, __LINE__);
155  exit(1);
156  }
157  /* Close the file */
158  if (nc_close(ncid) != NC_NOERR) {
159  fprintf(stderr, "-E- %s line %d: error closing %s.\n",
160  __FILE__, __LINE__, sstfile);
161  exit(1);
162  }
163 
164  } else {
165  fprintf(stderr, "-E- %s line %d: error reading %s.\n",
166  __FILE__, __LINE__, sstfile);
167  exit(1);
168  }
169  }
170 
171  /* rotate 180-deg and add wrapping border to simplify interpolation */
172  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
173 
174  for (j = 0; j < ny; j++) {
175  for (i = 0; i < nx; i++) {
176  ii = (i < nx / 2) ? i + nx / 2 : i - nx / 2;
177  if (ssttmp[j][i] > -999)
178  sstref[j + 1][ii + 1] = ssttmp[j][i] * slope + offset;
179  else
180  sstref[j + 1][ii + 1] = sstbad;
181  }
182  sstref[j + 1][0] = sstref[j + 1][nx];
183  sstref[j + 1][nx + 1] = sstref[j + 1][1];
184  }
185  for (i = 0; i < nx + 2; i++) {
186  sstref[0] [i] = sstref[1][i];
187  sstref[ny + 1][i] = sstref[ny][i];
188  }
189 
190  free(ssttmp);
191  }
192 
193 
194  /* locate LL position within reference grid */
195  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), OI4NX + 1), 0);
196  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), OI4NY + 1), 0);
197 
198  /* compute longitude and latitude of that grid element */
199  xx = i * dx - 180.0 - dx / 2;
200  yy = j * dy - 90.0 - dy / 2;
201 
202  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
203  t = (lon - xx) / dx;
204  u = (lat - yy) / dy;
205 
206  ftmp = 0.0;
207  ntmp = 0;
208  if (sstref[j ][i ] > sstbad + 1) {
209  ftmp += sstref[j ][i ];
210  ntmp++;
211  }
212  if (sstref[j ][i + 1] > sstbad + 1) {
213  ftmp += sstref[j ][i + 1];
214  ntmp++;
215  }
216  if (sstref[j + 1][i + 1] > sstbad + 1) {
217  ftmp += sstref[j + 1][i + 1];
218  ntmp++;
219  }
220  if (sstref[j + 1][i ] > sstbad + 1) {
221  ftmp += sstref[j + 1][i ];
222  ntmp++;
223  }
224  if (ntmp > 0) {
225  ftmp /= ntmp;
226  reftmp[0][0] = (sstref[j ][i ] > sstbad + 1 ? sstref[j ][i ] : ftmp);
227  reftmp[0][1] = (sstref[j ][i + 1] > sstbad + 1 ? sstref[j ][i + 1] : ftmp);
228  reftmp[1][1] = (sstref[j + 1][i + 1] > sstbad + 1 ? sstref[j + 1][i + 1] : ftmp);
229  reftmp[1][0] = (sstref[j + 1][i ] > sstbad + 1 ? sstref[j + 1][i ] : ftmp);
230 
231  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
232 
233  } else
234  sst = sstbad;
235 
236  return (sst);
237 }
238 
239 /* ----------------------------------------------------------------------------------- */
240 /* get_cmcsst() - read and interpolate CMC SST netcdf files v2 and v3 */
241 /* ----------------------------------------------------------------------------------- */
242 float get_cmcsst(char *sstfile, float lon, float lat, int isKelvin) {
243  static int firstCall = 1;
244  static size_t nx = 0;
245  static size_t ny = 0;
246  static float dx = 0.1;
247  static float dy = 0.1;
248 
249  static float **sstref;
250  static float *latitudes;
251  static float *longitudes;
252 
253  float sst = sstbad;
254  int i, j;
255  int ntmp;
256  float xx, yy;
257  float t, u;
258  float reftmp[2][2];
259  float ftmp;
260 
261  if (firstCall) {
262  short **ssttmp;
263 
264  char sdsname[H4_MAX_NC_NAME] = "";
265  int32_t ncid, ndims, nvars, ngatts, unlimdimid;
266  int32_t latDimID, lonDimID;
267  int32 sds_id;
268  float slope;
269  float offset;
270 
271  firstCall = 0;
272 
273  printf("Loading SST reference from %s\n", sstfile);
274  printf("\n");
275 
276  /* Open the file */
277  if (nc_open(sstfile, NC_NOWRITE, &ncid) == NC_NOERR) {
278  nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
279 
280  strcpy(sdsname, "analysed_sst");
281  if (nc_inq_varid(ncid, sdsname, &sds_id) != NC_NOERR){
282  strcpy(sdsname, "sst");
283  if (nc_inq_varid(ncid, sdsname, &sds_id) != NC_NOERR){
284  printf("Whoops! something is wrong reading the SST reference file: %s\n", sstfile);
285  exit(EXIT_FAILURE);
286  }
287  }
288 
289  if (nc_inq_dimid(ncid, "lat", &latDimID) != NC_NOERR) {
290  printf("Whoops! something is wrong reading the SST reference file: %s\n", sstfile);
291  exit(EXIT_FAILURE);
292  }
293  nc_inq_dimlen(ncid, latDimID, &ny);
294  dy = 180.0 / ny;
295  if (nc_inq_dimid(ncid, "lon", &lonDimID) != NC_NOERR) {
296  printf("Whoops! something is wrong reading the SST reference file: %s\n", sstfile);
297  exit(EXIT_FAILURE);
298  }
299  nc_inq_dimlen(ncid, lonDimID, &nx);
300  dx = 360.0 / nx;
301 
302  /* size the arrays */
303  ssttmp = allocate2d_short(ny, nx);
304  sstref = allocate2d_float(ny+2, nx+2);
305 
306  latitudes = (float *) calloc(ny,sizeof(float));
307  longitudes = (float *) calloc(nx,sizeof(float));
308 
309  /* Read the data. */
310  if (nc_get_var_short(ncid, sds_id, &ssttmp[0][0]) != NC_NOERR) {
311  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
312  __FILE__, __LINE__, sdsname, sstfile);
313  exit(1);
314  }
315 
316  if (nc_get_att_float(ncid, sds_id, "scale_factor", &slope) != NC_NOERR) {
317  fprintf(stderr, "-E- %s line %d: error reading scale factor.\n",
318  __FILE__, __LINE__);
319  exit(1);
320  }
321 
322  if (nc_get_att_float(ncid, sds_id, "add_offset", &offset) != NC_NOERR) {
323  fprintf(stderr, "-E- %s line %d: error reading scale offset.\n",
324  __FILE__, __LINE__);
325  exit(1);
326  }
327  strcpy(sdsname, "lat");
328  nc_inq_varid(ncid, sdsname, &sds_id);
329  if (nc_get_var_float(ncid, sds_id, latitudes) != NC_NOERR) {
330  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
331  __FILE__, __LINE__, sdsname, sstfile);
332  exit(1);
333  }
334  strcpy(sdsname, "lon");
335  nc_inq_varid(ncid, sdsname, &sds_id);
336  if (nc_get_var_float(ncid, sds_id, longitudes) != NC_NOERR) {
337  fprintf(stderr, "-E- %s line %d: Error reading %s from %s.\n",
338  __FILE__, __LINE__, sdsname, sstfile);
339  exit(1);
340  }
341  /* Close the file */
342  if (nc_close(ncid) != NC_NOERR) {
343  fprintf(stderr, "-E- %s line %d: error closing %s.\n",
344  __FILE__, __LINE__, sstfile);
345  exit(1);
346  }
347  } else {
348  fprintf(stderr, "-E- %s line %d: error opening %s.\n",
349  __FILE__, __LINE__, sstfile);
350  exit(1);
351  }
352 
353  for (j = ny-1; j >= 0; j--) {
354  for (i = 0; i < nx; i++) {
355  if (ssttmp[j][i] > -999)
356  sstref[j + 1][i + 1] = ssttmp[j][i] * slope + offset;
357  else
358  sstref[j + 1][i + 1] = sstbad;
359  }
360  sstref[j + 1][0] = sstref[j + 1][nx];
361  sstref[j + 1][nx + 1] = sstref[j + 1][1];
362  }
363  for (i = 0; i < nx + 2; i++) {
364  sstref[0] [i] = sstref[1][i];
365  sstref[ny + 1][i] = sstref[ny][i];
366  }
367 
368  free2d_short(ssttmp);
369  }
370 
371  /* locate position within reference grid */
372  float sdiff = 999.0;
373  int lonidx = 0;
374  for (i=0; i < nx; i++){
375  float cdiff = fabs(longitudes[i] - lon);
376  if (cdiff > sdiff)
377  break;
378  if (cdiff < sdiff){
379  sdiff = cdiff;
380  lonidx = i;
381  }
382  }
383  sdiff = 999.0;
384  int latidx = 0;
385  for (i=0; i < ny; i++){
386  float cdiff = fabs(latitudes[i] - lat);
387  if (cdiff > sdiff)
388  break;
389  if (cdiff < sdiff){
390  sdiff = cdiff;
391  latidx = i;
392  }
393  }
394 
395  /* compute longitude and latitude of that grid element */
396  xx = longitudes[lonidx];
397  yy = latitudes[latidx];
398 
399  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
400  t = (lon - xx) / dx;
401  u = (lat - yy) / dy;
402 
403  ftmp = 0.0;
404  ntmp = 0;
405  if (sstref[latidx ][lonidx ] > sstbad + 1) {
406  ftmp += sstref[latidx ][lonidx ];
407  ntmp++;
408  }
409  if (sstref[latidx][lonidx + 1] > sstbad + 1) {
410  ftmp += sstref[latidx ][lonidx + 1];
411  ntmp++;
412  }
413  if (sstref[latidx + 1][lonidx + 1] > sstbad + 1) {
414  ftmp += sstref[latidx + 1][lonidx + 1];
415  ntmp++;
416  }
417  if (sstref[latidx + 1][lonidx ] > sstbad + 1) {
418  ftmp += sstref[latidx + 1][lonidx ];
419  ntmp++;
420  }
421  if (ntmp > 0) {
422  ftmp /= ntmp;
423  reftmp[0][0] = (sstref[latidx ][lonidx ] > sstbad + 1 ? sstref[latidx ][lonidx ] : ftmp);
424  reftmp[0][1] = (sstref[latidx ][lonidx + 1] > sstbad + 1 ? sstref[latidx ][lonidx + 1] : ftmp);
425  reftmp[1][1] = (sstref[latidx + 1][lonidx + 1] > sstbad + 1 ? sstref[latidx + 1][lonidx + 1] : ftmp);
426  reftmp[1][0] = (sstref[latidx + 1][lonidx ] > sstbad + 1 ? sstref[latidx + 1][lonidx ] : ftmp);
427 
428  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
429  if (isKelvin)
430  sst -= CtoK;
431 
432  } else
433  sst = sstbad;
434 
435  return (sst);
436 }
437 
438 
439 /* ----------------------------------------------------------------------------------- */
440 /* get_oisst() - read and interpolate flat binary Reynolds OI SST */
441 /* */
442 /* The files were written in IEEE binary (big-endian). Each file contains four FORTRAN */
443 /* records described as follows: */
444 /* */
445 /* rec 1: date and version number (8 4-byte integer words) */
446 /* rec 2: gridded sst values in degC (360*180 4-byte real words) */
447 /* rec 3: normalized error variance (360*180 4-byte real words) */
448 /* rec 4: gridded ice concentration (360*180 1-byte integer words) */
449 /* */
450 /* B. Franz, SAIC, May 2004. */
451 /* ----------------------------------------------------------------------------------- */
452 
453 #define OINX 360
454 #define OINY 180
455 
456 float get_oisst(char *sstfile, float lon, float lat) {
457  static int firstCall = 1;
458  static int nx = OINX;
459  static int ny = OINY;
460  static float dx = 360.0 / OINX;
461  static float dy = 180.0 / OINY;
462 
463  typedef float sstref_t[OINX + 2];
464  static sstref_t* sstref;
465  //static float sstref[OINY+2][OINX+2];
466 
467  float sst = BAD_FLT;
468  int i, j, ii;
469  float xx, yy;
470  float t, u;
471 
472  if (firstCall) {
473  sstref = (sstref_t*) malloc((OINY + 2) * sizeof (sstref_t));
474 
475  typedef float ssttmp_t[OINX];
476  ssttmp_t *ssttmp = (ssttmp_t*) malloc(OINY * sizeof (ssttmp_t));
477  //float ssttmp[OINY][OINX];
478 
479  FILE *fp = NULL;
480 
481  int32_t syear, smon, sday;
482  int32_t eyear, emon, eday;
483  int32_t ndays, version;
484 
485  firstCall = 0;
486 
487  if ((fp = fopen(sstfile, "r")) == NULL) {
488  printf("Error opening SST reference file %s for reading.\n", sstfile);
489  exit(1);
490  }
491 
492  if (fseek(fp, 4, SEEK_SET) < 0) {
493  printf("Error reading SST reference file %s.\n", sstfile);
494  exit(1);
495  }
496  if (fread(&syear, sizeof (int32_t), 1, fp) != 1) {
497  printf("Error reading SST reference file %s.\n", sstfile);
498  exit(1);
499  }
500  fread(&smon, sizeof (int32_t), 1, fp);
501  fread(&sday, sizeof (int32_t), 1, fp);
502  fread(&eyear, sizeof (int32_t), 1, fp);
503  fread(&emon, sizeof (int32_t), 1, fp);
504  fread(&eday, sizeof (int32_t), 1, fp);
505  fread(&ndays, sizeof (int32_t), 1, fp);
506  fread(&version, sizeof (int32_t), 1, fp);
507  fseek(fp, 4, SEEK_CUR);
508 
509  if (endianess() == 1) {
510  swapc_bytes((char *) &syear, sizeof (int32_t), 1);
511  swapc_bytes((char *) &smon, sizeof (int32_t), 1);
512  swapc_bytes((char *) &sday, sizeof (int32_t), 1);
513  swapc_bytes((char *) &eyear, sizeof (int32_t), 1);
514  swapc_bytes((char *) &emon, sizeof (int32_t), 1);
515  swapc_bytes((char *) &eday, sizeof (int32_t), 1);
516  swapc_bytes((char *) &ndays, sizeof (int32_t), 1);
517  swapc_bytes((char *) &version, sizeof (int32_t), 1);
518  }
519 
520  printf("Loading Weekly 1-degree OI Reynolds SST reference from %s\n", sstfile);
521  printf(" file start date (y/m/d): %d / %d / %d\n", syear, smon, sday);
522  printf(" file end date (y/m/d): %d / %d / %d\n", eyear, emon, eday);
523  printf(" days composited: %d\n", ndays);
524  printf(" file version: %d\n", version);
525  printf("\n");
526 
527  if (fseek(fp, 4, SEEK_CUR) < 0) {
528  printf("Error reading SST reference file %s.\n", sstfile);
529  exit(1);
530  }
531  if (fread(&ssttmp[0][0], sizeof (float), nx * ny, fp) != nx * ny) {
532  printf("Error reading SST reference file %s.\n", sstfile);
533  exit(1);
534  }
535  fclose(fp);
536 
537  if (endianess() == 1)
538  swapc_bytes((char *) &ssttmp[0][0], 4, nx * ny);
539 
540  /* rotate 180-deg and add wrapping border to simplify interpolation */
541  /* new grid is -180.5,180.5 in i=0,361 and -90.5,90.5 in j=0,181 */
542 
543  for (j = 0; j < ny; j++) {
544  for (i = 0; i < nx; i++) {
545  ii = (i < nx / 2) ? i + nx / 2 : i - nx / 2;
546  sstref[j + 1][ii + 1] = ssttmp[j][i];
547  }
548  sstref[j + 1][0] = sstref[j + 1][nx];
549  sstref[j + 1][nx + 1] = sstref[j + 1][1];
550  }
551  for (i = 0; i < nx + 2; i++) {
552  sstref[0] [i] = sstref[1][i];
553  sstref[ny + 1][i] = sstref[ny][i];
554  }
555  free(ssttmp);
556  }
557 
558 
559  /* locate LL position within reference grid */
560  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), OINX + 1), 0);
561  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), OINY + 1), 0);
562 
563  /* compute longitude and latitude of that grid element */
564  xx = i * dx - 180.0 - dx / 2;
565  yy = j * dy - 90.0 - dy / 2;
566 
567  /* bilinearly interpolate */
568  t = (lon - xx) / dx;
569  u = (lat - yy) / dy;
570 
571  sst = (1 - t)*(1 - u) * sstref[j ][i ]
572  + t * (1 - u) * sstref[j ][i + 1]
573  + t * u * sstref[j + 1][i + 1]
574  + (1 - t) * u * sstref[j + 1][i ];
575 
576  /*
577  sst = sstref[j][i];
578  */
579 
580  return (sst);
581 }
582 
583 /* ----------------------------------------------------------------------------------- */
584 /* get_ntev2() - read and interpolate Reynolds 0.25-deg daily NTEV2 file */
585 /* */
586 /* B. Franz, SAIC, July 2008. */
587 /* ----------------------------------------------------------------------------------- */
588 
589 #define OI4NX 1440
590 #define OI4NY 720
591 
592 float get_ntev2(char *sstfile, float lon, float lat, int32_t year, int32_t day) {
593  static int firstCall = 1;
594  static int nx = OI4NX;
595  static int ny = OI4NY;
596  static float dx = 360.0 / OI4NX;
597  static float dy = 180.0 / OI4NY;
598 
599  typedef float sstref_t[OI4NX + 2];
600  static sstref_t* sstref;
601  //static float sstref[OI4NY+2][OI4NX+2];
602  //static float ssttmp[OI4NY][OI4NX];
603 
604 
605  float sst = sstbad;
606  int i, j, ii;
607  int ntmp;
608  int32_t syear, smon, sday;
609  float xx, yy;
610  float t, u;
611  float reftmp[2][2];
612  float ftmp;
613 
614  if (firstCall) {
615  sstref = (sstref_t*) malloc((OI4NY + 2) * sizeof (sstref_t));
616 
617  typedef float ssttmp_t[OI4NX];
618  ssttmp_t *ssttmp;
619  ssttmp = (ssttmp_t*) malloc(OI4NY * sizeof (ssttmp_t));
620 
621  /* daily file for years 2006..2009, 1440x720 */
622  FILE *fp = NULL;
623  firstCall = 0;
624 
625  /* just need first 1440x720 real*4 array of data from the daily file */
626  /* iyr,imo,ida,refval */
627 
628  if (year < 2006 || year > 2009) {
629  printf("ntev2 data only for years 2006 to 2009, not %d\n", year);
630  exit(1);
631  }
632  printf("Loading ntev2 daily field from %s\n", sstfile);
633  printf("\n");
634 
635  /* Open the file */
636  if ((fp = fopen(sstfile, "r")) == NULL) {
637  printf("-E- %s line %d: error opening SST reference file %s for reading.\n",
638  __FILE__, __LINE__, sstfile);
639  exit(1);
640  }
641 
642  if (fseek(fp, 4, SEEK_SET) < 0) {
643  printf("Error seeking SST NTEV2 reference file %s.\n", sstfile);
644  exit(1);
645  }
646  if (fread(&syear, sizeof (int32_t), 1, fp) != 1) {
647  printf("Error reading SST NTEV2 reference file %s.\n", sstfile);
648  exit(1);
649  }
650  fread(&smon, sizeof (int32_t), 1, fp);
651  fread(&sday, sizeof (int32_t), 1, fp);
652 
653  if (fread(&ssttmp[0][0], sizeof (float), nx * ny, fp) != nx * ny) {
654  printf("Error reading SST NTEV2 reference file %s.\n", sstfile);
655  exit(1);
656  }
657  fclose(fp);
658 
659  if (endianess() == 1) {
660  swapc_bytes((char *) &syear, sizeof (int32_t), 1);
661  swapc_bytes((char *) &smon, sizeof (int32_t), 1);
662  swapc_bytes((char *) &sday, sizeof (int32_t), 1);
663  swapc_bytes((char *) &ssttmp[0][0], sizeof (float), nx * ny);
664  }
665 
666  printf("Loading Daily 0.25 Deg NTEV2 Reynolds SST reference from %s\n", sstfile);
667  printf(" file date (y/m/d): %d / %d / %d\n", syear, smon, sday);
668  printf("\n");
669 
670 
671  /* rotate 180-deg and add wrapping border to simplify later interpolation between points */
672  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
673 
674  for (j = 0; j < ny; j++) {
675  for (i = 0; i < nx; i++) {
676  ii = (i < nx / 2) ? i + nx / 2 : i - nx / 2;
677  sstref[j + 1][ii + 1] = ssttmp[j][i];
678  }
679  sstref[j + 1][0] = sstref[j + 1][nx];
680  sstref[j + 1][nx + 1] = sstref[j + 1][1];
681  }
682  for (i = 0; i < nx + 2; i++) {
683  sstref[0] [i] = sstref[1][i];
684  sstref[ny + 1][i] = sstref[ny][i];
685  }
686  free(ssttmp);
687  }
688 
689 
690  /* locate LL position within reference grid */
691  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), OI4NX + 1), 0);
692  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), OI4NY + 1), 0);
693 
694  /* compute longitude and latitude of that grid element */
695  xx = i * dx - 180.0 - dx / 2;
696  yy = j * dy - 90.0 - dy / 2;
697 
698  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
699  t = (lon - xx) / dx;
700  u = (lat - yy) / dy;
701 
702  ftmp = 0.0;
703  ntmp = 0;
704  if (sstref[j ][i ] > sstbad + 1) {
705  ftmp += sstref[j ][i ];
706  ntmp++;
707  }
708  if (sstref[j ][i + 1] > sstbad + 1) {
709  ftmp += sstref[j ][i + 1];
710  ntmp++;
711  }
712  if (sstref[j + 1][i + 1] > sstbad + 1) {
713  ftmp += sstref[j + 1][i + 1];
714  ntmp++;
715  }
716  if (sstref[j + 1][i ] > sstbad + 1) {
717  ftmp += sstref[j + 1][i ];
718  ntmp++;
719  }
720  if (ntmp > 0) {
721  ftmp /= ntmp;
722  reftmp[0][0] = (sstref[j ][i ] > sstbad + 1 ? sstref[j ][i ] : ftmp);
723  reftmp[0][1] = (sstref[j ][i + 1] > sstbad + 1 ? sstref[j ][i + 1] : ftmp);
724  reftmp[1][1] = (sstref[j + 1][i + 1] > sstbad + 1 ? sstref[j + 1][i + 1] : ftmp);
725  reftmp[1][0] = (sstref[j + 1][i ] > sstbad + 1 ? sstref[j + 1][i ] : ftmp);
726 
727  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
728 
729  } else
730  sst = sstbad;
731 
732  return (sst);
733 }
734 
735 /* ----------------------------------------------------------------------------------- */
736 /* get_astr() - read and interpolate Reynolds 0.25-deg monthly ATSR file */
737 /* */
738 /* B. Franz, SAIC, July 2008. */
739 /* ----------------------------------------------------------------------------------- */
740 
741 #define OI4NX 1440
742 #define OI4NY 720
743 
744 float get_atsr(char *sstfile, float lon, float lat, int32_t year, int32_t day, int32_t minatsrcnt) {
745  static int firstCall = 1;
746  static int nx = OI4NX;
747  static int ny = OI4NY;
748  static int month;
749  static float dx = 360.0 / OI4NX;
750  static float dy = 180.0 / OI4NY;
751 
752  typedef float sstref_t[OI4NX + 2];
753  static sstref_t* sstref;
754  //static float sstref[OI4NY+2][OI4NX+2];
755  //static float ssttmpb[2][OI4NY][OI4NX];
756  //static float sstcnt[2][OI4NY][OI4NX];
757 
758  float sst = sstbad;
759  int i, j, ii;
760  int ntmp;
761  float xx, yy;
762  float t, u;
763  float reftmp[2][2];
764  float ftmp;
765 
766  if (firstCall) {
767  sstref = (sstref_t*) malloc((OI4NY + 2) * sizeof (sstref_t));
768 
769  typedef float ssttmp_t[OI4NY][OI4NX];
770  ssttmp_t* ssttmpb = (ssttmp_t*) malloc(2 * sizeof (ssttmp_t));
771  ssttmp_t* sstcnt = (ssttmp_t*) malloc(2 * sizeof (ssttmp_t));
772 
773  /* years 2006..2009, months 1..12, 1440x720 */
774  FILE *fp = NULL;
775  int otherleap;
776  int endofyear;
777  int numdays;
778  int32_t foff[2];
779  int32 status;
780  int leap, day1, day2;
781  int32_t reclen;
782  float w1, w2;
783  firstCall = 0;
784 
785  /* 4*12*3*(4+4+4+4+(4*1440*720)+4) */
786  /* years 2006..2009; months 1..12; mean,std,count */
787  /* 4 bytes of fortran record length; year,month,count (3 integers) */
788  /* floats 1440x720 data; 4 byte fortran record end */
789  /* use seek_set to find first month wanted? */
790  /* use seek_cur to jump to next months mean? */
791 
792  if (year < 2006 || year > 2009) {
793  printf("ATSR data only for years 2006 to 2009, not %d\n", year);
794  exit(1);
795  }
796  printf("Loading ATSR monthly fields from %s\n", sstfile);
797  printf("\n");
798  leap = (isleap(year) == TRUE ? 1 : 0);
799  for (month = 11; month >= 0; month--) {
800  if (day > MiddleOfMonth[leap][month]) {
801  break;
802  }
803  }
804  /* each "record" is recbeg,iyr,imo,ndct,data,recend */
805  /* which is 4+4+4+4+4*1440*720+4 = 4147220 bytes */
806  /* years go from 2006 to 2009; months from 1 to 12 */
807  /* mean is first of 3 data records so read one,skip two,read one */
808  /* Jan 2006 mean is record 1, Feb 2006 mean is record 4, ... */
809 
810  /* read in two months and interpolate */
811  reclen = 4147220; /* number of bytes in one record */
812 
813  /* linearly interpolate the reference values from the two months around day */
814  /* w1 is the weight to use with the ref value from the earlier month */
815  /* w2 is the weight to use with the ref value from the later month */
816  /* The next months mean is in the next record after the previous months counts so set the 2nd offset */
817  /* to skip the fortran end of record byte count, beg of rec byte count, year, month, count before the means */
818  /* foff[1] will be set to -1 when day is in the first half of the first month or the */
819  /* last half of the last month and therefore there is no previous or next month to read */
820  foff[1] = 0; /* last read is now counts, so next mean is next */
821  if (month == -1) {
822  /* day is in first half of Jan, interpolate with previous Dec */
823  if (year == 2006) {
824  /* There is no "previous Dec." so use first set as is */
825  w1 = 1.0;
826  w2 = 0.0;
827  /* foff[0] is offset to means for first month whose center is before day */
828  /* in this case it's the first data record so no offset */
829  /* except the beg,year,mon,count that's before the data */
830  foff[0] = 0;
831  /* there is no previous month so set 2nd offset to be zero so no read is done */
832  foff[1] = -1;
833  } else {
834  /* day is in first half of Jan, interpolate with previous Dec */
835  otherleap = (isleap(year - 1) == TRUE ? 1 : 0);
836  endofyear = (otherleap == 1 ? 366 : 365);
837  numdays = ((endofyear - MiddleOfMonth[otherleap][11]) +
838  MiddleOfMonth[leap][0]);
839  w1 = (float) (MiddleOfMonth[leap][0] - day) / (float) numdays;
840  w2 = 1.0 - w1;
841  /* foff[0] is offset to means for first month whose center is before day */
842  /* in this case it's Dec of the previous year */
843  /* data is after 16 bytes of beg+year+month+count */
844  foff[0] = ((year - 1 - 2006)*12 + 11) * reclen * 3;
845  }
846  } else if (month == 11) {
847  /* day is in last half of Dec, interpolate with next Jan */
848  if (year == 2009) {
849  /* there is no "next Jan." so just use Dec. coeffs as is */
850  w1 = 1.0;
851  w2 = 0.0;
852  /* foff[0] is offset to means for first month whose center is before day */
853  /* in this case it's Dec of the last year, so there is no next month */
854  /* plus the stuff at the beginning of each record: beg,year,month,count */
855  foff[0] = ((year - 2006)*12 + 11) * reclen * 3;
856  foff[1] = -1;
857  } else {
858  /* day is in last half of Dec, interpolate with next Jan */
859  otherleap = (isleap(year + 1) == TRUE ? 1 : 0);
860  endofyear = (leap == 1 ? 366 : 365);
861  numdays = ((endofyear - MiddleOfMonth[leap][11]) +
862  MiddleOfMonth[otherleap][0]);
863  w2 = (float) (day - MiddleOfMonth[leap][11]) / (float) numdays;
864  w1 = 1.0 - w2;
865  /* foff[0] is offset to means for first month whose center is before day */
866  /* in this case the first month to read is this Dec */
867  /* plus fortran record count,year,month,count=16 bytes before data starts */
868  foff[0] = ((year - 2006)*12 + 11) * reclen * 3;
869  }
870  } else {
871  /* interpolate with next month */
872  day1 = MiddleOfMonth[leap][month];
873  day2 = MiddleOfMonth[leap][month + 1];
874  w2 = (float) (day - day1) / (float) (day2 - day1);
875  w1 = 1.0 - w2;
876  /* foff[0] is offset to means for first month whose center is before day */
877  /* in this case it's a "normal" month in the "middle" of a year */
878  /* plus fortran record count,year,month,count=16 bytes before data starts */
879  foff[0] = ((year - 2006)*12 + month) * reclen * 3;
880  }
881  /* Open the file */
882  if ((fp = fopen(sstfile, "r")) == NULL) {
883  printf("-E- %s line %d: error opening SST reference file %s for reading.\n",
884  __FILE__, __LINE__, sstfile);
885  exit(1);
886  }
887  /* skip to the first month (and skip fortran begin record count,year,month,and day count) */
888  if (fseek(fp, foff[0] + 16, SEEK_SET) < 0) {
889  printf("Error seeking SST reference file %s.\n", sstfile);
890  exit(1);
891  }
892 
893  if ((status = fread(&ssttmpb[0][0][0], sizeof (float), nx * ny, fp)) != nx * ny) {
894  printf("Wrong atsr data read size: want %d got %d\n", nx*ny, status);
895  exit(1);
896  }
897  /* skip fortran end of record, std. dev., fortran beg of record to read counts,year,month,day count */
898  if (fseek(fp, 4 + reclen + 16, SEEK_CUR) < 0) {
899  printf("Error seeking first counts from SST reference file %s.\n", sstfile);
900  exit(1);
901  }
902  if ((status = fread(&sstcnt[0][0][0], sizeof (float), nx * ny, fp)) != nx * ny) {
903  printf("Wrong atsr counts data read size: want %d got %d\n", nx*ny, status);
904  exit(1);
905  }
906  if (foff[1] == 0) {
907  /* read next months mean */
908  /* skip fortran end of record count, begin of record count,year,month, and day count */
909  fseek(fp, 20, SEEK_CUR);
910  if ((status = fread(&ssttmpb[1][0][0], sizeof (float), nx * ny, fp)) != nx * ny) {
911  printf("Wrong 2nd atsr data read size: want %d got %d\n", nx*ny, status);
912  exit(1);
913  }
914  /* skip fortran end of record count for means, std. dev., begin of record count,year,month,day count to read counts */
915  if (fseek(fp, 4 + reclen + 16, SEEK_CUR) < 0) {
916  printf("Error seek to second counts from SST reference file %s.\n", sstfile);
917  exit(1);
918  }
919  if ((status = fread(&sstcnt[1][0][0], sizeof (float), nx * ny, fp)) != nx * ny) {
920  printf("Wrong 2nd atsr counts data read size: want %d got %d\n", nx*ny, status);
921  exit(1);
922  }
923  }
924 
925  if (endianess() == 1) {
926  swapc_bytes((char *) &ssttmpb[0][0][0], sizeof (float), 2 * nx * ny);
927  swapc_bytes((char *) &sstcnt[0][0][0], sizeof (float), 2 * nx * ny);
928  }
929 
930  fclose(fp);
931  /* interpolate between the two months in ssttmpb */
932  /* rotate 180-deg and add wrapping border to simplify later interpolation between points */
933  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
934 
935  for (j = 0; j < ny; j++) {
936  for (i = 0; i < nx; i++) {
937  ii = (i < nx / 2) ? i + nx / 2 : i - nx / 2;
938  if (ssttmpb[0][j][i] > -999.0 && ssttmpb[1][j][i] > -999.0 &&
939  sstcnt[0][j][i] >= minatsrcnt && sstcnt[1][j][i] >= minatsrcnt)
940  sstref[j + 1][ii + 1] = w1 * ssttmpb[0][j][i] + w2 * ssttmpb[1][j][i];
941  else {
942  if (ssttmpb[0][j][i] > -999.0 && sstcnt[0][j][i] >= minatsrcnt)
943  sstref[j + 1][ii + 1] = ssttmpb[0][j][i];
944  else if (ssttmpb[1][j][i] > -999.0 && sstcnt[1][j][i] >= minatsrcnt)
945  sstref[j + 1][ii + 1] = ssttmpb[1][j][i];
946  else
947  sstref[j + 1][ii + 1] = sstbad;
948  }
949  }
950  sstref[j + 1][0] = sstref[j + 1][nx];
951  sstref[j + 1][nx + 1] = sstref[j + 1][1];
952  }
953  for (i = 0; i < nx + 2; i++) {
954  sstref[0] [i] = sstref[1][i];
955  sstref[ny + 1][i] = sstref[ny][i];
956  }
957  free(ssttmpb);
958  free(sstcnt);
959  }
960 
961 
962  /* locate LL position within reference grid */
963  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), OI4NX + 1), 0);
964  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), OI4NY + 1), 0);
965 
966  /* compute longitude and latitude of that grid element */
967  xx = i * dx - 180.0 - dx / 2;
968  yy = j * dy - 90.0 - dy / 2;
969 
970  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
971  t = (lon - xx) / dx;
972  u = (lat - yy) / dy;
973 
974  ftmp = 0.0;
975  ntmp = 0;
976  if (sstref[j ][i ] > sstbad + 1) {
977  ftmp += sstref[j ][i ];
978  ntmp++;
979  }
980  if (sstref[j ][i + 1] > sstbad + 1) {
981  ftmp += sstref[j ][i + 1];
982  ntmp++;
983  }
984  if (sstref[j + 1][i + 1] > sstbad + 1) {
985  ftmp += sstref[j + 1][i + 1];
986  ntmp++;
987  }
988  if (sstref[j + 1][i ] > sstbad + 1) {
989  ftmp += sstref[j + 1][i ];
990  ntmp++;
991  }
992  if (ntmp > 0) {
993  ftmp /= ntmp;
994  reftmp[0][0] = (sstref[j ][i ] > sstbad + 1 ? sstref[j ][i ] : ftmp);
995  reftmp[0][1] = (sstref[j ][i + 1] > sstbad + 1 ? sstref[j ][i + 1] : ftmp);
996  reftmp[1][1] = (sstref[j + 1][i + 1] > sstbad + 1 ? sstref[j + 1][i + 1] : ftmp);
997  reftmp[1][0] = (sstref[j + 1][i ] > sstbad + 1 ? sstref[j + 1][i ] : ftmp);
998 
999  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
1000 
1001  } else
1002  sst = sstbad;
1003 
1004  return (sst);
1005 }
1006 
1007 /* ----------------------------------------------------------------------------------- */
1008 /* get_amsre() - read and interpolate AMSR-E data from SSMI */
1009 /* */
1010 /* S. Walsh, RSMAS, Aug 2010. */
1011 /* ----------------------------------------------------------------------------------- */
1012 
1013 #define OI4NX 1440
1014 #define OI4NY 720
1015 
1016 float get_amsre(char *sstfile, float lon, float lat, float solz, int32_t xsatid, int32_t sensorID) {
1017  static int firstCall = 1;
1018  static int nx = OI4NX;
1019  static int ny = OI4NY;
1020  static float dx = 360.0 / OI4NX;
1021  static float dy = 180.0 / OI4NY;
1022 
1023  typedef float sstref_t[OI4NX + 2][2];
1024  static sstref_t* sstref;
1025  //static float sstref[OI4NY+2][OI4NX+2][2];
1026 
1027  // char dbgfile[FILENAME_MAX] = "./sstref.dbg.file";
1028 
1029  float sst = sstbad;
1030  int jj, ii, kk, ij;
1031  int ad;
1032  int kkmax;
1033  int ntmp;
1034  int32_t foff[2];
1035  float xx, yy;
1036  float t, u;
1037  float reftmp[2][2];
1038  float ftmp;
1039 
1040  if (firstCall) {
1041  sstref = (sstref_t*) malloc((OI4NY + 2) * sizeof (sstref_t));
1042 
1043  FILE *fp = NULL;
1044  unsigned char ssttmp[OI4NX];
1045  int32 status;
1046 
1047  firstCall = 0;
1048 
1049  if (format == AMSRE3DAY)
1050  printf("Loading 3-Day 0.25-deg AMSR-E SST reference from %s\n", sstfile);
1051  else if (format == AMSRE3DN)
1052  printf("Loading 3-Day Day or Night 0.25-deg AMSR-E SST reference from %s\n", sstfile);
1053  else
1054  printf("Loading Daily 0.25-deg AMSR-E SST reference from %s\n", sstfile);
1055  printf("\n");
1056 
1057  /* Open the file */
1058  if ((fp = fopen(sstfile, "r")) == NULL) {
1059  printf("-E- %s line %d: error opening SST reference file %s for reading.\n",
1060  __FILE__, __LINE__, sstfile);
1061  exit(1);
1062  }
1063 
1064  if (format == AMSRE3DAY) {
1065  /* 3-Day - first field is 3 day avg of sst asc and dsc */
1066  foff[0] = 0;
1067  /* only one set for both ascending and descending */
1068  kkmax = 1;
1069  } else if (format == AMSRE3DN) {
1070  /* 3-Day - Day or Night - first field is 3 day avg of sst asc (day for aqua/amsr) and 2nd field is dsc night */
1071  foff[0] = 0;
1072  foff[1] = (int32_t) (nx * ny);
1073  /* one set of data for ascending, and one set for descending */
1074  kkmax = 2;
1075  } else {
1076  /* Daily - second field is 3 day avg of sst asc and dsc */
1077  /* v5 has 6 fields for day and 6 for night so fields 2 and 8 are sst */
1078  /* v7 has 7 fields for day and 7 for night so fields 2 and 9 are sst */
1079  foff[0] = (int32_t) (nx * ny);
1080  foff[1] = (int32_t) (nx * ny * 8); /* was 7 for v5 */
1081  /* one set of data for ascending, and one set for descending */
1082  kkmax = 2;
1083  }
1084 
1085  for (kk = 0; kk < kkmax; kk = kk + 1) {
1086  /* only read sst */
1087  /* ascending (day) is first, then descending (night) */
1088  if (fseek(fp, foff[kk], SEEK_SET) < 0) {
1089  fprintf(stderr, "-E- %s line %d: error reading SST reference file %s.\n",
1090  __FILE__, __LINE__, sstfile);
1091  exit(1);
1092  }
1093  for (jj = 0; jj < ny; jj = jj + 1) {
1094  if ((status = fread(&ssttmp[0], sizeof (char), nx, fp)) != nx) {
1095  printf("Wrong AMSR-E data read size: want %d got %d\n", nx, status);
1096  exit(1);
1097  }
1098  /* 0-250 are valid values */
1099  for (ii = 0; ii < nx; ii = ii + 1) {
1100  /* rotate 180-deg and add wrapping border to simplify interpolation */
1101  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
1102 
1103  ij = (ii < nx / 2) ? ii + nx / 2 : ii - nx / 2;
1104  if (ssttmp[ii] <= 250) {
1105  sstref[jj + 1][ij + 1][kk] = ssttmp[ii] * 0.15 - 3.0;
1106  } else {
1107  sstref[jj + 1][ij + 1][kk] = sstbad;
1108  }
1109  /* if d3d file put values in both asc and dsc slots */
1110  if (kkmax == 1) {
1111  sstref[jj + 1][ij + 1][1] = sstref[jj + 1][ij + 1][0];
1112  }
1113  }
1114 
1115  /* wrap edges */
1116  sstref[jj + 1][0] [kk] = sstref[jj + 1][nx][kk];
1117  sstref[jj + 1][nx + 1][kk] = sstref[jj + 1][1] [kk];
1118  }
1119  /* copy edges for interpolation */
1120  for (ii = 0; ii < nx + 2; ii = ii + 1) {
1121  sstref[0] [ii][kk] = sstref[1] [ii][kk];
1122  sstref[ny + 1][ii][kk] = sstref[ny][ii][kk];
1123  }
1124  }
1125  fclose(fp);
1126 
1127  }
1128 
1129  if (solz < 90.0) {
1130  ad = 0; /* day - ascending for amsr on aqua */
1131  } else {
1132  ad = 1; /* night - descending for amsr on aqua */
1133  }
1134 
1135  // morning satellites, daytime pixels use night amsr
1136  // afternoon satellites, daytime pixels use day amsr
1137 
1138  // /* some day may want to swap day/night for satellites that
1139  // * descend in the daytime so that the swath are the same direction
1140  // * amsr was on aqua so day is 1:30 pm local time, afternoon, lots of daytime warming
1141  // * so morning satellites, whose "day" pixels are before daytime warning,
1142  // * should use amsr night for their reference field instead of amsr (afternoon) day
1143  // */
1144 
1145  if ((format == AMSREDAY || format == AMSRE3DAY || format == AMSRE3DN) &&
1146  ((sensorID == MODIST ||
1147  (sensorID == AVHRR && (xsatid == NO06 || xsatid == NO08 ||
1148  xsatid == NO10 || xsatid == NO12 || xsatid == NO15 || xsatid == NO17))))) {
1149  /* These satellites descend in day, opposite of aqua */
1150  /* and their "day" pixels are 9:30 am, before daytime warming, so use amsr night field */
1151  ad = 1 - ad;
1152  }
1153 
1154  /* locate LL position within reference grid */
1155  ii = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
1156  jj = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
1157 
1158  /* compute longitude and latitude of that grid element */
1159  xx = ii * dx - 180.0 - dx / 2;
1160  yy = jj * dy - 90.0 - dy / 2;
1161 
1162  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
1163  t = (lon - xx) / dx;
1164  u = (lat - yy) / dy;
1165 
1166  ftmp = 0.0;
1167  ntmp = 0;
1168  if (sstref[jj ][ii ][ad] > sstbad + 1) {
1169  ftmp += sstref[jj ][ii ][ad];
1170  ntmp++;
1171  }
1172  if (sstref[jj ][ii + 1][ad] > sstbad + 1) {
1173  ftmp += sstref[jj ][ii + 1][ad];
1174  ntmp++;
1175  }
1176  if (sstref[jj + 1][ii + 1][ad] > sstbad + 1) {
1177  ftmp += sstref[jj + 1][ii + 1][ad];
1178  ntmp++;
1179  }
1180  if (sstref[jj + 1][ii ][ad] > sstbad + 1) {
1181  ftmp += sstref[jj + 1][ii ][ad];
1182  ntmp++;
1183  }
1184  if (ntmp > 0) {
1185  ftmp /= ntmp;
1186  reftmp[0][0] = (sstref[jj ][ii ][ad] > sstbad + 1 ? sstref[jj ][ii ][ad] : ftmp);
1187  reftmp[0][1] = (sstref[jj ][ii + 1][ad] > sstbad + 1 ? sstref[jj ][ii + 1][ad] : ftmp);
1188  reftmp[1][1] = (sstref[jj + 1][ii + 1][ad] > sstbad + 1 ? sstref[jj + 1][ii + 1][ad] : ftmp);
1189  reftmp[1][0] = (sstref[jj + 1][ii ][ad] > sstbad + 1 ? sstref[jj + 1][ii ][ad] : ftmp);
1190 
1191  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
1192 
1193  } else
1194  sst = sstbad;
1195 
1196  return (sst);
1197 }
1198 
1199 /* ----------------------------------------------------------------------------------- */
1200 /* get_windsat() - read and interpolate WindSat data from SSMI */
1201 /* */
1202 /* S. Walsh, RSMAS, Aug 2010. */
1203 /* ----------------------------------------------------------------------------------- */
1204 
1205 #define OI4NX 1440
1206 #define OI4NY 720
1207 
1208 float get_windsat(char *sstfile, float lon, float lat, float solz, int32_t xsatid, int32_t sensorID) {
1209  static int firstCall = 1;
1210  static int nx = OI4NX;
1211  static int ny = OI4NY;
1212  static float dx = 360.0 / OI4NX;
1213  static float dy = 180.0 / OI4NY;
1214 
1215  typedef float sstref_t[OI4NX + 2][2];
1216  static sstref_t* sstref;
1217  //static float sstref[OI4NY+2][OI4NX+2][2];
1218 
1219  // char dbgfile[FILENAME_MAX] = "./sstref.dbg.file";
1220 
1221  float sst = sstbad;
1222  int jj, ii, kk, ij;
1223  int ad;
1224  int kkmax;
1225  int ntmp;
1226  int32_t foff[2];
1227  float xx, yy;
1228  float t, u;
1229  float reftmp[2][2];
1230  float ftmp;
1231 
1232  if (firstCall) {
1233  sstref = (sstref_t*) malloc((OI4NY + 2) * sizeof (sstref_t));
1234 
1235  FILE *fp = NULL;
1236  unsigned char ssttmp[OI4NX];
1237  float fssttmp[OI4NX];
1238  int32 status;
1239 
1240  firstCall = 0;
1241 
1242  if (format == WINDSAT3DAY)
1243  printf("Loading 3-Day 0.25-deg WindSat SST reference from %s\n", sstfile);
1244  else if (format == WINDSAT3DN)
1245  printf("Loading 3-Day Day or Night 0.25-deg WindSat SST reference from %s\n", sstfile);
1246  else
1247  printf("Loading Daily 0.25-deg WindSat SST reference from %s\n", sstfile);
1248  printf("\n");
1249 
1250  /* Open the file */
1251  if ((fp = fopen(sstfile, "r")) == NULL) {
1252  printf("-E- %s line %d: error opening SST reference file %s for reading.\n",
1253  __FILE__, __LINE__, sstfile);
1254  exit(1);
1255  }
1256 
1257  if (format == WINDSAT3DAY) {
1258  /* 3-Day - first field is 3 day avg of sst dsc (day) and asc (night) */
1259  foff[0] = 0;
1260  /* only one set for both descending and ascending */
1261  kkmax = 1;
1262  } else if (format == WINDSAT3DN) {
1263  /* 3-Day - Day or Night - first field is 3 day avg of sst dsc and 2nd field is asc night */
1264  /* values are floats */
1265  foff[0] = 0;
1266  foff[1] = (int32_t) (nx * ny * sizeof (float));
1267  /* one set of data for descending, and one set for ascending */
1268  kkmax = 2;
1269  } else {
1270  /* Daily - second field is 3 day avg of sst desc and asc */
1271  /* 9 fields for desc, 9 for asc, sst's are bands 2 and 11 (one based) */
1272  /* values are bytes */
1273  foff[0] = (int32_t) (nx * ny);
1274  foff[1] = (int32_t) (nx * ny * 10);
1275  /* seek_set starts at 0 each time, so skip fields 1..10 */
1276  /* one set of data for descending, and one set for ascending */
1277  kkmax = 2;
1278  }
1279 
1280  for (kk = 0; kk < kkmax; kk = kk + 1) {
1281  /* only read sst */
1282  /* descending (morning) is first, then ascending (evening) */
1283  if (fseek(fp, foff[kk], SEEK_SET) < 0) {
1284  fprintf(stderr, "-E- %s line %d: error reading SST reference file %s.\n",
1285  __FILE__, __LINE__, sstfile);
1286  exit(1);
1287  }
1288  for (jj = 0; jj < ny; jj = jj + 1) {
1289  if (format == WINDSAT3DN) {
1290  if ((status = fread(&fssttmp[0], sizeof (float), nx, fp)) != nx) {
1291  printf("Wrong WINDSAT data read size: want %d got %d\n", nx, status);
1292  exit(1);
1293  }
1294  // if (jj == 0) {
1295  // printf(" first windsat values: %f %f %f %f %f %f %f\n", fssttmp[0],fssttmp[1],fssttmp[2],fssttmp[3],fssttmp[4],fssttmp[5],fssttmp[6]);
1296  // }
1297  if (endianess() == 1) {
1298  swapc_bytes((char *) &fssttmp[0], sizeof (float), nx);
1299  }
1300  // if (jj == 220) {
1301  // printf(" fssttmp[1164]=%f fssttmp[1165]=%f\n",fssttmp[1164],fssttmp[1165]);
1302  // printf(" fssttmp[444]=%f fssttmp[445]=%f\n",fssttmp[444],fssttmp[445]);
1303  // }
1304  // if (jj == 0) {
1305  // printf(" first windsat values: %f %f %f %f %f %f %f\n", fssttmp[0],fssttmp[1],fssttmp[2],fssttmp[3],fssttmp[4],fssttmp[5],fssttmp[6]);
1306  // }
1307  } else {
1308  if ((status = fread(&ssttmp[0], sizeof (char), nx, fp)) != nx) {
1309  printf("Wrong WINDSAT data read size: want %d got %d\n", nx, status);
1310  exit(1);
1311  }
1312  }
1313  /* 0-250 are valid values */
1314  for (ii = 0; ii < nx; ii = ii + 1) {
1315  /* rotate 180-deg and add wrapping border to simplify interpolation */
1316  /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */
1317 
1318  ij = (ii < nx / 2) ? ii + nx / 2 : ii - nx / 2;
1319  if (format == WINDSAT3DN) {
1320  /* for our 3 day average, bad values are -9999.0 */
1321  if (fssttmp[ii] > -9998.0) {
1322  sstref[jj + 1][ij + 1][kk] = fssttmp[ii];
1323  } else {
1324  sstref[jj + 1][ij + 1][kk] = sstbad;
1325  }
1326  } else {
1327  if (ssttmp[ii] <= 250) {
1328  sstref[jj + 1][ij + 1][kk] = ssttmp[ii] * 0.15 - 3.0;
1329  } else {
1330  sstref[jj + 1][ij + 1][kk] = sstbad;
1331  }
1332  }
1333  /* if d3d file put values in both asc and dsc slots */
1334  if (kkmax == 1) {
1335  sstref[jj + 1][ij + 1][1] = sstref[jj + 1][ij + 1][0];
1336  }
1337  }
1338 
1339  /* wrap edges */
1340  sstref[jj + 1][0] [kk] = sstref[jj + 1][nx][kk];
1341  sstref[jj + 1][nx + 1][kk] = sstref[jj + 1][1] [kk];
1342  }
1343  /* copy edges for interpolation */
1344  for (ii = 0; ii < nx + 2; ii = ii + 1) {
1345  sstref[0] [ii][kk] = sstref[1] [ii][kk];
1346  sstref[ny + 1][ii][kk] = sstref[ny][ii][kk];
1347  }
1348  }
1349  fclose(fp);
1350 
1351 
1352  }
1353 
1354  /* default is for afternoon satellites whose day pixels are around 1:30 pm local time */
1355  /* their day pixels should be compared to windsat evening pixels, and their night to windsat morning */
1356  if (solz < 90.0) {
1357  ad = 1; /* day pixel, use windsat evening - aescending */
1358  } else {
1359  ad = 0; /* night pixel, use windsat morning - descending */
1360  }
1361 
1362  // windsat fields are morning, descending (08:45) and evening, ascending (20:45)
1363  // viirs is descending at 13:45 and ascending at 01:45
1364  // so for viirs descending day, use evening ascending windsat
1365  // morning fields are best for night pixels from afternoon satellites because no daytime warming has occurred yet
1366  // evening fields are best for afternoon daytime pixels
1367 
1368  // morning satellites, daytime pixels use morning windsat
1369  // afternoon satellites, daytime pixels use evening windsat
1370 
1371  // /* some day may want to swap day/night for satellites that
1372  // * descend in the daytime so that the swath are the same direction
1373  // */
1374  if ((format == WINDSATDAY || format == WINDSAT3DAY || format == WINDSAT3DN) &&
1375  ((sensorID == MODIST ||
1376  (sensorID == AVHRR && (xsatid == NO06 || xsatid == NO08 ||
1377  xsatid == NO10 || xsatid == NO12 || xsatid == NO15 || xsatid == NO17))))) {
1378  /* These satellites descend in day, opposite of aqua */
1379  /* and their "day" pixels are 9:30 am, before daytime warming, so use amsr night field */
1380  ad = 1 - ad;
1381  }
1382 
1383  /* locate LL position within reference grid */
1384  ii = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), nx + 1), 0);
1385  jj = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), ny + 1), 0);
1386 
1387  /* compute longitude and latitude of that grid element */
1388  xx = ii * dx - 180.0 - dx / 2;
1389  yy = jj * dy - 90.0 - dy / 2;
1390 
1391  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
1392  t = (lon - xx) / dx;
1393  u = (lat - yy) / dy;
1394 
1395  ftmp = 0.0;
1396  ntmp = 0;
1397  if (sstref[jj ][ii ][ad] > sstbad + 1) {
1398  ftmp += sstref[jj ][ii ][ad];
1399  ntmp++;
1400  }
1401  if (sstref[jj ][ii + 1][ad] > sstbad + 1) {
1402  ftmp += sstref[jj ][ii + 1][ad];
1403  ntmp++;
1404  }
1405  if (sstref[jj + 1][ii + 1][ad] > sstbad + 1) {
1406  ftmp += sstref[jj + 1][ii + 1][ad];
1407  ntmp++;
1408  }
1409  if (sstref[jj + 1][ii ][ad] > sstbad + 1) {
1410  ftmp += sstref[jj + 1][ii ][ad];
1411  ntmp++;
1412  }
1413  if (ntmp > 0) {
1414  ftmp /= ntmp;
1415  reftmp[0][0] = (sstref[jj ][ii ][ad] > sstbad + 1 ? sstref[jj ][ii ][ad] : ftmp);
1416  reftmp[0][1] = (sstref[jj ][ii + 1][ad] > sstbad + 1 ? sstref[jj ][ii + 1][ad] : ftmp);
1417  reftmp[1][1] = (sstref[jj + 1][ii + 1][ad] > sstbad + 1 ? sstref[jj + 1][ii + 1][ad] : ftmp);
1418  reftmp[1][0] = (sstref[jj + 1][ii ][ad] > sstbad + 1 ? sstref[jj + 1][ii ][ad] : ftmp);
1419 
1420  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
1421 
1422  } else
1423  sst = sstbad;
1424 
1425  return (sst);
1426 }
1427 
1428 /* ----------------------------------------------------------------------------------- */
1429 /* get_atsrday() - read and interpolate ATSR 0.1-deg daily netcdf files */
1430 /* */
1431 /* S. Walsh, RSMAS, Feb 2011. */
1432 /* ----------------------------------------------------------------------------------- */
1433 
1434 #define OI1NX 3600
1435 #define OI1NY 1800
1436 
1437 float get_atsrday(char *sstfile, float lon, float lat) {
1438  static int firstCall = 1;
1439  static int nx = OI1NX;
1440  static int ny = OI1NY;
1441  static int32_t slen = 0;
1442  static float *sref = NULL;
1443  static float dx = 360.0 / OI1NX;
1444  static float dy = 180.0 / OI1NY;
1445 
1446  float sst = sstbad;
1447  int i, j, jl, nl;
1448  int ntmp;
1449  float xx, yy;
1450  float t, u;
1451  float reftmp[2][2];
1452  float ftmp;
1453 
1454  if (firstCall) {
1455 
1456  char name [H4_MAX_NC_NAME] = "";
1457  char sdsname[H4_MAX_NC_NAME] = "";
1458  int ll, jj, jjl;
1459  int32 sd_id;
1460  int32 sds_id;
1461  int32 rank;
1462  int32 nt;
1463  int32 dims[H4_MAX_VAR_DIMS];
1464  int32 nattrs;
1465  int32 start[4] = {0, 0, 0, 0};
1466  int32 status;
1467  int32_t stlen = 0;
1468  float *stmp = NULL;
1469 
1470  firstCall = 0;
1471 
1472  printf("Loading Daily 0.1-deg ATSR SST reference from %s\n", sstfile);
1473  printf("\n");
1474 
1475  /* allocate arrays */
1476  slen = (OI1NY + 2) * (OI1NX + 2);
1477  if ((sref = (float *) malloc(slen * sizeof (float))) == NULL) {
1478  fprintf(stderr, "-E- %s line %d: Unable to allocate sstref array.\n",
1479  __FILE__, __LINE__);
1480  exit(1);
1481  }
1482 
1483  stlen = OI1NY * OI1NX;
1484  if ((stmp = (float *) malloc(stlen * sizeof (float))) == NULL) {
1485  fprintf(stderr, "-E- %s line %d: Unable to allocate ssttmp array.\n",
1486  __FILE__, __LINE__);
1487  exit(1);
1488  }
1489 
1490  /* Open the file */
1491  sd_id = SDstart(sstfile, DFACC_RDONLY);
1492  if (sd_id == FAIL) {
1493  fprintf(stderr, "-E- %s line %d: error reading %s.\n",
1494  __FILE__, __LINE__, sstfile);
1495  exit(1);
1496  }
1497 
1498  strcpy(sdsname, "sst_skin");
1499  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1500  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1501  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) stmp);
1502  if (status != 0) {
1503  fprintf(stderr, "-E- %s Line %d: Error reading SDS %s from %s.\n",
1504  __FILE__, __LINE__, sdsname, sstfile);
1505  exit(1);
1506  }
1507  status = SDendaccess(sds_id);
1508  status = SDend(sd_id);
1509 
1510  /* ATSR grid is lon centers -179.95..179.95 in 1..3600; lat centers 89.95..-89.95 */
1511  /* first pix represents -180..-179.9,90..89.9 */
1512  /* flip and add wrapping border to simplify interpolation */
1513  /* new grid is -180.1,180.1 in i=0,3601 and -90.1,90.1 in j=0,1801 */
1514 
1515  for (j = 0; j < ny; j++) {
1516  /* stmp points to ny x nx array read in from file */
1517  jl = j * (nx + 2); /* offset to line j in stmp */
1518  /* sref points to ny+2 * nx+2 array */
1519  /* flip so need to calculate line in sref */
1520  /* stmp line 0 goes into sref line ny */
1521  /* stmp line ny-1 goes into sref line 1 */
1522  jj = ny - j;
1523  jjl = jj * (nx + 2); /* offset to line jj in sref */
1524  for (i = 0; i < nx; i++) {
1525  if (stmp[jl + i] > -998.0) {
1526  sref[jjl + i + 1] = stmp[jl + i];
1527  } else {
1528  sref[jjl + i + 1] = sstbad;
1529  }
1530  }
1531  sref[jl] = sref[jl + nx];
1532  sref[jl + nx + 2 + nx + 1] = sref[jl + nx + 2 + 1];
1533  }
1534  ll = (ny + 1)*(nx + 2); /* start of last line */
1535  nl = ny * (nx + 2); /* start of next to last line */
1536  for (i = 0; i < nx + 2; i++) {
1537  sref[i] = sref[nx + 2 + i];
1538  sref[ll + i] = sref[nl + i];
1539  }
1540  free(stmp);
1541  }
1542 
1543 
1544  /* locate LL position within reference grid */
1545  i = MAX(MIN((int) ((lon + 180.0 + dx / 2) / dx), OI1NX + 1), 0);
1546  j = MAX(MIN((int) ((lat + 90.0 + dy / 2) / dy), OI1NY + 1), 0);
1547 
1548  /* compute longitude and latitude of that grid element */
1549  xx = i * dx - 180.0 - dx / 2;
1550  yy = j * dy - 90.0 - dy / 2;
1551 
1552  /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
1553  t = (lon - xx) / dx;
1554  u = (lat - yy) / dy;
1555 
1556  ftmp = 0.0;
1557  ntmp = 0;
1558  jl = (j - 1)*(nx + 2); /* offset to start of line j */
1559  nl = j * (nx + 2); /* offset to start of line j+1 */
1560  if (sref[jl + i ] > sstbad + 1) {
1561  ftmp += sref[jl + i ];
1562  ntmp++;
1563  }
1564  if (sref[jl + i + 1] > sstbad + 1) {
1565  ftmp += sref[jl + i + 1];
1566  ntmp++;
1567  }
1568  if (sref[nl + i + 1] > sstbad + 1) {
1569  ftmp += sref[nl + i + 1];
1570  ntmp++;
1571  }
1572  if (sref[nl + i ] > sstbad + 1) {
1573  ftmp += sref[nl + i ];
1574  ntmp++;
1575  }
1576  if (ntmp > 0) {
1577  ftmp /= ntmp;
1578  reftmp[0][0] = (sref[jl + i ] > sstbad + 1 ? sref[jl + i ] : ftmp);
1579  reftmp[0][1] = (sref[jl + i + 1] > sstbad + 1 ? sref[jl + i + 1] : ftmp);
1580  reftmp[1][1] = (sref[nl + i + 1] > sstbad + 1 ? sref[nl + i + 1] : ftmp);
1581  reftmp[1][0] = (sref[nl + i ] > sstbad + 1 ? sref[nl + i ] : ftmp);
1582 
1583  sst = (1 - t)*(1 - u) * reftmp[0][0] + t * (1 - u) * reftmp[0][1] + t * u * reftmp[1][1] + (1 - t) * u * reftmp[1][0];
1584 
1585  } else
1586  sst = sstbad;
1587 
1588  return (sst);
1589 }
1590 
1591 /* ------------------------------------------------------------------------------------ */
1592 /* get_sstref() - retrieves reference sea surface temperature */
1593 /* ssttype has to be set for each pixel because climatology could be used */
1594 /* instead of sstfile */
1595 /* */
1596 /* B. Franz, SAIC, May 2004. */
1597 /* ---------------------------------------------------------------------------------------------- */
1598 #include "smi_climatology.h"
1599 #include "l1.h"
1600 
1601 float get_sstref(short reftyp, char *sstfile, l1str *l1rec, int32_t ip) {
1602  float lon = l1rec->lon[ip];
1603  float lat = l1rec->lat[ip];
1604  float solz = l1rec->solz[ip];
1605  int32_t xsatid = l1rec->l1file->subsensorID;
1606  int32_t sensorID = l1rec->l1file->sensorID;
1607  int16_t yr, dy;
1608  double sec;
1609 
1610  unix2yds(l1rec->scantime, &yr, &dy, &sec);
1611  int32_t year = (int32_t) yr;
1612  int32_t day = (int32_t) dy;
1613  float sst = sstbad;
1614  int32_t sd_id;
1615 
1616  if (format < 0) {
1617 
1618  /* Does the file exist? */
1619  if (access(sstfile, F_OK) || access(sstfile, R_OK)) {
1620  printf("-E- %s line %d: SST input file '%s' does not exist or cannot open.\n",
1621  __FILE__, __LINE__, sstfile);
1622  exit(1);
1623  }
1624 
1625  /* What is it? */
1626  if (reftyp == 8) {
1627  /* It's a 3-day Day or Night WindSat file from avgwindsat */
1628  format = WINDSAT3DN;
1629 
1630  } else if (reftyp == 7) {
1631  /* It's a 3-day WindSat file from SSMI */
1632  format = WINDSAT3DAY;
1633 
1634  } else if (reftyp == 6) {
1635  /* It's a daily WindSat file from SSMI */
1636  format = WINDSATDAY;
1637 
1638  } else if (reftyp == 5) {
1639  /* It's a 3-day Day or Night AMSR-E file from SSMI */
1640  format = AMSRE3DN;
1641 
1642  } else if (reftyp == 4) {
1643  /* It's a Daily NTEV2 file from Reynolds */
1644  format = NTEV2;
1645 
1646  } else if (reftyp == 3) {
1647  /* It's a Monthly ATSR file from Reynolds */
1648  format = ATSR;
1649 
1650  } else if (reftyp == 2) {
1651  /* It's a 3-day AMSR-E file from SSMI */
1652  format = AMSRE3DAY;
1653 
1654  } else if (reftyp == 1) {
1655  /* It's a daily AMSR-E file from SSMI */
1656  format = AMSREDAY;
1657 
1658  } else {
1659 
1660  char title[255] = "";
1661 
1662  // SDstart seg faults on the mac, so need to protect the
1663  // call with Hishdf()
1664  if (Hishdf(sstfile))
1665  sd_id = SDstart(sstfile, DFACC_RDONLY);
1666  else
1667  sd_id = FAIL;
1668 
1669  if (sd_id != FAIL) {
1670  /* Format is HDF-like */
1671  char title[255] = "";
1672  if (SDreadattr(sd_id, SDfindattr(sd_id, "title"), (VOIDP) title) == 0) {
1673  if (strstr(title, "Daily-OI-V2") != NULL) {
1674  format = OISSTV2D;
1675  } else if (strstr(title, "ARC Sea Surface Temperature") != NULL) {
1676  format = ATSRDAY;
1677 
1678  } else {
1679  printf("-E- %s line %d: Unable to initialize SST file\n", __FILE__, __LINE__);
1680  printf("%s %d\n", sstfile, day);
1681  exit(1);
1682  }
1683  } else if (SDreadattr(sd_id, SDfindattr(sd_id, "Title"), (VOIDP) title) == 0) {
1684  if (strstr(title, "SST Climatology") != NULL) {
1685  format = PATHCLIM;
1686  if (smi_climatology_init(sstfile, day, SST) != 0) {
1687  printf("-E- %s line %d: Unable to initialize SST file\n", __FILE__, __LINE__);
1688  printf("%s %d\n", sstfile, day);
1689  exit(1);
1690  }
1691  }
1692  }
1693  SDend(sd_id);
1694 
1695  } else {
1696  if (nc_open(sstfile, NC_NOWRITE, &sd_id) == NC_NOERR) {
1697  /* Format is netCDF */
1698  if (nc_get_att_text(sd_id, NC_GLOBAL, "title", title) == NC_NOERR) {
1699  if (strstr(title, "Daily-OI") ||
1700  strstr(title, "NOAA/NCEI 1/4 Degree Daily Optimum Interpolation Sea Surface Temperature (OISST) Analysis")) {
1701  format = OISSTV2D;
1702  } else if (strstr(title, "CMC") != NULL) {
1703  format = CMC;
1704  } else if (strstr(title, "Standard Mapped Image") != NULL) {
1705  format = L3M;
1706  }else {
1707  printf("-E- %s : Unable to initialize SST file\n", __FILE__);
1708  printf("%s %d\n", sstfile, day);
1709  exit(1);
1710  }
1711  }
1712  nc_close(sd_id);
1713  } else {
1714 
1715  /* Format is not HDF or netCDF. Assuming OISST Binary. */
1716  format = OISSTBIN;
1717  }
1718  }
1719  }
1720 
1721  /* need a back-up climatology for NRT products */
1722  if (format != PATHCLIM) {
1723  char *tmp_str;
1724  char file [FILENAME_MAX] = "";
1725  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
1726  printf("OCDATAROOT environment variable is not defined.\n");
1727  exit(1);
1728  }
1729  strcpy(file, tmp_str);
1730  strcat(file, "/common/sst_climatology.hdf");
1731  if (smi_climatology_init(file, day, SST) != 0) {
1732  printf("-E- %s line %d: Unable to initialize SST file\n", __FILE__, __LINE__);
1733  printf("%s %d\n", file, day);
1734  exit(1);
1735  }
1736  }
1737  }
1738 
1739  switch (format) {
1740  case PATHCLIM:
1741  sst = smi_climatology(lon, lat, SST);
1742  l1rec->ssttype[ip] = 0; /* SST Climatology */
1743  /* don't know if this is skin or bulk, but sst doesn't use these so don't care */
1744  break;
1745  case OISSTBIN:
1746  sst = get_oisst(sstfile, lon, lat);
1747  l1rec->ssttype[ip] = 1; /* Reynolds Weekly 1 deg OI */
1748  /* All Reynolds reference fields are bulk, sst calculation needs skin */
1749  sst = sst - 0.17;
1750  break;
1751  case OISSTV2D:
1752  sst = get_oisstv2d(sstfile, lon, lat);
1753  l1rec->ssttype[ip] = 2; /* Reynolds Daily 0.25-deg OI */
1754  /* All Reynolds reference fields are bulk, sst calculation needs skin */
1755  sst = sst - 0.17;
1756  break;
1757  case AMSRE3DAY:
1758  sst = get_amsre(sstfile, lon, lat, solz, xsatid, sensorID);
1759  l1rec->ssttype[ip] = 3; /* AMSR-E 3-Day 0.25-deg */
1760  /* AMSR is microwave skin, needs correction to skin */
1761  sst = sst - 0.15;
1762  break;
1763  case AMSREDAY:
1764  sst = get_amsre(sstfile, lon, lat, solz, xsatid, sensorID);
1765  l1rec->ssttype[ip] = 4; /* AMSR-E Daily 0.25-deg */
1766  /* AMSR is microwave skin, needs correction to skin */
1767  sst = sst - 0.15;
1768  break;
1769  case ATSR:
1770  sst = get_atsr(sstfile, lon, lat, year, day, 1);
1771  l1rec->ssttype[ip] = 5; /* ATSR Monthly 0.25-deg */
1772  /* atsr is skin, no correction needed */
1773  break;
1774  case NTEV2:
1775  sst = get_ntev2(sstfile, lon, lat, year, day);
1776  l1rec->ssttype[ip] = 6; /* NTEV2 Daily 0.25-deg */
1777  /* All Reynolds reference fields are bulk, sst calculation needs skin */
1778  sst = sst - 0.17;
1779  break;
1780  case AMSRE3DN:
1781  sst = get_amsre(sstfile, lon, lat, solz, xsatid, sensorID);
1782  l1rec->ssttype[ip] = 7; /* AMSR-E 3-Day Day or Night 0.25-deg */
1783  /* AMSR is microwave skin, needs correction to skin */
1784  if (sst > sstbad + 1) {
1785  sst = sst - 0.15;
1786  }
1787  break;
1788  case ATSRDAY:
1789  sst = get_atsrday(sstfile, lon, lat);
1790  l1rec->ssttype[ip] = 8; /* ATSR Daily 0.1-deg */
1791  /* We use the sst_skin value from the ATSR files */
1792  /* convert Kelvin to C */
1793  sst = sst - 273.15;
1794  break;
1795  case WINDSAT3DAY:
1796  sst = get_windsat(sstfile, lon, lat, solz, xsatid, sensorID);
1797  l1rec->ssttype[ip] = 9; /* WindSat 3-Day 0.25-deg */
1798  /* WindSat is microwave skin, needs correction to skin */
1799  sst = sst - 0.15;
1800  break;
1801  case WINDSATDAY:
1802  sst = get_windsat(sstfile, lon, lat, solz, xsatid, sensorID);
1803  l1rec->ssttype[ip] = 10; /* WindSat Daily 0.25-deg */
1804  /* WindSat is microwave skin, needs correction to skin */
1805  sst = sst - 0.15;
1806  break;
1807  case WINDSAT3DN:
1808  sst = get_windsat(sstfile, lon, lat, solz, xsatid, sensorID);
1809  l1rec->ssttype[ip] = 11; /* WindSat 3-Day Day or Night 0.25-deg */
1810  /* WindSat is microwave skin, needs correction to skin */
1811  if (sst > sstbad + 1) {
1812  sst = sst - 0.15;
1813  }
1814  break;
1815  case CMC:
1816  sst = get_cmcsst(sstfile, lon, lat, 1); // units of Kelvin
1817  l1rec->ssttype[ip] = 12; /* CMC analyzed SST */
1818  /* CMC analyzed SST fields are "foundation" i.e. bulk, sst calculation needs skin correction*/
1819  sst = sst - 0.17;
1820  break;
1821  case L3M:
1822  sst = get_cmcsst(sstfile, lon, lat, 0);// units of degress_C
1823  l1rec->ssttype[ip] = 13; /* L3M SST */
1824  break;
1825  default:
1826  printf("-E- %s line %d: unknown SST input file format for %s.\n", __FILE__, __LINE__, sstfile);
1827  break;
1828  }
1829 
1830  if (sst < (sstbad + 1)) {
1831  sst = smi_climatology(lon, lat, SST);
1832  l1rec->ssttype[ip] = 0; /* SST Climatology */
1833  /* don't know if this is skin or bulk, but sst doesn't use these so don't care */
1834  }
1835 
1836  return (sst);
1837 }
OINX
#define OINX
Definition: sstref.c:453
MAX
#define MAX(A, B)
Definition: swl0_utils.h:25
WINDSATDAY
#define WINDSATDAY
Definition: sstref.c:29
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
eday
int16 eday
Definition: l1_czcs_hdf.c:17
smi_climatology_init
int smi_climatology_init(char *file, int day, int prodID)
Definition: smi_climatology.c:13
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
NO06
#define NO06
Definition: l1_aci.h:5
nl
int32_t nl
Definition: atrem_corl1.h:132
AVHRR
#define AVHRR
Definition: sensorDefs.h:15
OINY
#define OINY
Definition: sstref.c:454
CMC
#define CMC
Definition: sstref.c:31
get_windsat
float get_windsat(char *sstfile, float lon, float lat, float solz, int32_t xsatid, int32_t sensorID)
Definition: sstref.c:1208
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
NULL
#define NULL
Definition: decode_rs.h:63
get_cmcsst
float get_cmcsst(char *sstfile, float lon, float lat, int isKelvin)
Definition: sstref.c:242
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
AMSRE3DAY
#define AMSRE3DAY
Definition: sstref.c:22
OI4NY
#define OI4NY
Definition: sstref.c:1206
TRUE
#define TRUE
Definition: rice.h:165
MODIST
#define MODIST
Definition: sensorDefs.h:18
get_sstref
float get_sstref(short reftyp, char *sstfile, l1str *l1rec, int32_t ip)
Definition: sstref.c:1601
WINDSAT3DAY
#define WINDSAT3DAY
Definition: sstref.c:28
eyear
int16 eyear
Definition: l1_czcs_hdf.c:17
syear
int syear
Definition: l1_czcs_hdf.c:15
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
free2d_short
void free2d_short(short **p)
Free a two-dimensional array created by allocate2d_short.
Definition: allocate2d.c:96
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
endianess
int endianess(void)
determine endianess
Definition: endianess.c:10
swapc_bytes
int swapc_bytes(char *in, int nbyte, int ntime)
Definition: swapc_bytes.c:4
get_oisstv2d
float get_oisstv2d(char *sstfile, float lon, float lat)
Definition: sstref.c:52
get_atsrday
float get_atsrday(char *sstfile, float lon, float lat)
Definition: sstref.c:1437
l12_proto.h
sday
int sday
Definition: l1_czcs_hdf.c:15
PATHCLIM
#define PATHCLIM
Definition: sstref.c:19
OI1NY
#define OI1NY
Definition: sstref.c:1435
run_local.lat
lat
Definition: run_local.py:69
ATSR
#define ATSR
Definition: sstref.c:24
run_local.lon
lon
Definition: run_local.py:69
OISSTV2D
#define OISSTV2D
Definition: sstref.c:21
get_oisst
float get_oisst(char *sstfile, float lon, float lat)
Definition: sstref.c:456
NO12
#define NO12
Definition: l1_aci.h:11
ATSRDAY
#define ATSRDAY
Definition: sstref.c:27
NTEV2
#define NTEV2
Definition: sstref.c:25
CtoK
#define CtoK
Definition: sstref.c:34
slope
float32 slope[]
Definition: l2lists.h:30
NO17
#define NO17
Definition: l1_aci.h:15
get_atsr
float get_atsr(char *sstfile, float lon, float lat, int32_t year, int32_t day, int32_t minatsrcnt)
Definition: sstref.c:744
getHDFattr
int getHDFattr(int32_t fileID, const char attrname[], const char sdsname[], void *data)
NO08
#define NO08
Definition: l1_aci.h:7
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
l12_parms.h
isleap
int isleap(int year)
Definition: isleap.c:3
smi_climatology.h
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
WINDSAT3DN
#define WINDSAT3DN
Definition: sstref.c:30
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
smi_climatology
float smi_climatology(float lon, float lat, int prodID)
Definition: smi_climatology.c:126
u
data_t u
Definition: decode_rs.h:74
OISSTBIN
#define OISSTBIN
Definition: sstref.c:20
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
AMSREDAY
#define AMSREDAY
Definition: sstref.c:23
AMSRE3DN
#define AMSRE3DN
Definition: sstref.c:26
L3M
#define L3M
Definition: sstref.c:32
allocate2d_short
short ** allocate2d_short(size_t h, size_t w)
Allocate a two-dimensional array of type short of a given size.
Definition: allocate2d.c:79
SST
The driver class to interface with the pure legacy C code.
Definition: sst_cloud_mask.cpp:1826
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
get_ntev2
float get_ntev2(char *sstfile, float lon, float lat, int32_t year, int32_t day)
Definition: sstref.c:592
offset
l2prod offset
Definition: HOWTO_Add_a_product.txt:70
leap
logical function leap(YEAR)
Definition: leap.f:10
i
int i
Definition: decode_rs.h:71
get_amsre
float get_amsre(char *sstfile, float lon, float lat, float solz, int32_t xsatid, int32_t sensorID)
Definition: sstref.c:1016
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:91
OI4NX
#define OI4NX
Definition: sstref.c:1205
l1_aci.h
NO10
#define NO10
Definition: l1_aci.h:9
MDN.setup.version
version
Definition: setup.py:15
NO15
#define NO15
Definition: l1_aci.h:13
OI1NX
#define OI1NX
Definition: sstref.c:1434


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