ocssw V2020
sst.c
Go to the documentation of this file.
1 #include <sys/types.h>
2 #include <unistd.h>
3 
4 /* ============================================================================ */
5 /* sst.c - functions for retrieval of sea surface temperature */
6 /* */
7 /* Written By: B. Franz, NASA/SIMBIOS, August 2003. */
8 /* */
9 /* ============================================================================ */
10 
11 /*
12  * nightcube input parameter eliminated - hardcoded to 1 in calls to the sses
13  * functions. This is due to the note in comp_sses_sst:
14  * "GHRSST sses defintion is that the sses are based on night stats only"
15  * Since the SSES products are only for GHRSST, silly to make it an option...
16  * But July 2017, GHRSST definition has changed and viirs now has day and night sses
17  * so hardcoding it in (for now? Why was there an option before?)
18  */
19 #include "l12_proto.h"
20 #include "flags_sst.h"
21 #include "d3940tref.h"
22 #include "l1_aci_hdf.h"
23 #include "sst_dsdi.h"
24 
25 static float sstbad = BAD_FLT;
26 static int32 evalmask = 0;
27 
28 
29 typedef struct stat_struct {
30  float min;
31  float max;
32  float avg;
33  float med;
34  float sqr;
35  float sd;
36  float cen;
37  int cnt;
38 } statstr;
39 
40 #define CtoK 273.15 /* conversion between degree C and Kelvin*/
41 
42 /* sst sensor-specific error table (SSES) v6 viirs */
43 #define NSSTMAXv6v 7 /* was 8 */
44 #define NDAYMAXv6v 2 /* viirs sst now has day and night sses */
45 #define NDAYMAXv6v3 1 /* viirs sst3 is night only */
46 #define NQUARMAXv6v 4
47 #define NSENZMAXv6v 4
48 #define NDIFFMAXv6v 4
49 #define NDIFFMAXv6v3 4
50 #define NLATMAXv6v 6 /* was 5 */
51 #define NQUALMAXv6v 5 /* was 4, then 5 */
52 #define NSSESDIMv6v 7
53 /* sst sensor-specific error table (SSES) v6 modis */
54 #define NSSTMAXv6m 7 /* was 8 */
55 #define NDAYMAXv6m 2 /* sst now has day and night sses */
56 #define NDAYMAXv6m4 1 /* sst4 is night only */
57 #define NQUARMAXv6m 4
58 #define NSENZMAXv6m 4
59 #define NDIFFMAXv6m 4
60 #define NDIFFMAXv6m4 4
61 #define NLATMAXv6m 7 /* was 5 */
62 #define NQUALMAXv6m 5 /* was 4, then 5 */
63 #define NSSESDIMv6m 7
64 /* sst sensor-specific error table (SSES) v6 avhrr */
65 #define NSSTMAXv6a 7
66 #define NQUARMAXv6a 4
67 #define NSENZMAXv6a 4
68 #define NDIFFMAXv6a 4
69 #define NLATMAXv6a 6
70 #define NQUALMAXv6a 9
71 #define NSSESDIMv6a 6
72 /* avhrr sensor-specific error table (SSES)
73  * -2 to 3C, 3+ to 8C, 8+ to 13C, 13+ to 18C, 18+ to 23C, 23+ to 28C, > 28C
74  * set hdf limits to: -2, 3, 8, 13, 18, 23, 28
75  * Day, Night
76  * 1Q, 2Q, 3Q, 4Q
77  * 0 to 30 deg, 30+ to 40 deg, 40+ to 50 deg, 50+ deg
78  * set hdf limits to : 0, 30, 40, 50
79  * sst: < 0.0C, 0.0 to 0.7C, 0.7+ to 2.0C, > 2.0C
80  * sst: set hdf limits to -1000, 0.0, 0.7, 2.0
81  * sst4:< -0.5C, -0.5 to 0.0C, 0.0+ to 0.5C, > 0.5C
82  * sst4: set hdf limits to -1000, -0.5, 0.0, 0.5
83  * 90S to 40S, 40S+ to 20S, 20S+ to Eq, Eq+ to 20N, 20N+ to 40N, 40N+ to 90N
84  * set hdf lower limits to -90 -40, -20, 0, 20, 40
85  * avhrr: 0, 1, 2, 3, 4, 5, 6, 7, 8, modis: 0, 1, 2, 3, 4
86  * set hdf lower limits to avhrr: 0, 1, 2, 3, 4, modis: 0, 1, 2, 3, 4
87  * use bias=-10, stdv=5 for avhrr qual 8, and modis qual 4
88  */
89 
90 /* this v6 structure is for viirs */
91 
92 typedef struct ssestab_structv6v {
93  int nsst;
94  int nday;
95  int nquar;
96  int nsenz;
97  int ndiff;
98  int nlat;
99  int nqual;
100 
101  float sst[NSSTMAXv6v];
102  float senz[NSENZMAXv6v];
103  float diff[NDIFFMAXv6v];
104  float lat[NLATMAXv6v];
105  float qual[NQUALMAXv6v];
106 
107  float bias[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v][NSSTMAXv6v];
108  float stdv[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v][NSSTMAXv6v];
109  float bias_mean[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v][NSSTMAXv6v];
110  int16 counts[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v][NSSTMAXv6v];
111 
112 } ssestabstrv6v;
113 
114 /* structure for viirs sst3 stats, with different NDAY and NDIFF sizes */
115 typedef struct ssestab_structv6v3 {
116  int nsst;
117  int nday;
118  int nquar;
119  int nsenz;
120  int ndiff;
121  int nlat;
122  int nqual;
123 
124  float sst[NSSTMAXv6v];
125  float senz[NSENZMAXv6v];
126  float diff[NDIFFMAXv6v3];
127  float lat[NLATMAXv6v];
128  float qual[NQUALMAXv6v];
129 
130  float bias[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v3][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v3][NSSTMAXv6v];
131  float stdv[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v3][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v3][NSSTMAXv6v];
132  float bias_mean[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v3][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v3][NSSTMAXv6v];
133  int16 counts[NQUALMAXv6v][NLATMAXv6v][NDIFFMAXv6v3][NSENZMAXv6v][NQUARMAXv6v][NDAYMAXv6v3][NSSTMAXv6v];
134 
135 } ssestabstrv6v3;
136 
137 static ssestabstrv6v sses_sstv6v;
138 static ssestabstrv6v3 sses_sst3v6v3;
139 
140 /* this v6 structure is for modis */
141 
142 typedef struct ssestab_structv6m {
143  int nsst;
144  int nday;
145  int nquar;
146  int nsenz;
147  int ndiff;
148  int nlat;
149  int nqual;
150 
151  float sst[NSSTMAXv6m];
152  float senz[NSENZMAXv6m];
153  float diff[NDIFFMAXv6m];
154  float lat[NLATMAXv6m];
155  float qual[NQUALMAXv6m];
156 
157  float bias[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m][NSSTMAXv6m];
158  float stdv[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m][NSSTMAXv6m];
159  float bias_mean[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m][NSENZMAXv6m][NQUARMAXv6m][NSSTMAXv6m];
160  int16 counts[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m][NSENZMAXv6m][NQUARMAXv6m][NSSTMAXv6m];
161 
162 } ssestabstrv6m;
163 
164 /* structure for modis sst4 stats, with different NDAY and NDIFF sizes */
165 typedef struct ssestab_structv6m4 {
166  int nsst;
167  int nday;
168  int nquar;
169  int nsenz;
170  int ndiff;
171  int nlat;
172  int nqual;
173 
174  float sst[NSSTMAXv6m];
175  float senz[NSENZMAXv6m];
176  float diff[NDIFFMAXv6m4];
177  float lat[NLATMAXv6m];
178  float qual[NQUALMAXv6m];
179 
180  float bias[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m4][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m4][NSSTMAXv6m];
181  float stdv[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m4][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m4][NSSTMAXv6m];
182  float bias_mean[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m4][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m4][NSSTMAXv6m];
183  int16 counts[NQUALMAXv6m][NLATMAXv6m][NDIFFMAXv6m4][NSENZMAXv6m][NQUARMAXv6m][NDAYMAXv6m4][NSSTMAXv6m];
184 
185 } ssestabstrv6m4;
186 static ssestabstrv6m sses_sstv6m;
187 static ssestabstrv6m4 sses_sst4v6m;
188 
189 /* this sses structure is for avhhr which has more quality levels and modis and viirs */
190 typedef struct ssestab_structv6a {
191  int nsst;
192  int nquar;
193  int nsenz;
194  int ndiff;
195  int nlat;
196  int nqual;
197 
198  float sst[NSSTMAXv6a];
199  float senz[NSENZMAXv6a];
200  float diff[NDIFFMAXv6a];
201  float lat[NLATMAXv6a];
202  float qual[NQUALMAXv6a];
203 
204  float bias[NQUALMAXv6a][NLATMAXv6a][NDIFFMAXv6a][NSENZMAXv6a][NQUARMAXv6a][NSSTMAXv6a];
205  float stdv[NQUALMAXv6a][NLATMAXv6a][NDIFFMAXv6a][NSENZMAXv6a][NQUARMAXv6a][NSSTMAXv6a];
206  float bias_mean[NQUALMAXv6a][NLATMAXv6a][NDIFFMAXv6a][NSENZMAXv6a][NQUARMAXv6a][NSSTMAXv6a];
207  int16 counts[NQUALMAXv6a][NLATMAXv6a][NDIFFMAXv6a][NSENZMAXv6a][NQUARMAXv6a][NSSTMAXv6a];
208 
209 } ssestabstrv6a;
210 
211 static ssestabstrv6a sses_sstv6a;
212 
213 /* scans of computed quantities */
214 typedef float s_array[MAXPIX];
215 static s_array *sstq, *sst4q, *sst3q;
216 //static float *sst = NULL;
217 //static float *sst4 = NULL;
218 //static float *sst3 = NULL;
219 static float *treesum = NULL;
220 static int16 *flags_sst = NULL;
221 static int16 *flags_sst4 = NULL;
222 static int16 *flags_sst3 = NULL;
223 static int8 *qual_sst = NULL;
224 static int8 *qual_sst4 = NULL;
225 static int8 *qual_sst3 = NULL;
226 static float *bias_sst = NULL;
227 static float *bias_sst4 = NULL;
228 static float *bias_sst3 = NULL;
229 static float *stdv_sst = NULL;
230 static float *stdv_sst4 = NULL;
231 static float *stdv_sst3 = NULL;
232 static float *bias_mean_sst = NULL;
233 static float *bias_mean_sst4 = NULL;
234 static float *bias_mean_sst3 = NULL;
235 static int16 *bias_counts_sst = NULL;
236 static int16 *bias_counts_sst4 = NULL;
237 static int16 *bias_counts_sst3 = NULL;
238 static float *dsdi_correction = NULL;
239 
240 static float *d3940ref = NULL;
241 
242 // set up arrays for box
243 static float *LtRED_maxmin = NULL;
244 static float *Bt11_maxmin = NULL;
245 static float *Bt11_max = NULL;
246 static float *Bt11_min = NULL;
247 static float *Bt11_stdev = NULL;
248 static float *Bt12_maxmin = NULL;
249 static float *Bt12_min = NULL;
250 static float *Bt37_maxmin = NULL;
251 static float *Bt37_stdev = NULL;
252 static float *Bt39_maxmin = NULL;
253 static float *Bt40_maxmin = NULL;
254 //static float *Bt40_avg = NULL;
255 static float *Bt40_stdev = NULL;
256 static float *Bt85_min = NULL;
257 static float *Bt73_max = NULL;
258 static float *rhoCirrus_maxmin = NULL;
259 static float *rhoCirrus_min = NULL;
260 static float *rhoCirrus_max = NULL;
261 static float *rhotRED_maxmin = NULL;
262 static float *rhotRED_min = NULL;
263 static float *rhotRED_max = NULL;
264 static float *rhotNIR7_min = NULL;
265 static float *rhot16_min = NULL;
266 static float *rhotRED = NULL;
267 static float *rhotNIR7 = NULL;
268 static float *rhot16 = NULL;
269 static float *sst_stdev = NULL;
270 
271 /* precomputed indices */
272 static int32_t recnumSST = -1;
273 static int haveSST4 = 0;
274 static int haveSST = 0;
275 static int haveSSES = 1;
276 static int haveRed = 0;
277 static int ib07 = -1;
278 static int ib08 = -1;
279 static int ib16 = -1;
280 static int ib37 = -1;
281 static int ib39 = -1;
282 static int ib40 = -1;
283 static int ib67 = -1;
284 static int ib73 = -1;
285 static int ib85 = -1;
286 static int ib11 = -1;
287 static int ib12 = -1;
288 static int ibred = -1;
289 static int nbvis = -1;
290 static int nbir = -1;
291 static float satred;
292 static int isV5 = 0;
293 
294 /* quality test parameters */
295 static int fullscanpix = 1354; // intialize to modis, will get set in init_sst for others
296 static int32_t cldbox = 3;
297 static int32_t cldboxv = 5;
298 static int sstboxcscan = -1;
299 static float cldthresh = 0.01; /* modis */
300 static float cldthreshv = 0.04; /* viirs */
301 /* Sue read the comment before changing btbox for aqua !!!!!!!!*/
302 static int32_t btbox = 3; /* SUE: read the next comment!!!!! */
303 /* if btbox changes to 5 for aqua then change modisa/msl12_filter.dat
304  * to keep 7 lines for btavg (Bt40 detector zero replacement) to work
305  */
306 static int32_t btboxv = 5;
307 static int32_t csstbox = -1;
308 static float hisenz = 55.0;
309 static float hisenza = 45.0;
310 static float vhisenz = 75.0;
311 static float vhisenza = 55.0; /* make this higher? ask Bob */
312 static float vhisenzv2 = 65.0; /* for VIIRS v6.4 sst2b May 2015 */
313 static float solznight = SOLZNIGHT; /* becomes SOLZNIGHTA for AVHRR */
314 static float Btmin = -4.0;
315 static float Btmina = -10.0;
316 //static float Btminv = -4.0;//+CtoK;
317 static float Btmax = 37.0; /* pre Nov 2012 was 33.0 */
318 static float Btmaxa = 37.0; /* pre Nov 2012 was 35.0 */
319 //static float Btmaxv = 37.0;//+CtoK; /* pre Nov 2012 was 33.0+CtoK */
320 static float Btmax40 = 35.0;
321 //static float Btmax40v = 35.0;//+CtoK;
322 static float SSTmin = -1.8; /* Peter Minnett, 2017: should be -1.8 instead of -2.0 */
323 static float SSTmax = 40.0; /* pre Nov 2012 was 45.0 */
324 static float SSTmaxa = 40.0; /* pre Nov 2012 was 35.0 */
325 static float SSTmaxn = 37.0; /* night sst max */
326 static float glintmax = 0.005;
327 static float dBtmin = 0.0;
328 static float dBtmax = 3.6;
329 static float dBt4min = 0.0;
330 static float dBt4max = 8.0;
331 static float SSTdiffa = 2.0;
332 static float SSTdiff = 3.0;
333 /* Sue: test Nov 2016 - use sstvdiff instead of sstdiff */
334 /*
335  * Kay discovered an issue...
336  * 1. The two thresholds, SST4diff1,and SST4diff2, should be negative, -0.8, and -1.0,
337 respectively
338  * 2. The comparisons should be less than eg. dSST< SST4diff1 and dSST < SSTdiff2
339  */
340 //static float SSTdiff = 5.0;
341 static float SSTvdiff = 5.0;
342 static float SST4diff1 = -0.8;
343 static float SST4diff2 = -1.0;
344 static float SST3diff1 = 0.8;
345 static float SST3diff2 = 1.0;
346 /* Bt11unif1 changes in run_avhrr_sst because it varies by satellite */
347 static float Bt11unif1 = 0.7;
348 static float Bt12unif1 = 0.7;
349 static float Bt11unif2 = 1.2;
350 static float Bt12unif2 = 1.2;
351 static float Bt37unif1 = 0.7;
352 static float Bt37unif2 = 1.2;
353 static float Bt39unif1 = 0.7;
354 static float Bt40unif1 = 0.7;
355 static float Bt39unif2 = 1.2;
356 static float Bt40unif2 = 1.2;
357 static float dBtrefmin = -1.1;
358 static float dBtrefmax = 10.0;
359 /* tighter test between 10S and 30N and -105 to 105 longitude */
360 static float equatorialNorth = 30.0;
361 static float equatorialSouth = -10.0;
362 static float equatorialWest = -105.0;
363 static float equatorialEast = 105.0;
364 
365 /* NOTE: flag bit settings are now in flags_sst.h */
366 /* flag bit settings */
367 
368 static float latwin = 2.5; /* half the overlap for lat band coeffs */
369 
370 static int32 tmonth = -1;
371 
372 static int StartOfMonth[2][12] = {
373  { 0, 31, 59, 90, 120, 151, 181, 212, 243,
374  273, 304, 334},
375  { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}
376 };
377 
378 /* ----------------------------------------------------------------------------------- */
379 /* btrefdiffv6() - scan-dependent test on 4um BT temps to reference (RSMAS) */
380 /* */
381 /* B. Franz, SAIC, August 2005. */
382 
383 /* ----------------------------------------------------------------------------------- */
384 float btrefdiffv6(int32_t ip, float BT39, float BT40, l2str *l2rec) {
385  /* Gui's new table is by senz, which goes from - to + */
386  /* sza and trefv6 are in the d3940tref.h include file */
387 
388  float diff;
389  float satzdir;
390  float senz;
391  float tref;
392 
393  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
394  /* viirs has 3200 pixels per scan, pixnum's start at zero
395  * so we want pixels 0..1599 in the first half and 1600..3199 in the 2nd */
396  satzdir = (l2rec->l1rec->pixnum[ip] < fullscanpix / 2) ? -1.0 : 1.0;
397  senz = l2rec->l1rec->senz[ip] * satzdir;
398  tref = linterp(sza, trefv6, NSENZ, senz);
399  diff = BT39 - BT40 - tref;
400 
401  return (diff);
402 }
403 
404 /* ----------------------------------------------------------------------------------------------- */
405 /* load_sses_sstv6v() - loads the specified SSES (sensor-specific error stats) table for VIIRS sst */
406 
407 /* ----------------------------------------------------------------------------------------------- */
408 void load_sses_sstv6v(int32_t sensorID, char *file, ssestabstrv6v *sses) {
409  int32 sd_id;
410  int32 sds_id;
411  int32 status;
412  int32 rank;
413  int32 nt;
414  int32 nattrs;
415  int32 dims[NSSESDIMv6v];
416  int32 start[NSSESDIMv6v];
417 
418  char name[H4_MAX_NC_NAME] = "";
419  char sdsname[H4_MAX_NC_NAME] = "";
420 
421  memset(start, 0, NSSESDIMv6v * sizeof (int32));
422 
423  if (strcmp(file, "") == 0) {
424  printf("\nNo SSES data provided for this sensor.\n");
425  haveSSES = 0;
426  return;
427  }
428 
429  printf("\nLoading SSES table from %s\n", file);
430 
431  /* open the file and initiate the SD interface */
432  sd_id = SDstart(file, DFACC_RDONLY);
433  if (sd_id == -1) {
434  printf("-E- %s line %d: Error opening file %s.\n", __FILE__, __LINE__,
435  file);
436  exit(1);
437  }
438 
439  /* read the bias and standard deviation */
440 
441  strcpy(sdsname, "bias");
442  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
443  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
444  if (rank != NSSESDIMv6v) {
445  printf(
446  "-E- %s line %d: Table dimensions for %s do not match expectation. got: %d expected %d\n",
447  __FILE__, __LINE__, sdsname, rank, NSSESDIMv6v);
448  exit(1);
449  }
450  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
451  || dims[2] != NDIFFMAXv6v || dims[3] != NSENZMAXv6v
452  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v
453  || dims[6] != NSSTMAXv6v) {
454  printf(
455  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
456  __FILE__, __LINE__, sdsname);
457  exit(1);
458  }
459  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
460  if (status != 0) {
461  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
462  __LINE__, sdsname, file);
463  exit(1);
464  }
465  status = SDendaccess(sds_id);
466 
467  strcpy(sdsname, "stdv");
468  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
469  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
470  if (rank != NSSESDIMv6v) {
471  printf(
472  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
473  __FILE__, __LINE__, sdsname);
474  exit(1);
475  }
476  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
477  || dims[2] != NDIFFMAXv6v || dims[3] != NSENZMAXv6v
478  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v
479  || dims[6] != NSSTMAXv6v) {
480  printf(
481  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
482  __FILE__, __LINE__, sdsname);
483  exit(1);
484  }
485  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
486  if (status != 0) {
487  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
488  __LINE__, sdsname, file);
489  exit(1);
490  }
491  status = SDendaccess(sds_id);
492 
493  /* new hypercubes have bias_mean and counts also */
494 
495  strcpy(sdsname, "bias_mean");
496  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
497  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
498  if (rank != NSSESDIMv6v) {
499  printf(
500  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
501  __FILE__, __LINE__, sdsname);
502  exit(1);
503  }
504  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
505  || dims[2] != NDIFFMAXv6v || dims[3] != NSENZMAXv6v
506  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v
507  || dims[6] != NSSTMAXv6v) {
508  printf(
509  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
510  __FILE__, __LINE__, sdsname);
511  exit(1);
512  }
513  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias_mean);
514  if (status != 0) {
515  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
516  __LINE__, sdsname, file);
517  exit(1);
518  }
519  status = SDendaccess(sds_id);
520 
521  strcpy(sdsname, "counts");
522  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
523  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
524  if (rank != NSSESDIMv6v) {
525  printf(
526  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
527  __FILE__, __LINE__, sdsname);
528  exit(1);
529  }
530  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
531  || dims[2] != NDIFFMAXv6v || dims[3] != NSENZMAXv6v
532  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v
533  || dims[6] != NSSTMAXv6v) {
534  printf(
535  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
536  __FILE__, __LINE__, sdsname);
537  exit(1);
538  }
539  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->counts);
540  if (status != 0) {
541  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
542  __LINE__, sdsname, file);
543  exit(1);
544  }
545  status = SDendaccess(sds_id);
546 
547  /* save the table dimensions */
548  sses->nqual = dims[0];
549  sses->nlat = dims[1];
550  sses->ndiff = dims[2];
551  sses->nsenz = dims[3];
552  sses->nquar = dims[4];
553  sses->nday = dims[5];
554  sses->nsst = dims[6];
555 
556  /* read the table indice ranges */
557 
558  strcpy(sdsname, "sst");
559  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
560  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
561  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
562  if (status != 0) {
563  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
564  __LINE__, sdsname, file);
565  exit(1);
566  }
567  status = SDendaccess(sds_id);
568  sses->nsst = dims[0];
569 
570  strcpy(sdsname, "senz");
571  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
572  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
573  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
574  if (status != 0) {
575  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
576  __LINE__, sdsname, file);
577  exit(1);
578  }
579  status = SDendaccess(sds_id);
580  sses->nsenz = dims[0];
581 
582  strcpy(sdsname, "BTdiff");
583  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
584  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
585  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
586  if (status != 0) {
587  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
588  __LINE__, sdsname, file);
589  exit(1);
590  }
591  status = SDendaccess(sds_id);
592  sses->ndiff = dims[0];
593 
594  strcpy(sdsname, "lat");
595  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
596  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
597  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
598  if (status != 0) {
599  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
600  __LINE__, sdsname, file);
601  exit(1);
602  }
603  status = SDendaccess(sds_id);
604  sses->nlat = dims[0];
605 
606  /* terminate access to the SD interface and close the file */
607  status = SDend(sd_id);
608 }
609 
610 /* ------------------------------------------------------------------------------------------------ */
611 /* load_sses_sstv6v3() - loads the specified SSES (sensor-specific error stats) table for VIIRS sst3 */
612 
613 /* ------------------------------------------------------------------------------------------------ */
614 void load_sses_sstv6v3(int32_t sensorID, char *file, ssestabstrv6v3 *sses) {
615  int32 sd_id;
616  int32 sds_id;
617  int32 status;
618  int32 rank;
619  int32 nt;
620  int32 nattrs;
621  int32 dims[NSSESDIMv6v];
622  int32 start[NSSESDIMv6v];
623 
624  char name[H4_MAX_NC_NAME] = "";
625  char sdsname[H4_MAX_NC_NAME] = "";
626 
627  memset(start, 0, NSSESDIMv6v * sizeof (int32));
628 
629  if (strcmp(file, "") == 0) {
630  printf("\nNo SSES data provided for this sensor.\n");
631  haveSSES = 0;
632  return;
633  }
634 
635  printf("\nLoading SSES table from %s\n", file);
636 
637  /* open the file and initiate the SD interface */
638  sd_id = SDstart(file, DFACC_RDONLY);
639  if (sd_id == -1) {
640  printf("-E- %s line %d: Error opening file %s.\n", __FILE__, __LINE__,
641  file);
642  exit(1);
643  }
644 
645  /* read the bias and standard deviation */
646 
647  strcpy(sdsname, "bias");
648  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
649  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
650  if (rank != NSSESDIMv6v) {
651  printf(
652  "-E- %s line %d: Table dimensions for %s do not match expectation. got: %d expected %d\n",
653  __FILE__, __LINE__, sdsname, rank, NSSESDIMv6v);
654  exit(1);
655  }
656  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
657  || dims[2] != NDIFFMAXv6v3 || dims[3] != NSENZMAXv6v
658  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v3
659  || dims[6] != NSSTMAXv6v) {
660  printf(
661  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
662  __FILE__, __LINE__, sdsname);
663  exit(1);
664  }
665  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
666  if (status != 0) {
667  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
668  __LINE__, sdsname, file);
669  exit(1);
670  }
671  status = SDendaccess(sds_id);
672 
673  strcpy(sdsname, "stdv");
674  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
675  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
676  if (rank != NSSESDIMv6v) {
677  printf(
678  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
679  __FILE__, __LINE__, sdsname);
680  exit(1);
681  }
682  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
683  || dims[2] != NDIFFMAXv6v3 || dims[3] != NSENZMAXv6v
684  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v3
685  || dims[6] != NSSTMAXv6v) {
686  printf(
687  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
688  __FILE__, __LINE__, sdsname);
689  exit(1);
690  }
691  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
692  if (status != 0) {
693  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
694  __LINE__, sdsname, file);
695  exit(1);
696  }
697  status = SDendaccess(sds_id);
698 
699  /* new hypercubes have bias_mean and counts also */
700 
701  strcpy(sdsname, "bias_mean");
702  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
703  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
704  if (rank != NSSESDIMv6v) {
705  printf(
706  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
707  __FILE__, __LINE__, sdsname);
708  exit(1);
709  }
710  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
711  || dims[2] != NDIFFMAXv6v3 || dims[3] != NSENZMAXv6v
712  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v3
713  || dims[6] != NSSTMAXv6v) {
714  printf(
715  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
716  __FILE__, __LINE__, sdsname);
717  exit(1);
718  }
719  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias_mean);
720  if (status != 0) {
721  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
722  __LINE__, sdsname, file);
723  exit(1);
724  }
725  status = SDendaccess(sds_id);
726 
727  strcpy(sdsname, "counts");
728  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
729  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
730  if (rank != NSSESDIMv6v) {
731  printf(
732  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
733  __FILE__, __LINE__, sdsname);
734  exit(1);
735  }
736  if (dims[0] != NQUALMAXv6v || dims[1] != NLATMAXv6v
737  || dims[2] != NDIFFMAXv6v3 || dims[3] != NSENZMAXv6v
738  || dims[4] != NQUARMAXv6v || dims[5] != NDAYMAXv6v3
739  || dims[6] != NSSTMAXv6v) {
740  printf(
741  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
742  __FILE__, __LINE__, sdsname);
743  exit(1);
744  }
745  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->counts);
746  if (status != 0) {
747  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
748  __LINE__, sdsname, file);
749  exit(1);
750  }
751  status = SDendaccess(sds_id);
752 
753  /* save the table dimensions */
754  sses->nqual = dims[0];
755  sses->nlat = dims[1];
756  sses->ndiff = dims[2];
757  sses->nsenz = dims[3];
758  sses->nquar = dims[4];
759  sses->nday = dims[5];
760  sses->nsst = dims[6];
761 
762  /* read the table indice ranges */
763 
764  strcpy(sdsname, "sst");
765  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
766  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
767  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
768  if (status != 0) {
769  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
770  __LINE__, sdsname, file);
771  exit(1);
772  }
773  status = SDendaccess(sds_id);
774  sses->nsst = dims[0];
775 
776  strcpy(sdsname, "senz");
777  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
778  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
779  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
780  if (status != 0) {
781  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
782  __LINE__, sdsname, file);
783  exit(1);
784  }
785  status = SDendaccess(sds_id);
786  sses->nsenz = dims[0];
787 
788  strcpy(sdsname, "BTdiff");
789  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
790  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
791  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
792  if (status != 0) {
793  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
794  __LINE__, sdsname, file);
795  exit(1);
796  }
797  status = SDendaccess(sds_id);
798  sses->ndiff = dims[0];
799 
800  strcpy(sdsname, "lat");
801  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
802  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
803  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
804  if (status != 0) {
805  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
806  __LINE__, sdsname, file);
807  exit(1);
808  }
809  status = SDendaccess(sds_id);
810  sses->nlat = dims[0];
811 
812  /* terminate access to the SD interface and close the file */
813  status = SDend(sd_id);
814 }
815 
816 /* ---------------------------------------------------------------------------------------- */
817 /* load_sses_sstv6m() - loads the specified SSES (sensor-specific error stats) table */
818 
819 /* ---------------------------------------------------------------------------------------- */
820 void load_sses_sstv6m(int32_t sensorID, char *file, ssestabstrv6m *sses) {
821  int32 sd_id;
822  int32 sds_id;
823  int32 status;
824  int32 rank;
825  int32 nt;
826  int32 nattrs;
827  int32 dims[NSSESDIMv6m];
828  int32 start[NSSESDIMv6m];
829 
830  char name[H4_MAX_NC_NAME] = "";
831  char sdsname[H4_MAX_NC_NAME] = "";
832 
833  memset(start, 0, NSSESDIMv6m * sizeof (int32));
834 
835  if (strcmp(file, "") == 0) {
836  printf("\nNo SSES data provided for this sensor.\n");
837  haveSSES = 0;
838  return;
839  }
840 
841  printf("\nLoading SSES table from %s\n", file);
842 
843  /* open the file and initiate the SD interface */
844  sd_id = SDstart(file, DFACC_RDONLY);
845  if (sd_id == -1) {
846  printf("-E- %s line %d: Error opening file %s.\n", __FILE__, __LINE__,
847  file);
848  exit(1);
849  }
850 
851  /* read the bias and standard deviation */
852 
853  strcpy(sdsname, "bias");
854  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
855  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
856  if (rank != NSSESDIMv6m) {
857  printf(
858  "-E- %s line %d: Table dimensions for %s do not match expectation. got: %d expected %d\n",
859  __FILE__, __LINE__, sdsname, rank, NSSESDIMv6m);
860  exit(1);
861  }
862  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
863  || dims[2] != NDIFFMAXv6m || dims[3] != NSENZMAXv6m
864  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m
865  || dims[6] != NSSTMAXv6m) {
866  printf("%d,%d,%d,%d,%d,%d,%d\n",dims[0],dims[1],dims[2],dims[3],dims[4],dims[5],dims[6]);
867  printf(
868  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
869  __FILE__, __LINE__, sdsname);
870  exit(1);
871  }
872  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
873  if (status != 0) {
874  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
875  __LINE__, sdsname, file);
876  exit(1);
877  }
878  status = SDendaccess(sds_id);
879 
880  strcpy(sdsname, "stdv");
881  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
882  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
883  if (rank != NSSESDIMv6m) {
884  printf(
885  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
886  __FILE__, __LINE__, sdsname);
887  exit(1);
888  }
889  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
890  || dims[2] != NDIFFMAXv6m || dims[3] != NSENZMAXv6m
891  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m
892  || dims[6] != NSSTMAXv6m) {
893  printf(
894  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
895  __FILE__, __LINE__, sdsname);
896  exit(1);
897  }
898  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
899  if (status != 0) {
900  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
901  __LINE__, sdsname, file);
902  exit(1);
903  }
904  status = SDendaccess(sds_id);
905 
906  /* new hypercubes have bias_mean and counts also */
907 
908  strcpy(sdsname, "bias_mean");
909  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
910  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
911  if (rank != NSSESDIMv6m) {
912  printf(
913  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
914  __FILE__, __LINE__, sdsname);
915  exit(1);
916  }
917  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
918  || dims[2] != NDIFFMAXv6m || dims[3] != NSENZMAXv6m
919  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m
920  || dims[6] != NSSTMAXv6m) {
921  printf(
922  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
923  __FILE__, __LINE__, sdsname);
924  exit(1);
925  }
926 // status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias_mean);
927 // if (status != 0) {
928 // printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
929 // __LINE__, sdsname, file);
930 // exit(1);
931 // }
932 // status = SDendaccess(sds_id);
933 //
934 // strcpy(sdsname, "counts");
935 // sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
936 // status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
937 // if (rank != NSSESDIMv6m) {
938 // printf(
939 // "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
940 // __FILE__, __LINE__, sdsname);
941 // exit(1);
942 // }
943 // if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
944 // || dims[2] != NDIFFMAXv6m || dims[3] != NSENZMAXv6m
945 // || dims[4] != NQUARMAXv6m || dims[5] != NDAYMAXv6m
946 // || dims[6] != NSSTMAXv6m) {
947 // printf(
948 // "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
949 // __FILE__, __LINE__, sdsname);
950 // exit(1);
951 // }
952 // status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->counts);
953 // if (status != 0) {
954 // printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
955 // __LINE__, sdsname, file);
956 // exit(1);
957 // }
958  status = SDendaccess(sds_id);
959 
960  /* save the table dimensions */
961  sses->nqual = dims[0];
962  sses->nlat = dims[1];
963  sses->ndiff = dims[2];
964  sses->nsenz = dims[3];
965  sses->nquar = dims[4];
966  sses->nday = dims[5];
967  sses->nsst = dims[6];
968 
969  /* read the table indice ranges */
970 
971  strcpy(sdsname, "sst");
972  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
973  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
974  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
975  if (status != 0) {
976  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
977  __LINE__, sdsname, file);
978  exit(1);
979  }
980  status = SDendaccess(sds_id);
981  sses->nsst = dims[0];
982 
983  strcpy(sdsname, "senz");
984  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
985  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
986  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
987  if (status != 0) {
988  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
989  __LINE__, sdsname, file);
990  exit(1);
991  }
992  status = SDendaccess(sds_id);
993  sses->nsenz = dims[0];
994 
995  strcpy(sdsname, "BTdiff");
996  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
997  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
998  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
999  if (status != 0) {
1000  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1001  __LINE__, sdsname, file);
1002  exit(1);
1003  }
1004  status = SDendaccess(sds_id);
1005  sses->ndiff = dims[0];
1006 
1007  strcpy(sdsname, "lat");
1008  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1009  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1010  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
1011  if (status != 0) {
1012  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1013  __LINE__, sdsname, file);
1014  exit(1);
1015  }
1016  status = SDendaccess(sds_id);
1017  sses->nlat = dims[0];
1018 
1019  strcpy(sdsname, "lat");
1020  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1021  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1022  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
1023  if (status != 0) {
1024  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1025  __LINE__, sdsname, file);
1026  exit(1);
1027  }
1028  status = SDendaccess(sds_id);
1029  sses->nlat = dims[0];
1030 
1031  /* terminate access to the SD interface and close the file */
1032  status = SDend(sd_id);
1033 }
1034 
1035 /* ---------------------------------------------------------------------------------------- */
1036 /* load_sses_sstv6m4() - loads the specified SST4 SSES (sensor-specific error stats) table */
1037 
1038 /* ---------------------------------------------------------------------------------------- */
1039 void load_sses_sstv6m4(int32_t sensorID, char *file, ssestabstrv6m4 *sses) {
1040  int32 sd_id;
1041  int32 sds_id;
1042  int32 status;
1043  int32 rank;
1044  int32 nt;
1045  int32 nattrs;
1046  int32 dims[NSSESDIMv6m];
1047  int32 start[NSSESDIMv6m];
1048 
1049  char name[H4_MAX_NC_NAME] = "";
1050  char sdsname[H4_MAX_NC_NAME] = "";
1051 
1052  memset(start, 0, NSSESDIMv6m * sizeof (int32));
1053 
1054  if (strcmp(file, "") == 0) {
1055  printf("\nNo SSES data provided for this sensor.\n");
1056  haveSSES = 0;
1057  return;
1058  }
1059 
1060  printf("\nLoading SSES table from %s\n", file);
1061 
1062  /* open the file and initiate the SD interface */
1063  sd_id = SDstart(file, DFACC_RDONLY);
1064  if (sd_id == -1) {
1065  printf("-E- %s line %d: Error opening file %s.\n", __FILE__, __LINE__,
1066  file);
1067  exit(1);
1068  }
1069 
1070  /* read the bias and standard deviation */
1071 
1072  strcpy(sdsname, "bias");
1073  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1074  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1075  if (rank != NSSESDIMv6m) {
1076  printf(
1077  "-E- %s line %d: Table dimensions for %s do not match expectation. got: %d expected %d\n",
1078  __FILE__, __LINE__, sdsname, rank, NSSESDIMv6m);
1079  exit(1);
1080  }
1081  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
1082  || dims[2] != NDIFFMAXv6m4 || dims[3] != NSENZMAXv6m
1083  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m4
1084  || dims[6] != NSSTMAXv6m) {
1085  printf("%d,%d,%d,%d,%d,%d,%d\n",dims[0],dims[1],dims[2],dims[3],dims[4],dims[5],dims[6]);
1086  printf(
1087  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1088  __FILE__, __LINE__, sdsname);
1089  exit(1);
1090  }
1091  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
1092  if (status != 0) {
1093  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1094  __LINE__, sdsname, file);
1095  exit(1);
1096  }
1097  status = SDendaccess(sds_id);
1098 
1099  strcpy(sdsname, "stdv");
1100  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1101  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1102  if (rank != NSSESDIMv6m) {
1103  printf(
1104  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1105  __FILE__, __LINE__, sdsname);
1106  exit(1);
1107  }
1108  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
1109  || dims[2] != NDIFFMAXv6m4 || dims[3] != NSENZMAXv6m
1110  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m4
1111  || dims[6] != NSSTMAXv6m) {
1112  printf(
1113  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1114  __FILE__, __LINE__, sdsname);
1115  exit(1);
1116  }
1117  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
1118  if (status != 0) {
1119  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1120  __LINE__, sdsname, file);
1121  exit(1);
1122  }
1123  status = SDendaccess(sds_id);
1124 
1125  /* new hypercubes have bias_mean and counts also */
1126 
1127  strcpy(sdsname, "bias_mean");
1128  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1129  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1130  if (rank != NSSESDIMv6m) {
1131  printf(
1132  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1133  __FILE__, __LINE__, sdsname);
1134  exit(1);
1135  }
1136  if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
1137  || dims[2] != NDIFFMAXv6m || dims[3] != NSENZMAXv6m
1138  || dims[4] != NQUARMAXv6m || dims[5] > NDAYMAXv6m
1139  || dims[6] != NSSTMAXv6m) {
1140  printf(
1141  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1142  __FILE__, __LINE__, sdsname);
1143  exit(1);
1144  }
1145 // status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias_mean);
1146 // if (status != 0) {
1147 // printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1148 // __LINE__, sdsname, file);
1149 // exit(1);
1150 // }
1151 // status = SDendaccess(sds_id);
1152 //
1153 // strcpy(sdsname, "counts");
1154 // sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1155 // status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1156 // if (rank != NSSESDIMv6m) {
1157 // printf(
1158 // "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1159 // __FILE__, __LINE__, sdsname);
1160 // exit(1);
1161 // }
1162 // if (dims[0] != NQUALMAXv6m || dims[1] != NLATMAXv6m
1163 // || dims[2] != NDIFFMAXv6m4 || dims[3] != NSENZMAXv6m
1164 // || dims[4] != NQUARMAXv6m || dims[5] != NDAYMAXv6m4
1165 // || dims[6] != NSSTMAXv6m) {
1166 // printf(
1167 // "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1168 // __FILE__, __LINE__, sdsname);
1169 // exit(1);
1170 // }
1171 // status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->counts);
1172 // if (status != 0) {
1173 // printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1174 // __LINE__, sdsname, file);
1175 // exit(1);
1176 // }
1177  status = SDendaccess(sds_id);
1178 
1179  /* save the table dimensions */
1180  sses->nqual = dims[0];
1181  sses->nlat = dims[1];
1182  sses->ndiff = dims[2];
1183  sses->nsenz = dims[3];
1184  sses->nquar = dims[4];
1185  sses->nday = dims[5];
1186  sses->nsst = dims[6];
1187 
1188  /* read the table indice ranges */
1189 
1190  strcpy(sdsname, "sst");
1191  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1192  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1193  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
1194  if (status != 0) {
1195  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1196  __LINE__, sdsname, file);
1197  exit(1);
1198  }
1199  status = SDendaccess(sds_id);
1200  sses->nsst = dims[0];
1201 
1202  strcpy(sdsname, "senz");
1203  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1204  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1205  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
1206  if (status != 0) {
1207  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1208  __LINE__, sdsname, file);
1209  exit(1);
1210  }
1211  status = SDendaccess(sds_id);
1212  sses->nsenz = dims[0];
1213 
1214  strcpy(sdsname, "BTdiff");
1215  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1216  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1217  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
1218  if (status != 0) {
1219  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1220  __LINE__, sdsname, file);
1221  exit(1);
1222  }
1223  status = SDendaccess(sds_id);
1224  sses->ndiff = dims[0];
1225 
1226  strcpy(sdsname, "lat");
1227  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1228  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1229  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
1230  if (status != 0) {
1231  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1232  __LINE__, sdsname, file);
1233  exit(1);
1234  }
1235  status = SDendaccess(sds_id);
1236  sses->nlat = dims[0];
1237 
1238  strcpy(sdsname, "lat");
1239  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1240  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1241  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
1242  if (status != 0) {
1243  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1244  __LINE__, sdsname, file);
1245  exit(1);
1246  }
1247  status = SDendaccess(sds_id);
1248  sses->nlat = dims[0];
1249 
1250  /* terminate access to the SD interface and close the file */
1251  status = SDend(sd_id);
1252 }
1253 
1254 /* ---------------------------------------------------------------------------------------- */
1255 /* load_sses_sstv6a() - loads the specified AVHRR SSES (sensor-specific error stats) table */
1256 
1257 /* ---------------------------------------------------------------------------------------- */
1258 void load_sses_sstv6a(int32_t sensorID, char *file, ssestabstrv6a *sses) {
1259  int32 sd_id;
1260  int32 sds_id;
1261  int32 status;
1262  int32 rank;
1263  int32 nt;
1264  int32 nattrs;
1265  int32 dims[NSSESDIMv6a];
1266  int32 start[NSSESDIMv6a];
1267 
1268  char name[H4_MAX_NC_NAME] = "";
1269  char sdsname[H4_MAX_NC_NAME] = "";
1270 
1271  memset(start, 0, NSSESDIMv6a * sizeof (int32));
1272 
1273  if (strcmp(file, "") == 0) {
1274  printf("\nNo SSES data provided for this sensor.\n");
1275  haveSSES = 0;
1276  return;
1277  }
1278 
1279  printf("\nLoading SSES table from %s\n", file);
1280 
1281  /* open the file and initiate the SD interface */
1282  sd_id = SDstart(file, DFACC_RDONLY);
1283  if (sd_id == -1) {
1284  printf("-E- %s line %d: Error opening file %s.\n", __FILE__, __LINE__,
1285  file);
1286  exit(1);
1287  }
1288 
1289  /* read the bias and standard deviation */
1290 
1291  strcpy(sdsname, "bias");
1292  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1293  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1294  if (rank != NSSESDIMv6a) {
1295  printf(
1296  "-E- %s line %d: Table dimensions for %s do not match expectation. got: %d expected %d\n",
1297  __FILE__, __LINE__, sdsname, rank, NSSESDIMv6a);
1298  exit(1);
1299  }
1300  if (dims[0] != NQUALMAXv6a || dims[1] != NLATMAXv6a
1301  || dims[2] != NDIFFMAXv6a || dims[3] != NSENZMAXv6a
1302  || dims[4] != NQUARMAXv6a || dims[5] != NSSTMAXv6a) {
1303  printf(
1304  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1305  __FILE__, __LINE__, sdsname);
1306  exit(1);
1307  }
1308  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
1309  if (status != 0) {
1310  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1311  __LINE__, sdsname, file);
1312  exit(1);
1313  }
1314  status = SDendaccess(sds_id);
1315 
1316  strcpy(sdsname, "stdv");
1317  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1318  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1319  if (rank != NSSESDIMv6a) {
1320  printf(
1321  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1322  __FILE__, __LINE__, sdsname);
1323  exit(1);
1324  }
1325  if (dims[0] != NQUALMAXv6a || dims[1] != NLATMAXv6a
1326  || dims[2] != NDIFFMAXv6a || dims[3] != NSENZMAXv6a
1327  || dims[4] != NQUARMAXv6a || dims[5] != NSSTMAXv6a) {
1328  printf(
1329  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1330  __FILE__, __LINE__, sdsname);
1331  exit(1);
1332  }
1333  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
1334  if (status != 0) {
1335  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1336  __LINE__, sdsname, file);
1337  exit(1);
1338  }
1339  status = SDendaccess(sds_id);
1340 
1341  /* new hypercubes have bias_mean and counts also */
1342 
1343  strcpy(sdsname, "bias_mean");
1344  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1345  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1346  if (rank != NSSESDIMv6a) {
1347  printf(
1348  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1349  __FILE__, __LINE__, sdsname);
1350  exit(1);
1351  }
1352  if (dims[0] != NQUALMAXv6a || dims[1] != NLATMAXv6a
1353  || dims[2] != NDIFFMAXv6a || dims[3] != NSENZMAXv6a
1354  || dims[4] != NQUARMAXv6a || dims[5] != NSSTMAXv6a) {
1355  printf(
1356  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1357  __FILE__, __LINE__, sdsname);
1358  exit(1);
1359  }
1360  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias_mean);
1361  if (status != 0) {
1362  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1363  __LINE__, sdsname, file);
1364  exit(1);
1365  }
1366  status = SDendaccess(sds_id);
1367 
1368  strcpy(sdsname, "counts");
1369  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1370  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1371  if (rank != NSSESDIMv6a) {
1372  printf(
1373  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1374  __FILE__, __LINE__, sdsname);
1375  exit(1);
1376  }
1377  if (dims[0] != NQUALMAXv6a || dims[1] != NLATMAXv6a
1378  || dims[2] != NDIFFMAXv6a || dims[3] != NSENZMAXv6a
1379  || dims[4] != NQUARMAXv6a || dims[5] != NSSTMAXv6a) {
1380  printf(
1381  "-E- %s line %d: Table dimensions for %s do not match expectation.\n",
1382  __FILE__, __LINE__, sdsname);
1383  exit(1);
1384  }
1385  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->counts);
1386  if (status != 0) {
1387  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1388  __LINE__, sdsname, file);
1389  exit(1);
1390  }
1391  status = SDendaccess(sds_id);
1392 
1393  /* save the table dimensions */
1394  sses->nqual = dims[0];
1395  sses->nlat = dims[1];
1396  sses->ndiff = dims[2];
1397  sses->nsenz = dims[3];
1398  sses->nquar = dims[4];
1399  sses->nsst = dims[5];
1400 
1401  /* read the table indice ranges */
1402 
1403  strcpy(sdsname, "sst");
1404  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1405  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1406  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
1407  if (status != 0) {
1408  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1409  __LINE__, sdsname, file);
1410  exit(1);
1411  }
1412  status = SDendaccess(sds_id);
1413  sses->nsst = dims[0];
1414 
1415  strcpy(sdsname, "senz");
1416  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1417  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1418  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
1419  if (status != 0) {
1420  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1421  __LINE__, sdsname, file);
1422  exit(1);
1423  }
1424  status = SDendaccess(sds_id);
1425  sses->nsenz = dims[0];
1426 
1427  strcpy(sdsname, "BTdiff");
1428  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1429  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1430  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
1431  if (status != 0) {
1432  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1433  __LINE__, sdsname, file);
1434  exit(1);
1435  }
1436  status = SDendaccess(sds_id);
1437  sses->ndiff = dims[0];
1438 
1439  strcpy(sdsname, "lat");
1440  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
1441  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
1442  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
1443  if (status != 0) {
1444  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
1445  __LINE__, sdsname, file);
1446  exit(1);
1447  }
1448  status = SDendaccess(sds_id);
1449  sses->nlat = dims[0];
1450 
1451  /* terminate access to the SD interface and close the file */
1452  status = SDend(sd_id);
1453 }
1454 
1455 void isV5coef(l2str *l2rec) {
1456  FILE *fp;
1457  char line [200] = "";
1458  char mission[5];
1459  int datechk;
1460  int32_t ncid;
1461 
1462  /* Open the coefficient file */
1463  if (nc_open(input->sstcoeffile, NC_NOWRITE, &ncid) != NC_NOERR) {
1464 
1465  if ((fp = fopen(input->sstcoeffile, "r")) == NULL) {
1466  fprintf(stderr,
1467  "-E- %s line %d: unable to open sst coef file %s for reading\n",
1468  __FILE__, __LINE__, input->sstcoeffile);
1469  exit(1);
1470  }
1471 
1472  while (fgets(line, 200, fp)) {
1473  if (line[0] == '#') {
1474  continue;
1475  } else {
1476  sscanf(line, "%4s %d", mission, &datechk);
1477  break;
1478  }
1479  }
1480  fclose(fp);
1481  /*
1482  * If it's not a VIIRS file, and the datechk is true, it's a V5 table
1483  * VIIRS still uses dates, but we don't have v5 coef for VIIRS, so for now,
1484  * this simple check should work -
1485  */
1486  if ((strcmp(mission, "VIIR") != 0) && (datechk > 1000)) {
1487  isV5 = 1;
1488  if (want_verbose)
1489  printf("Using V5 coefficients to compute SST\n");
1490  }
1491  } else {
1492  nc_close(ncid);
1493  }
1494 }
1495 
1496 /* ----------------------------------------------------------------------------------- */
1497 /* init_sst() - Initialize for SST processing */
1498 /* B. Franz, SAIC, August 2005 */
1499 
1500 /* ----------------------------------------------------------------------------------- */
1501 void init_sst(l2str *l2rec) {
1502  extern l1qstr l1que;
1503  int32_t npix = l2rec->l1rec->npix;
1504 
1505  evalmask = input->evalmask;
1506 
1507  nbvis = l2rec->l1rec->l1file->nbands;
1508  nbir = NBANDSIR;
1509  ibred = bindex_get(678); /* 9 for modisa */
1510  ib07 = bindex_get(748); /* 10 for modisa */
1511  ib08 = bindex_get(850);
1512  if (ib08 < 0) {
1513  ib08 = bindex_get(865);
1514  }
1515  ib16 = bindex_get(1640); /* 14 for modisa */
1516  ib37 = bindex_get(3750) - l2rec->l1rec->l1file->nbands; /* 0 for modisa */
1517  ib39 = bindex_get(3959) - l2rec->l1rec->l1file->nbands; /* 1 for modisa */
1518  ib40 = bindex_get(4050) - l2rec->l1rec->l1file->nbands; /* 2 for modisa */
1519  ib67 = bindex_get(6715) - l2rec->l1rec->l1file->nbands;
1520  ib73 = bindex_get(7325) - l2rec->l1rec->l1file->nbands;
1521  ib85 = bindex_get(8550) - l2rec->l1rec->l1file->nbands; /* 5 for modisa */
1522  ib11 = bindex_get(11000) - l2rec->l1rec->l1file->nbands; /* 6 for modisa */
1523  ib12 = bindex_get(12000) - l2rec->l1rec->l1file->nbands; /* 7 for modisa */
1524 
1525  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
1526  if (strncmp(l2rec->l1rec->l1file->spatialResolution, "4.6km", 5) == 0) {
1527  fullscanpix = 409;
1528  } else {
1529  fullscanpix = 2048;
1530  }
1531  /* avhrr wavelenths are 630, 855, 3700, 11000, 12000 */
1532  /* or, are they different for different satellites? */
1533  /* ibred should be index 0 for channel 1 data */
1534  ibred = bindex_get(630);
1535  //Hide ib39, and ib40 as 3700 is put in here from above...
1536  ib39 = -1;
1537  ib40 = -1;
1538  /* channel 3a, 1610, exists for daytime NO16, NO17, NO18, NO19 */
1539  /* all other cases it's 3700 so msl12_sensor.dat has 3700 for channel 3 */
1540  // ib16 = ib08 + 1;
1541  /* some limits are different for AVHRR */
1542  solznight = SOLZNIGHTA;
1543  hisenz = hisenza;
1544  vhisenz = vhisenza;
1545  Btmin = Btmina;
1546  Btmax = Btmaxa;
1547  SSTmax = SSTmaxa;
1548  SSTdiff = SSTdiffa;
1549  switch (l2rec->l1rec->l1file->subsensorID) {
1550  case NO07:
1551  case NO11:
1552  case NO14:
1553  case NO15:
1554  Bt11unif1 = 0.6;
1555  break;
1556  case NO09:
1557  Bt11unif1 = 0.7;
1558  break;
1559  default:
1560  Bt11unif1 = 1.0;
1561  break;
1562  }
1563  }
1564 
1565  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
1566  fullscanpix = 3200;
1567  cldbox = cldboxv;
1568  btbox = btboxv;
1569  // Red and NIR bands in VIIRS is far enough from MODIS to require a separate call
1570  ibred = bindex_get(672); /* SVM05 */
1571  ib16 = bindex_get(1601); /* SVM10 */
1572 
1573  //Hide ib39, as 3700 is put in here from above...
1574  ib39 = -1;
1575 
1576  cldthresh = cldthreshv;
1577  /* viirs Bt's are in K, not Deg C (not so using OBPG reader - we read the radiances and create the BTs)*/
1578  // Btmin = Btminv;
1579  // Btmax = Btmaxv;
1580  // Btmax40 = Btmax40v;
1581  }
1582 
1583  if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
1584  if (input->resolution == 250) {
1585  fullscanpix = 5416;
1586  } else if (input->resolution == 500) {
1587  fullscanpix = 2708;
1588  }
1589  }
1590 
1591  if (ib11 < 0 || ib12 < 0)
1592  haveSST = 0;
1593  else
1594  haveSST = 1;
1595 
1596  if (ib39 < 0 || ib40 < 0)
1597  haveSST4 = 0;
1598  else
1599  haveSST4 = 1;
1600 
1601  if (!haveSST && !haveSST4) {
1602  fprintf(stderr, "-E- %s line %d: no SST bands found.\n",
1603  __FILE__, __LINE__);
1604  exit(1);
1605  }
1606 
1607  /* this is the value that a saturated radiance would be assigned */
1608  if (ibred < 0)
1609  haveRed = 0;
1610  else {
1611  haveRed = 1;
1612  // satred = 1000.0/input->gain[ibred]-1.0;
1613  satred = 1000.0 - 1.0;
1614  }
1615  /* make sure l1que.nq is not smaller than btbox, can only calc sst's if have Bt's */
1616  if (l1que.nq < btbox) {
1617  fprintf(stderr, "-E- %s line %d: filter queue (l1que) is too small. Have %d need %d.\n",
1618  __FILE__, __LINE__, l1que.nq, btbox);
1619  exit(1);
1620  }
1621 
1622  /* allocate private arrays for a single scan line */
1623  if (haveSST) {
1624  sstq = (s_array*) calloc(l1que.nq, sizeof (s_array));
1625  flags_sst = (int16*) calloc(npix, sizeof (int16));
1626  qual_sst = (int8*) calloc(npix, sizeof (int8));
1627  bias_sst = (float*) calloc(npix, sizeof (float));
1628  stdv_sst = (float*) calloc(npix, sizeof (float));
1629  bias_mean_sst = (float*) calloc(npix, sizeof (float));
1630  bias_counts_sst = (int16*) calloc(npix, sizeof (int16));
1631  dsdi_correction = (float*) calloc(npix, sizeof (float));
1632 
1633  treesum = (float*) calloc(npix, sizeof (float));
1634 
1635  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
1636  load_sses_sstv6a(l2rec->l1rec->l1file->sensorID, input->sstssesfile,
1637  &sses_sstv6a);
1638  } else {
1639  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
1640  load_sses_sstv6v(l2rec->l1rec->l1file->sensorID, input->sstssesfile,
1641  &sses_sstv6v);
1642  } else {
1643  /* modis */
1644  load_sses_sstv6m(l2rec->l1rec->l1file->sensorID, input->sstssesfile,
1645  &sses_sstv6m);
1646  }
1647  }
1648 
1649  // VIIRSN can do a triple window sst algorithm...
1650  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
1651  sst3q = (s_array*) calloc(l1que.nq, sizeof (s_array));
1652  flags_sst3 = (int16*) calloc(npix, sizeof (int16));
1653  qual_sst3 = (int8*) calloc(npix, sizeof (int8));
1654  bias_sst3 = (float*) calloc(npix, sizeof (float));
1655  stdv_sst3 = (float*) calloc(npix, sizeof (float));
1656  bias_mean_sst3 = (float*) calloc(npix, sizeof (float));
1657  bias_counts_sst3 = (int16*) calloc(npix, sizeof (int16));
1658  load_sses_sstv6v3(l2rec->l1rec->l1file->sensorID, input->sst3ssesfile, &sses_sst3v6v3);
1659 
1660  }
1661 
1662  }
1663 
1664  if (haveSST4) {
1665  sst4q = (s_array*) calloc(l1que.nq, sizeof (s_array));
1666  flags_sst4 = (int16*) calloc(npix, sizeof (int16));
1667  qual_sst4 = (int8*) calloc(npix, sizeof (int8));
1668  bias_sst4 = (float*) calloc(npix, sizeof (float));
1669  stdv_sst4 = (float*) calloc(npix, sizeof (float));
1670  bias_mean_sst4 = (float*) calloc(npix, sizeof (float));
1671  bias_counts_sst4 = (int16*) calloc(npix, sizeof (int16));
1672  load_sses_sstv6m4(l2rec->l1rec->l1file->sensorID, input->sst4ssesfile,
1673  &sses_sst4v6m);
1674 
1675  d3940ref = (float*) calloc(npix, sizeof (float));
1676  }
1677  rhoCirrus_maxmin = (float*) calloc(npix, sizeof (float));
1678  rhoCirrus_min = (float*) calloc(npix, sizeof (float));
1679  rhoCirrus_max = (float*) calloc(npix, sizeof (float));
1680 
1681  if (haveRed) {
1682  LtRED_maxmin = (float*) calloc(npix, sizeof (float));
1683  rhotRED_maxmin = (float*) calloc(npix, sizeof (float));
1684  rhotRED_min = (float*) calloc(npix, sizeof (float));
1685  rhotRED_max = (float*) calloc(npix, sizeof (float));
1686  rhotRED = (float*) calloc(npix, sizeof (float));
1687  }
1688  if (ib11 >= 0) {
1689  Bt11_maxmin = (float*) calloc(npix, sizeof (float));
1690  Bt11_max = (float*) calloc(npix, sizeof (float));
1691  Bt11_min = (float*) calloc(npix, sizeof (float));
1692  Bt11_stdev = (float*) calloc(npix, sizeof (float));
1693  }
1694  if (ib12 >= 0) {
1695  Bt12_maxmin = (float*) calloc(npix, sizeof (float));
1696  Bt12_min = (float*) calloc(npix, sizeof (float));
1697  }
1698  if (ib37 >= 0) {
1699  Bt37_maxmin = (float*) calloc(npix, sizeof (float));
1700  Bt37_stdev = (float*) calloc(npix, sizeof (float));
1701  }
1702  if (ib73 >= 0)
1703  Bt73_max = (float*) calloc(npix, sizeof (float));
1704 
1705  if (ib85 >= 0)
1706  Bt85_min = (float*) calloc(npix, sizeof (float));
1707 
1708  if (ib07 >= 0) {
1709  rhotNIR7_min = (float*) calloc(npix, sizeof (float));
1710  rhotNIR7 = (float*) calloc(npix, sizeof (float));
1711  }
1712  if (ib16 >= 0) {
1713  rhot16_min = (float*) calloc(npix, sizeof (float));
1714  rhot16 = (float*) calloc(npix, sizeof (float));
1715  }
1716  if (ib39 >= 0)
1717  Bt39_maxmin = (float*) calloc(npix, sizeof (float));
1718 
1719  if (ib40 >= 0) {
1720  Bt40_maxmin = (float*) calloc(npix, sizeof (float));
1721  Bt40_stdev = (float*) calloc(npix, sizeof (float));
1722  // if (l2rec->l1rec->l1file->sensorID == MODISA)
1723  // Bt40_avg = (float*) calloc(npix, sizeof(float));
1724  }
1725 
1726  sst_stdev = (float*) calloc(npix, sizeof (float));
1727 
1728  csstbox = l1que.nq / 2; /* csstbox is the center of the box of sst's */
1729  isV5coef(l2rec);
1730 }
1731 
1732 /* ----------------------------------------------------------------------------------- */
1733 /* sst_ran() - Determine if we have already run sst for this scan line */
1734 
1735 /* ----------------------------------------------------------------------------------- */
1736 int sst_ran(int recnum) {
1737  if (recnum == recnumSST)
1738  return 1;
1739  else
1740  return 0;
1741 }
1742 
1743 /* ----------------------------------------------------------------------------------- */
1744 /* read_v5_sst_coeff() - reads sst coefficients for the specified date */
1745 /* */
1746 /* B. Franz, SAIC, 11 August 2005 */
1747 
1748 /* ----------------------------------------------------------------------------------- */
1749 void read_v5_sst_coeff(l2str *l2rec, float **bounds, float **coef, float *sstrefoffday, float *sstrefoffnight) {
1750  char mission[5] = "";
1751  char mission2[5] = "";
1752 
1753  FILE *fp;
1754  char line [200] = "";
1755  char odates [8] = ""; /* yyyyddd */
1756  char sdates [8] = "";
1757  char edates [8] = "";
1758  char name [80];
1759  char value [80];
1760  char *p;
1761  char *p1;
1762  char *p2;
1763 
1764  int32 found = 0;
1765  int32 indx = 0;
1766  int32 sensorID = l2rec->l1rec->l1file->sensorID;
1767  double pasutime = l2rec->l1rec->scantime;
1768  int16_t year, day;
1769  double sec;
1770  unix2yds(pasutime, &year, &day, &sec);
1771 
1772  switch (sensorID) {
1773  case MODISA:
1774  strcpy(mission, "AQUA");
1775  break;
1776  case MODIST:
1777  strcpy(mission, "TERR");
1778  break;
1779  case AVHRR:
1780  strcpy(mission, xsatid2name(l2rec->l1rec->l1file->subsensorID));
1781  break;
1782  default:
1783  strcpy(mission, "AQUA");
1784  break;
1785  }
1786 
1787  /* get v5 coeffs for AVHRR or MODIS */
1788  fp = fopen(input->sstcoeffile, "r");
1789  /* Form date string */
1790 
1791 uselastyear:
1792  sprintf(odates, "%4d%03d", year, day);
1793 
1794  /* Loop through to find bounding times, for 2 sets of coeffs */
1795 
1796  indx = 0;
1797  while (fgets(line, 200, fp)) {
1798  if (line[0] == '#') {
1799  /* look for lines with: # variable = value */
1800  if (!(p = strchr(line, '=')))
1801  continue;
1802  p1 = line + 1; /* look for white space after # and before variable name */
1803  while (isspace(*p1))
1804  p1++;
1805  p2 = p - 1; /* look for white space before = and after variable name */
1806  while (isspace(*p2))
1807  p2--;
1808  /* get variable name from input line */
1809  strncpy(name, p1, p2 - p1 + 1);
1810  name[p2 - p1 + 1] = '\0';
1811 
1812  /*
1813  * Parse parameter value string
1814  */
1815  /* start at character after = and ignore white space before value */
1816  p1 = p + 1;
1817  while (isspace(*p1))
1818  p1++;
1819  p2 = p1;
1820  /* look for white space to determine end of value */
1821  while (!isspace(*p2))
1822  p2++;
1823  /* get value from input line */
1824  strncpy(value, p1, p2 - p1);
1825  value[p2 - p1] = '\0';
1826 
1827  } else {
1828  /* read sst v5 coeffs for AVHRR */
1829  if (strncmp(line, mission, 4) == 0) {
1830  sscanf(line, "%4s %7s %7s %f %f %f %f",
1831  mission2, sdates, edates, &coef[indx][0], &coef[indx][1], &coef[indx][2], &coef[indx][3]);
1832  coef[indx][4] = 0.0;
1833  if (strcmp(odates, sdates) >= 0 && (strcmp(odates, edates) <= 0
1834  || strcmp(edates, "0000000") == 0
1835  || strcmp(edates, "") == 0)) {
1836  indx++;
1837  if (indx == 2) {
1838  found = 1;
1839  break;
1840  }
1841  }
1842  }
1843  }
1844  }
1845 
1846  if (found == 0 && year > 2004) {
1847  printf("Warning: No SST coefficients available for %s, reverting to previous year.\n", odates);
1848  year--;
1849  rewind(fp);
1850  goto uselastyear;
1851  }
1852 
1853 
1854  fclose(fp);
1855 
1856  if (found == 1) {
1857  printf("Loading SST coefficients from %s:\n", input->sstcoeffile);
1858  printf("%s %s %6.3f %6.3f %6.3f %6.3f %6.3f\n", sdates, edates, coef[0][0], coef[0][1], coef[0][2], coef[0][3], coef[0][4]);
1859  printf("%s %s %6.3f %6.3f %6.3f %6.3f %6.3f\n\n", sdates, edates, coef[1][0], coef[1][1], coef[1][2], coef[1][3], coef[1][4]);
1860 
1861  printf(" sst reference day offset = %f\n", *sstrefoffday);
1862  printf(" sst reference night offset = %f\n", *sstrefoffnight);
1863  } else {
1864  fprintf(stderr,
1865  "-E- %s line %d: unable to locate valid SST coefficients for %s in %s\n",
1866  __FILE__, __LINE__, odates, input->sstcoeffile);
1867  exit(1);
1868  }
1869 
1870  return;
1871 }
1872 
1873 /* ----------------------------------------------------------------------------------- */
1874 /* read_sst_coeff() - reads sst coefficients for the specified date */
1875 /* */
1876 /* B. Franz, SAIC, 11 August 2005 */
1877 
1878 /* ----------------------------------------------------------------------------------- */
1879 void read_sst_coeff(l2str *l2rec, float **bounds, float **coef,
1880  float *sstrefoffday, float *sstrefoffnight) {
1881  char mission[5] = "";
1882  char mission2[5] = "";
1883 
1884  FILE *fp;
1885  char line[200] = "";
1886  char odatel[14] = ""; /* yyyydddhhmmss */
1887  char sdatel[14] = ""; /* yyyydddhhmmss */
1888  char edatel[14] = ""; /* yyyydddhhmmss */
1889  char odates[8] = ""; /* yyyyddd */
1890  char sdates[8] = "";
1891  char edates[8] = "";
1892  char stime[7] = "";
1893  char etime[7] = "";
1894  char *coeflabel[] ={"day dry ", "day moist ", "night dry ", "night moist "};
1895  char dorn[2] = ""; /* day or night, we assume DDNN, so don't really need it */
1896  char *ztime;
1897  char name[80];
1898  char value[80];
1899  char *p;
1900  char *p1;
1901  char *p2;
1902  int32 gotsstrefoffday = 0;
1903  int32 gotsstrefoffnight = 0;
1904  int32 found = 0;
1905  int32 indx = 0;
1906  // int32 tmonth = 0;
1907  int32 month;
1908  int32 leap;
1909  int32 sensorID = l2rec->l1rec->l1file->sensorID;
1910  double pasutime = l2rec->l1rec->scantime;
1911  int16_t year, day;
1912  double sec;
1913  unix2yds(pasutime, &year, &day, &sec);
1914 
1915 // int32_t year = *l2rec->year;
1916 // int32_t day = *l2rec->day;
1917 // int32_t msec = *l2rec->msec;
1918  int32 tmp1;
1919  float tmp2, tmp3;
1920  //float64 pasutime;
1921 
1922  switch (sensorID) {
1923  case MODISA:
1924  strcpy(mission, "AQUA");
1925  break;
1926  case MODIST:
1927  strcpy(mission, "TERR");
1928  break;
1929  case AVHRR:
1930  strcpy(mission, xsatid2name(l2rec->l1rec->l1file->subsensorID));
1931  break;
1932  case VIIRSN:
1933  strcpy(mission, "VIIR");
1934  break;
1935  default:
1936  strcpy(mission, "AQUA");
1937  break;
1938  }
1939 
1940  /* Open the file */
1941  if ((fp = fopen(input->sstcoeffile, "r")) == NULL) {
1942  fprintf(stderr,
1943  "-E- %s line %d: unable to open sst coef file %s for reading\n",
1944  __FILE__, __LINE__, input->sstcoeffile);
1945  exit(1);
1946  }
1947  /* Find month */
1948  leap = (isleap(year) == TRUE ? 1 : 0);
1949  if (tmonth < 0) {
1950  for (tmonth = 11; tmonth >= 0; tmonth--) {
1951  /* day is one based, StartOfMonth is zero based */
1952  if (day > StartOfMonth[leap][tmonth]) {
1953  break;
1954  }
1955  }
1956  }
1957 
1958  if (sensorID == VIIRSN) {
1959  /* VIIRS latband coeffs are this format: */
1960  /* sat, start yyyyddd, start hhmmss, end yyyyddd, end hhmmss, min lat, max lat, c0, c1, c2, c3 */
1961  /* VIIR 2012001 000000 2012069 235900 -90 -40 -4.202194 1.0196764 0.002352195 1.206611 */
1962 
1963  /* Form date string */
1964 
1965  /* msec is the scan line start time */
1966  /* should calculate the actual pixel time (extrapolate from msec) for field 57? */
1967  /* modis, viirs, and avhhr: year, day are per scan line, not just time of first scan */
1968  // pasutime = yds2unix(year, day, ((double) (msec)) / 1000.0);
1969  ztime = ydhmsf(pasutime, 'G');
1970  /* ztime is yyyydddhhmmssfff */
1971  strncpy(odatel, ztime, 13);
1972  odatel[13] = '\0';
1973 
1974  /* Loop through to find bounding times */
1975 
1976  indx = 0;
1977  while (fgets(line, 200, fp)) {
1978  if (line[0] == '#') {
1979  /* look for lines with: # variable = value */
1980  if (!(p = strchr(line, '=')))
1981  continue;
1982  p1 = line + 1; /* look for white space after # and before variable name */
1983  while (isspace(*p1))
1984  p1++;
1985  p2 = p - 1; /* look for white space before = and after variable name */
1986  while (isspace(*p2))
1987  p2--;
1988  /* get variable name from input line */
1989  strncpy(name, p1, p2 - p1 + 1);
1990  name[p2 - p1 + 1] = '\0';
1991 
1992  /*
1993  * Parse parameter value string
1994  */
1995  /* start at character after = and ignore white space before value */
1996  p1 = p + 1;
1997  while (isspace(*p1))
1998  p1++;
1999  p2 = p1;
2000  /* look for white space to determine end of value */
2001  while (!isspace(*p2))
2002  p2++;
2003  /* get value from input line */
2004  strncpy(value, p1, p2 - p1);
2005  value[p2 - p1] = '\0';
2006  /*
2007  * Copy value to appropriate variable
2008  */
2009  if (strcmp(name, "sstref_day_offset") == 0) {
2010  *sstrefoffday = (float) atof(value);
2011  gotsstrefoffday = 1;
2012  }
2013  if (strcmp(name, "sstref_night_offset") == 0) {
2014  *sstrefoffnight = (float) atof(value);
2015  gotsstrefoffnight = 1;
2016  }
2017 
2018  } else {
2019  /* read viirs sst latband coeffs */
2020  if (strncmp(line, mission, 4) == 0) {
2021  if (input->viirsnosisaf == 1) {
2022  sscanf(line, "%4s %7s %7s %1s %f %f %f %f %f %f %f",
2023  mission2, sdates, edates, dorn,
2024  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2025  &coef[indx][3], &coef[indx][4], &coef[indx][5],
2026  &coef[indx][6]);
2027  } else if (input->viirsnv7 >= 0) {
2028  /* all latband versions except v6.4.1 */
2029  sscanf(line,
2030  "%4s %7s %6s %7s %6s %f %f %f %f %f %f %f %f %f",
2031  mission2, sdates, stime, edates, etime,
2032  &bounds[indx][0], &bounds[indx][1],
2033  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2034  &coef[indx][3], &coef[indx][4], &coef[indx][5],
2035  &coef[indx][6]);
2036  sprintf(sdatel, "%s%s", sdates, stime);
2037  sprintf(edatel, "%s%s", edates, etime);
2038  if (strcmp(odatel, sdatel) >= 0
2039  && (strcmp(odatel, edatel) <= 0
2040  || strcmp(edates, "0000000") == 0
2041  || strcmp(edates, "") == 0)) {
2042  indx++;
2043  if (indx == 7) {
2044  found = 1;
2045  break;
2046  }
2047  }
2048  } else {
2049  /* viirs nlsst v6.4.1 has extra satz terms, but not mirror side */
2050  /* and is by month, not start and end dates */
2051  sscanf(line, "%4s %d %f %f %f %f %f %f %f %f",
2052  mission2, &month, &bounds[indx][0], &bounds[indx][1],
2053  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2054  &coef[indx][3], &coef[indx][5], &coef[indx][6]);
2055  /* no mirror side for viirs */
2056  coef[indx][4] = 0.0;
2057  if (month == tmonth + 1) {
2058  indx++;
2059  if (indx == 7) {
2060  found = 1;
2061  break;
2062  }
2063  }
2064  }
2065 
2066  }
2067  }
2068  }
2069 
2070  } else {
2071  /* get latband coeffs for AVHRR or MODIS */
2072 
2073  /* Loop through to find 6 sets of coeffs for this month */
2074  indx = 0;
2075  while (fgets(line, 200, fp)) {
2076  if (line[0] == '#') {
2077  /* look for lines with: # variable = value */
2078  if (!(p = strchr(line, '=')))
2079  continue;
2080  p1 = line + 1; /* look for white space after # and before variable name */
2081  while (isspace(*p1))
2082  p1++;
2083  p2 = p - 1; /* look for white space before = and after variable name */
2084  while (isspace(*p2))
2085  p2--;
2086  /* get variable name from input line */
2087  strncpy(name, p1, p2 - p1 + 1);
2088  name[p2 - p1 + 1] = '\0';
2089 
2090  /*
2091  * Parse parameter value string
2092  */
2093  /* start at character after = and ignore white space before value */
2094  p1 = p + 1;
2095  while (isspace(*p1))
2096  p1++;
2097  p2 = p1;
2098  /* look for white space to determine end of value */
2099  while (!isspace(*p2))
2100  p2++;
2101  /* get value from input line */
2102  strncpy(value, p1, p2 - p1);
2103  value[p2 - p1] = '\0';
2104  /*
2105  * Copy value to appropriate variable
2106  */
2107  if (strcmp(name, "sstref_day_offset") == 0) {
2108  *sstrefoffday = (float) atof(value);
2109  gotsstrefoffday = 1;
2110  }
2111  if (strcmp(name, "sstref_night_offset") == 0) {
2112  *sstrefoffnight = (float) atof(value);
2113  gotsstrefoffnight = 1;
2114  }
2115 
2116  } else {
2117  /* read sst latband coeffs for AVHRR or MODIS */
2118  if (strncmp(line, mission, 4) == 0) {
2119  if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
2120  /* terra and aqua nlsst has extra mirror side and satz terms */
2121  sscanf(line,
2122  "%4s %d %f %f %f %f %f %f %f %f %f %d %f %f",
2123  mission2, &month, &bounds[indx][0],
2124  &bounds[indx][1], &coef[indx][0],
2125  &coef[indx][1], &coef[indx][2], &coef[indx][3],
2126  &coef[indx][4], &coef[indx][5], &coef[indx][6],
2127  &tmp1, &tmp2, &tmp3);
2128  } else {
2129  /* don't have extra terms for avhrr nlsst */
2130  sscanf(line, "%4s %d %f %f %f %f %f %f", mission2,
2131  &month, &bounds[indx][0], &bounds[indx][1],
2132  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2133  &coef[indx][3]);
2134  coef[indx][4] = 0.0;
2135  coef[indx][5] = 0.0;
2136  coef[indx][6] = 0.0;
2137  }
2138  if (month == tmonth + 1) {
2139  indx++;
2140  if (indx == 7){
2141  found = 1;
2142  break;
2143  }
2144  }
2145  }
2146  }
2147  }
2148  }
2149 
2150  fclose(fp);
2151 
2152  if (found == 1) {
2153 
2154  if (sensorID == VIIRSN
2155  && (gotsstrefoffday == 0 || gotsstrefoffnight == 0)) {
2156  fprintf(stderr,
2157  "-E- %s line %d: Day and night sst reference offsets not found in %s\n",
2158  __FILE__, __LINE__, input->sstcoeffile);
2159  exit(1);
2160  }
2161 
2162  printf("Loading SST lat band coefficients from %s:\n",
2163  input->sstcoeffile);
2164 
2165  if (sensorID == VIIRSN && input->viirsnosisaf == 1) {
2166  /* viirs osi-saf coefs are not latband */
2167  for (indx = 0; indx < 4; indx++) {
2168  printf("%s %s %s %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2169  coeflabel[indx], sdates, edates, coef[indx][0], coef[indx][1], coef[indx][2],
2170  coef[indx][3], coef[indx][4], coef[indx][5], coef[indx][6]);
2171  }
2172  } else {
2173  for (indx = 0; indx < 7; indx++) {
2174  if (sensorID == VIIRSN && input->viirsnv7 >= 1) {
2175  /* v7 viirs have start and end dates, not months */
2176  printf(
2177  "%s %s %s %s %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2178  sdates, stime, edates, etime, bounds[indx][0], bounds[indx][1],
2179  coef[indx][0], coef[indx][1], coef[indx][2], coef[indx][3],
2180  coef[indx][4], coef[indx][5], coef[indx][6]);
2181  } else {
2182  /* most latband coeffs are by month */
2183  printf(
2184  "%d %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2185  month, bounds[indx][0], bounds[indx][1], coef[indx][0],
2186  coef[indx][1], coef[indx][2], coef[indx][3],
2187  coef[indx][4], coef[indx][5], coef[indx][6]);
2188  }
2189  }
2190  }
2191  printf(" sst reference day offset = %f\n", *sstrefoffday);
2192  printf(" sst reference night offset = %f\n", *sstrefoffnight);
2193 
2194  } else {
2195  fprintf(stderr,
2196  "-E- %s line %d: unable to locate valid SST coefficients for %s in %s\n",
2197  __FILE__, __LINE__, odates, input->sstcoeffile);
2198  exit(1);
2199  }
2200 
2201  return;
2202 }
2203 
2204 /* ----------------------------------------------------------------------------------- */
2205 /* read_sst4_coeff() - reads 4um sst coefficients for the specified date */
2206 /* */
2207 /* B. Franz, SAIC, 11 August 2005 */
2208 
2209 /* ----------------------------------------------------------------------------------- */
2210 void read_sst4_coeff(int32 sensorID, char *filename, double scantime,
2211  float **bounds, float **coef) {
2212  char mission[5];
2213  char mission2[5];
2214 
2215  FILE *fp;
2216  char *line;
2217  line = (char *) calloc(200, sizeof (char));
2218  char odate[8];
2219  int found = 0;
2220  // int tmonth;
2221  int16_t year, day;
2222  double sec;
2223  unix2yds(scantime, &year, &day, &sec);
2224 
2225  int month;
2226  int indx;
2227  // char tmp1[200] = "\0"; /* random length to slurp up whatevers after the coefs */
2228  int32 leap;
2229 
2230  switch (sensorID) {
2231  case MODISA:
2232  strcpy(mission, "AQUA");
2233  break;
2234  case MODIST:
2235  strcpy(mission, "TERR");
2236  break;
2237  case VIIRSN:
2238  strcpy(mission, "VIIR");
2239  break;
2240  default:
2241  strcpy(mission, "AQUA");
2242  break;
2243  }
2244 
2245  /* Open the file */
2246  if ((fp = fopen(filename, "r")) == NULL) {
2247  fprintf(stderr, "-E- %s line %d: unable to open %s for reading\n",
2248  __FILE__, __LINE__, filename);
2249  exit(1);
2250  }
2251 
2252  sprintf(odate, "%4d%03d", year, day);
2253 
2254  /* Find month */
2255  leap = (isleap(year) == TRUE ? 1 : 0);
2256  if (tmonth < 0) {
2257  for (tmonth = 11; tmonth >= 0; tmonth--) {
2258  if (day > StartOfMonth[leap][tmonth]) {
2259  break;
2260  }
2261  }
2262  }
2263  /* Loop through to find 6 sets of coeffs for this month */
2264  indx = 0;
2265  fprintf(stderr, " looking for month %d mission %s\n", tmonth + 1, mission);
2266  while (fgets(line, 200, fp)) {
2267  /* read sst4 latband coeffs */
2268  if (strncmp(line, mission, 4) == 0) {
2269  sscanf(line, "%4s %d %f %f %f %f %f %f %f %f %f", &mission2[0],
2270  &month, &bounds[indx][0], &bounds[indx][1], &coef[indx][0],
2271  &coef[indx][1], &coef[indx][2], &coef[indx][3],
2272  &coef[indx][4], &coef[indx][5], &coef[indx][6]);
2273  if (month == tmonth + 1) {
2274  indx++;
2275  if (indx == 7) {
2276  found = 1;
2277  break;
2278  }
2279  }
2280  }
2281  }
2282 
2283  fclose(fp);
2284 
2285  if (found == 1) {
2286 
2287  printf("Loading SST4 lat band coefficients from %s:\n", filename);
2288  for (indx = 0; indx < 7; indx++) {
2289  printf("%d %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2290  month, bounds[indx][0], bounds[indx][1], coef[indx][0], coef[indx][1], coef[indx][2],
2291  coef[indx][3], coef[indx][4], coef[indx][5], coef[indx][6]);
2292  }
2293 
2294  } else {
2295  fprintf(stderr,
2296  "-E- %s line %d: unable to locate valid SST4 coefficients for %s in %s\n",
2297  __FILE__, __LINE__, odate, filename);
2298  exit(1);
2299  }
2300  free(line);
2301 
2302  return;
2303 }
2304 
2305 /* ----------------------------------------------------------------------------------- */
2306 /* sstcloud() - uses red band to test homogeneity in nx x ny box (RSMAS) */
2307 /* */
2308 /* B. Franz, SAIC, August 2005. */
2309 
2310 /* ----------------------------------------------------------------------------------- */
2311 int sstcloud(int32_t ip, int32_t nx, int32_t ny, float thresh) {
2312  extern l1qstr l1que;
2313 
2314  int32_t nscan = l1que.nq;
2315  int32_t npix = l1que.r[0].npix;
2316  int32_t is = nscan / 2;
2317  int32_t ip1, ip2;
2318  int32_t is1, is2;
2319  int32_t i, j;
2320  int32_t ipb;
2321  float rhot;
2322  int flag = 0;
2323  int cnt = 0;
2324  float maxv = BAD_FLT;
2325  float minv = -1.0 * BAD_FLT; /* initial minimum has to be large positive number */
2326 
2327  if (!haveRed)
2328  return (flag);
2329 
2330  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2331  if (nscan < ny) {
2332  printf(
2333  "-E- %s line %d: L1 queue size of %d is too small for requested homogeneity test of %d x %d.\n",
2334  __FILE__, __LINE__, l1que.nq, nx, ny);
2335  exit(1);
2336  }
2337 
2338  /* algorithm is only valid for daytime and in non glint area */
2339  if (l1que.r[is].solz[ip] >= solznight
2340  || l1que.r[is].glint_coef[ip] > glintmax)
2341  return (flag);
2342 
2343  /* compute pix/scan limits for the Row Of Interest (ROI) */
2344  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2345  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2346  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2347  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2348  // filt_dist = -nx / 2;
2349  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2350  // ipmod -= l1que.r[is].spix;
2351  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2352  //
2353  // filt_dist = nx / 2;
2354  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2355  // ipmod -= l1que.r[is].spix;
2356  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2357  // } else {
2358  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2359  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2360  // }
2361 
2362  /* compute max and min normalized reflectance in window */
2363  /* l1_hmodis_hdf.c converts visible bands, i.e. red, to radiance */
2364  /* so we need to convert it to reflectance */
2365  /* atmospheric transmittance and path radiance due to scattering and attenuation by
2366  * atmospheric gases and aerosols is a problem for ocean color, but not for sst
2367  * so, while we probably don't need the extra terms here, we definitely don't need
2368  * them for the cloud decision trees */
2369  for (j = is1; j <= is2; j++)
2370  for (i = ip1; i <= ip2; i++) {
2371  ipb = i * nbvis + ibred;
2372  if (l1que.r[j].Lt[ipb] > 0.0 && l1que.r[j].Lt[ipb] < satred) {
2373  rhot = PI * l1que.r[j].Lt[ipb] / l1que.r[j].Fo[ibred]
2374  / l1que.r[j].tg_sol[ipb] / l1que.r[j].tg_sen[ipb]
2375  / l1que.r[j].t_sol[ipb] / l1que.r[j].t_sen[ipb]
2376  / cos(l1que.r[j].solz[ip] / RADEG);
2377  maxv = MAX(maxv, rhot);
2378  minv = MIN(minv, rhot);
2379  cnt++;
2380  }
2381  }
2382 
2383  /* added for collect 5 (if all saturated but 1, then cloud) */
2384  if (cnt < 2 || (maxv - minv) > thresh) {
2385  flag = 1;
2386  }
2387 
2388  return (flag);
2389 }
2390 
2391 /* ----------------------------------------------------------------------------------- */
2392 /* rhoCboxstats() - test homogeneity in nx x ny box of rho_cirrus (RSMAS) */
2393 
2394 /* ----------------------------------------------------------------------------------- */
2395 int32_t rhoCboxstats(int32_t ip, int32_t nx, int32_t ny, statstr *statrec) {
2396  extern l1qstr l1que;
2397 
2398  static int32_t len = 0;
2399  int32_t nscan = l1que.nq;
2400  int32_t npix = l1que.r[0].npix;
2401  int32_t ip1, ip2;
2402  int32_t is1, is2;
2403  int32_t is, i, j, ix;
2404  float Bt;
2405  static float *x = NULL;
2406 
2407 
2408  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2409  if (nscan < ny) {
2410  printf("-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2411  __FILE__, __LINE__, l1que.nq, nx, ny);
2412  exit(1);
2413  }
2414 
2415  /* allocate sufficient workspace for the filter size */
2416  if (x == NULL || len < nx * ny) {
2417  len = nx*ny;
2418  if (x != NULL) free(x);
2419  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2420  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2421  __FILE__, __LINE__);
2422  exit(1);
2423  }
2424  }
2425 
2426  /* Compute queue scan limits for the Row Of Interest (ROI) */
2427  is = nscan / 2;
2428  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2429  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2430 
2431  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2432  /* for aggregated VIIRS, limits are aggregation zone dependent */
2433  // if( ( l1que.r[is].l1file->sensorID == VIIRS ) && ( l1que.r[is].scn_fmt == 0 ) ) {
2434  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2435  // filt_dist = -nx / 2;
2436  // viirs_pxcvt_agdel( ip_scan, filt_dist, &ipmod );
2437  // ipmod -= l1que.r[is].spix;
2438  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2439  //
2440  // filt_dist = nx / 2;
2441  // viirs_pxcvt_agdel( ip_scan, filt_dist, &ipmod );
2442  // ipmod -= l1que.r[is].spix;
2443  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2444  // } else {
2445  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2446  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2447  // }
2448 
2449  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2450  statrec->min = -1.0 * BAD_FLT;
2451  statrec->max = BAD_FLT;
2452  statrec->avg = BAD_FLT;
2453  statrec->med = BAD_FLT;
2454  statrec->cen = BAD_FLT;
2455  statrec->sqr = BAD_FLT;
2456  statrec->sd = BAD_FLT;
2457  statrec->cnt = 0;
2458  for (ix = 0; ix < len; ix++) {
2459  x[ix] = 0.0;
2460  }
2461 
2462  /* compute stats for window */
2463  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
2464  Bt = l1que.r[j].rho_cirrus[i];
2465  if (Bt > (BAD_FLT + 0.1)) {
2466  statrec->max = MAX(statrec->max, Bt);
2467  statrec->min = MIN(statrec->min, Bt);
2468  if (statrec->cnt == 0) {
2469  /* first value - start the summing for the average at zero */
2470  statrec->avg = 0.0;
2471  }
2472  statrec->avg += Bt;
2473  statrec->sqr += Bt * Bt;
2474  x[statrec->cnt] = Bt;
2475  statrec->cnt++;
2476  }
2477  }
2478  if (statrec->cnt > 2) {
2479  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2480  statrec->max = x[statrec->cnt - 1];
2481  statrec->min = x[0];
2482  if (statrec->cnt > 1) {
2483  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2484  if (statrec->sd > 0.0) {
2485  statrec->sd = sqrt(statrec->sd);
2486  } else {
2487  statrec->sd = BAD_FLT;
2488  }
2489  }
2490 
2491  statrec->med = x[statrec->cnt / 2];
2492  }
2493  if (statrec->cnt > 0)
2494  statrec->avg /= statrec->cnt;
2495 
2496  return (statrec->cnt);
2497 }
2498 
2499 /* ----------------------------------------------------------------------------------- */
2500 /* rhotboxstats() - test homogeneity in nx x ny box of rhot's (RSMAS) */
2501 /* */
2502 
2503 /* ----------------------------------------------------------------------------------- */
2504 int32_t rhotboxstats(int32_t ip, int ib, int nbands, int32_t nx, int32_t ny,
2505  statstr *statrec) {
2506  extern l1qstr l1que;
2507 
2508  static int32_t len = 0;
2509  int32_t nscan = l1que.nq;
2510  int32_t npix = l1que.r[0].npix;
2511  int32_t ip1, ip2;
2512  int32_t is1, is2;
2513  int32_t is, i, j, ix;
2514  int32_t ipb;
2515  // int32_t ip_scan, filt_dist, ipmod;
2516  float rhot;
2517  static float *x = NULL;
2518 
2519  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2520  if (nscan < ny) {
2521  printf(
2522  "-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2523  __FILE__, __LINE__, l1que.nq, nx, ny);
2524  exit(1);
2525  }
2526 
2527  /* allocate sufficient workspace for the filter size */
2528  if (x == NULL || len < nx * ny) {
2529  len = nx*ny;
2530  if (x != NULL) free(x);
2531  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2532  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2533  __FILE__, __LINE__);
2534  exit(1);
2535  }
2536  }
2537 
2538  /* Compute queue scan limits for the Row Of Interest (ROI) */
2539  is = nscan / 2;
2540  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2541  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2542 
2543  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2544  /* for aggregated VIIRS, limits are aggregation zone dependent */
2545  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2546  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2547  // filt_dist = -nx / 2;
2548  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2549  // ipmod -= l1que.r[is].spix;
2550  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2551  //
2552  // filt_dist = nx / 2;
2553  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2554  // ipmod -= l1que.r[is].spix;
2555  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2556  // } else {
2557  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2558  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2559  // }
2560 
2561  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2562  statrec->min = -1.0 * BAD_FLT;
2563  statrec->max = BAD_FLT;
2564  statrec->avg = BAD_FLT;
2565  statrec->med = BAD_FLT;
2566  statrec->cen = BAD_FLT;
2567  statrec->sqr = BAD_FLT;
2568  statrec->sd = BAD_FLT;
2569  statrec->cnt = 0;
2570  for (ix = 0; ix < len; ix++) {
2571  x[ix] = 0.0;
2572  }
2573 
2574  /* compute stats for window */
2575  for (j = is1; j <= is2; j++)
2576  for (i = ip1; i <= ip2; i++) {
2577  ipb = i * nbands + ib;
2578  /* saturated Lt is 1000.0 or sstbad=-32767 (class bad could be 0.0, but don't care?) */
2579  if (l1que.r[j].Lt[ipb] > sstbad + 1.0 && l1que.r[j].Lt[ipb] < 1000.0 - 1.0) {
2580  rhot = M_PI * l1que.r[j].Lt[ipb] / l1que.r[j].Fo[ib] / l1que.r[j].csolz[i];
2581  statrec->max = MAX(statrec->max, rhot);
2582  statrec->min = MIN(statrec->min, rhot);
2583  if (statrec->cnt == 0) {
2584  /* first value - start the summing for the average at zero */
2585  statrec->avg = 0.0;
2586  }
2587  statrec->avg += rhot;
2588  statrec->sqr += rhot * rhot;
2589  x[statrec->cnt] = rhot;
2590  statrec->cnt++;
2591  }
2592  }
2593  if (statrec->cnt > 2) {
2594  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2595  statrec->max = x[statrec->cnt - 1];
2596  statrec->min = x[0];
2597  if (statrec->cnt > 1) {
2598  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2599  if (statrec->sd > 0.0) {
2600  statrec->sd = sqrt(statrec->sd);
2601  } else {
2602  statrec->sd = BAD_FLT;
2603  }
2604  }
2605 
2606  statrec->med = x[statrec->cnt / 2];
2607  }
2608  if (statrec->cnt > 0)
2609  statrec->avg /= statrec->cnt;
2610 
2611  return (statrec->cnt);
2612 }
2613 
2614 /* ----------------------------------------------------------------------------------- */
2615 /* btboxstats() - test homogeneity in nx x ny box of Bt's (RSMAS) */
2616 /* */
2617 /* B. Franz, SAIC, August 2005. */
2618 
2619 /* ----------------------------------------------------------------------------------- */
2620 int32_t btboxstats(int32_t ip, int ib, int nbands, int32_t nx, int32_t ny,
2621  statstr *statrec) {
2622  extern l1qstr l1que;
2623 
2624  static int32_t len = 0;
2625  int32_t nscan = l1que.nq;
2626  int32_t npix = l1que.r[0].npix;
2627  int32_t ip1, ip2;
2628  int32_t is1, is2;
2629  int32_t is, i, j, ix;
2630  int32_t ipb;
2631  float Bt;
2632  static float *x = NULL;
2633 
2634  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2635  if (nscan < ny) {
2636  printf(
2637  "-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2638  __FILE__, __LINE__, l1que.nq, nx, ny);
2639  exit(1);
2640  }
2641 
2642  /* allocate sufficient workspace for the filter size */
2643  if (x == NULL || len < nx * ny) {
2644  len = nx*ny;
2645  if (x != NULL) free(x);
2646  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2647  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2648  __FILE__, __LINE__);
2649  exit(1);
2650  }
2651  }
2652 
2653  /* Compute queue scan limits for the Row Of Interest (ROI) */
2654  is = nscan / 2;
2655  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2656  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2657 
2658  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2659  /* for aggregated VIIRS, limits are aggregation zone dependent */
2660  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2661  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2662  // filt_dist = -nx / 2;
2663  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2664  // ipmod -= l1que.r[is].spix;
2665  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2666  //
2667  // filt_dist = nx / 2;
2668  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2669  // ipmod -= l1que.r[is].spix;
2670  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2671  // } else {
2672  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2673  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2674  // }
2675 
2676  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2677  statrec->min = -1.0 * BAD_FLT;
2678  statrec->max = BAD_FLT;
2679  statrec->avg = BAD_FLT;
2680  statrec->med = BAD_FLT;
2681  statrec->cen = BAD_FLT;
2682  statrec->sqr = BAD_FLT;
2683  statrec->sd = BAD_FLT;
2684  statrec->cnt = 0;
2685  for (ix = 0; ix < len; ix++) {
2686  x[ix] = 0.0;
2687  }
2688 
2689  /* compute stats for window */
2690  for (j = is1; j <= is2; j++) {
2691  /* don't use values from aqua detector zero Bt40 */
2692  /* Peter and Kay say to not use average, just skip detector zero */
2693  if ((l1que.r[j].l1file->sensorID != MODISA)
2694  || (ib != ib40)
2695  || (l1que.r[j].detnum != 0)) {
2696 
2697  for (i = ip1; i <= ip2; i++) {
2698  ipb = i * nbands + ib;
2699  Bt = l1que.r[j].Bt[ipb];
2700  if (Bt > BT_LO + 0.1 && Bt < BT_HI - 0.1) {
2701  statrec->max = MAX(statrec->max, Bt);
2702  statrec->min = MIN(statrec->min, Bt);
2703  if (statrec->cnt == 0) {
2704  /* first value - start the summing for the average at zero */
2705  statrec->avg = 0.0;
2706  statrec->sqr = 0.0;
2707  }
2708  statrec->avg += Bt;
2709  statrec->sqr += Bt * Bt;
2710  x[statrec->cnt] = Bt;
2711  statrec->cnt++;
2712  }
2713  }
2714  }
2715  }
2716  if (statrec->cnt > 2) {
2717  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2718  statrec->max = x[statrec->cnt - 1];
2719  statrec->min = x[0];
2720  if (statrec->cnt > 1) {
2721  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2722  if (statrec->sd > 0.0) {
2723  statrec->sd = sqrt(statrec->sd);
2724  } else {
2725  statrec->sd = BAD_FLT;
2726  }
2727  }
2728 
2729  statrec->med = x[statrec->cnt / 2];
2730  }
2731  if (statrec->cnt > 0) {
2732  statrec->avg /= statrec->cnt;
2733  }
2734 
2735  return (statrec->cnt);
2736 }
2737 
2738 /* ----------------------------------------------------------------------------------- */
2739 /* sstboxstats() - test homogeneity in nx x ny box of Bt's (RSMAS) */
2740 /* */
2741 /* B. Franz, SAIC, August 2005. */
2742 
2743 /* ----------------------------------------------------------------------------------- */
2744 int32_t sstboxstats(int32_t ip, int32_t nx, int32_t ny, statstr *statrec) {
2745  extern l1qstr l1que;
2746 
2747  static int32_t len = 0;
2748  int32_t nscan = l1que.nq;
2749  int32_t npix = l1que.r[0].npix;
2750  int32_t ip1, ip2;
2751  int32_t is1, is2;
2752  int32_t is, i, j, ix;
2753  // int32_t ip_scan, filt_dist, ipmod;
2754  float sst;
2755  static float *x = NULL;
2756 
2757  /* make sure code is not inconsistent */
2758  if (nscan < ny) {
2759  printf(
2760  "-E- %s line %d: sst queue size of %d is too small for %d x %d box stats.\n",
2761  __FILE__, __LINE__, nscan, nx, ny);
2762  exit(1);
2763  }
2764 
2765  /* allocate sufficient workspace for the filter size */
2766  if (x == NULL || len < nx * ny) {
2767  len = nx*ny;
2768  if (x != NULL) free(x);
2769  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2770  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2771  __FILE__, __LINE__);
2772  exit(1);
2773  }
2774  }
2775 
2776  /* Compute queue scan limits for the Row Of Interest (ROI) */
2777  is = nscan / 2;
2778  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2779  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2780 
2781  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2782  /* for aggregated VIIRS, limits are aggregation zone dependent */
2783  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2784  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2785  // filt_dist = -nx / 2;
2786  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2787  // ipmod -= l1que.r[is].spix;
2788  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2789  //
2790  // filt_dist = nx / 2;
2791  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2792  // ipmod -= l1que.r[is].spix;
2793  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2794  // } else {
2795  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2796  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2797  // }
2798 
2799  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2800  statrec->min = -1.0 * BAD_FLT;
2801  statrec->max = BAD_FLT;
2802  statrec->avg = BAD_FLT;
2803  statrec->med = BAD_FLT;
2804  statrec->cen = BAD_FLT;
2805  statrec->sqr = BAD_FLT;
2806  statrec->sd = BAD_FLT;
2807  statrec->cnt = 0;
2808  for (ix = 0; ix < len; ix++) {
2809  x[ix] = 0.0;
2810  }
2811 
2812  /* compute stats for window */
2813  for (j = is1; j <= is2; j++) {
2814  for (i = ip1; i <= ip2; i++) {
2815  sst = sstq[j][i];
2816  if (sst >= SSTmin && sst <= SSTmax) {
2817  statrec->max = MAX(statrec->max, sst);
2818  statrec->min = MIN(statrec->min, sst);
2819  if (statrec->cnt == 0) {
2820  /* first value - start the summing for the average at zero */
2821  statrec->avg = 0.0;
2822  statrec->sqr = 0.0;
2823  }
2824  statrec->avg += sst;
2825  statrec->sqr += sst * sst;
2826  x[statrec->cnt] = sst;
2827  statrec->cnt++;
2828  }
2829  }
2830  }
2831  if (statrec->cnt > 2) {
2832  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2833  statrec->max = x[statrec->cnt - 1];
2834  statrec->min = x[0];
2835  if (statrec->cnt > 1) {
2836  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2837  if (statrec->sd > 0.0) {
2838  statrec->sd = sqrt(statrec->sd);
2839  } else {
2840  statrec->sd = BAD_FLT;
2841  }
2842  }
2843 
2844  statrec->med = x[statrec->cnt / 2];
2845  }
2846  if (statrec->cnt > 0) {
2847  statrec->avg /= statrec->cnt;
2848  }
2849 
2850  return (statrec->cnt);
2851 }
2852 /* ----------------------------------------------------------------------------------- */
2853 /* run_sstboxstdev() - calculate the sst standard deviation */
2854 
2855 /* ----------------------------------------------------------------------------------- */
2856 void run_sstboxstdev(int npix, float *sst_stdev) {
2857 
2858  int ip;
2859  statstr statrec;
2860  for (ip = 0; ip < npix; ip++) {
2861  sst_stdev[ip] = sstbad;
2862  if (sstboxstats(ip, btbox, btbox, &statrec) > 0) {
2863  sst_stdev[ip] = statrec.sd;
2864  }
2865  }
2866 }
2867 
2868 /* ----------------------------------------------------------------------------------- */
2869 /* run_rhotboxmin() - initialize the static min arrays */
2870 
2871 /* ----------------------------------------------------------------------------------- */
2872 void run_rhotboxmin(int npix, int btidx, int nbands, float *minarr) {
2873 
2874  int ip;
2875  statstr statrec;
2876  for (ip = 0; ip < npix; ip++) {
2877  minarr[ip] = sstbad;
2878  if (rhotboxstats(ip, btidx, nbands, btbox, btbox, &statrec) > 0) {
2879  minarr[ip] = statrec.min;
2880  }
2881  }
2882 }
2883 
2884 /* ----------------------------------------------------------------------------------- */
2885 /* run_rhotboxmaxmin() - initialize the static maxmin arrays */
2886 
2887 /* ----------------------------------------------------------------------------------- */
2888 void run_rhotboxmaxmin(int npix, int btidx, int nbands, float *maxminarr, float *maxarr) {
2889 
2890  int ip;
2891  statstr statrec;
2892  for (ip = 0; ip < npix; ip++) {
2893  maxminarr[ip] = sstbad;
2894  maxarr[ip] = sstbad;
2895  if (rhotboxstats(ip, btidx, nbands, btbox, btbox, &statrec) > 0) {
2896  maxminarr[ip] = statrec.max - statrec.min;
2897  maxarr[ip] = statrec.max;
2898  }
2899  }
2900 }
2901 
2902 /* ----------------------------------------------------------------------------------- */
2903 /* run_btboxavg() - initialize the static avg arrays */
2904 
2905 /* ----------------------------------------------------------------------------------- */
2906 void run_btboxavg(int npix, int btidx, float *avgarr) {
2907 
2908  int ip;
2909  statstr statrec;
2910  for (ip = 0; ip < npix; ip++) {
2911  avgarr[ip] = sstbad;
2912  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2913  avgarr[ip] = statrec.avg;
2914  }
2915  }
2916 }
2917 
2918 /* ----------------------------------------------------------------------------------- */
2919 /* run_btboxmaxmin() - initialize the static maxmin arrays */
2920 
2921 /* ----------------------------------------------------------------------------------- */
2922 void run_btboxmaxmin(int npix, int btidx, float *maxminarr) {
2923 
2924  int ip;
2925  statstr statrec;
2926  for (ip = 0; ip < npix; ip++) {
2927  maxminarr[ip] = sstbad;
2928  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2929  maxminarr[ip] = statrec.max - statrec.min;
2930  }
2931  }
2932 }
2933 /* ----------------------------------------------------------------------------------- */
2934 /* run_btboxmin() - initialize the static min arrays */
2935 
2936 /* ----------------------------------------------------------------------------------- */
2937 void run_btboxmin(int npix, int btidx, float *minarr) {
2938 
2939  int ip;
2940  statstr statrec;
2941  for (ip = 0; ip < npix; ip++) {
2942  minarr[ip] = sstbad;
2943  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2944  minarr[ip] = statrec.min;
2945  }
2946  }
2947 }
2948 /* ----------------------------------------------------------------------------------- */
2949 /* run_btboxmax() - initialize the static max arrays */
2950 
2951 /* ----------------------------------------------------------------------------------- */
2952 void run_btboxmax(int npix, int btidx, float *maxarr) {
2953 
2954  int ip;
2955  statstr statrec;
2956  for (ip = 0; ip < npix; ip++) {
2957  maxarr[ip] = sstbad;
2958  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2959  maxarr[ip] = statrec.max;
2960  }
2961  }
2962 }
2963 
2964 void run_rhoCboxmaxmin(int npix, float *maxminarr, float *minarr, float *maxarr) {
2965 
2966  int ip;
2967  statstr statrec;
2968  for (ip = 0; ip < npix; ip++) {
2969  maxminarr[ip] = sstbad;
2970  minarr[ip] = sstbad;
2971  maxarr[ip] = sstbad;
2972  if (rhoCboxstats(ip, btbox, btbox, &statrec) > 0) {
2973  maxminarr[ip] = statrec.max - statrec.min;
2974  minarr[ip] = statrec.min;
2975  maxarr[ip] = statrec.max;
2976  }
2977  }
2978 }
2979 
2980 void run_btboxstdev(int npix, int btidx, float *stdevarr) {
2981  int ip;
2982  statstr statrec;
2983  for (ip = 0; ip < npix; ip++) {
2984  stdevarr[ip] = sstbad;
2985  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2986  stdevarr[ip] = statrec.sd;
2987  }
2988  }
2989 }
2990 /* ----------------------------------------------------------------------------------- */
2991 /* btavg() - Average detector 9 with detector 1 to replace aqua channel 23 detector 0 */
2992 /* */
2993 /* S. Walsh, RSMAS, August 2013. */
2994 
2995 /* ----------------------------------------------------------------------------------- */
2996 int32_t btavg(int32_t ip, int is, int ib, int nbands, statstr *statrec) {
2997  extern l1qstr l1que;
2998 
2999  int32_t nscan = l1que.nq;
3000  int32_t is1, is2;
3001  int32_t j;
3002  int32_t ipb;
3003  float Bt;
3004 
3005  /* this routine should only be used to fake Bt40's for aqua detector zero */
3006  if ((l1que.r[is].l1file->sensorID != MODISA)
3007  || (ib != ib40)
3008  || (l1que.r[is].detnum != 0)) {
3009  printf(
3010  "-E- %s line %d: btavg should only be used for AQUA detector zero Bt40.\n",
3011  __FILE__, __LINE__);
3012  exit(1);
3013  }
3014 
3015  /* make sure code is not inconsistent NQMIN in l12_parms.h */
3016  /* Right now modisa uses 3x3 boxes for aqua.
3017  * When we switch to 5x5 then modisa/msl12_filter.dat will have to change to keep a 7x7 box
3018  * so that we have the line before and after the 5x5 box to calculate the Bt40 average for
3019  * detector 0 in whichever line of the box is being calculated
3020  */
3021  if ((nscan < (btbox + 2)) && (is == 0 || is == btbox)) {
3022  printf(
3023  "-E- %s line %d: L1 queue size of %d is too small for 3 line Bt40 average around line %d.\n",
3024  __FILE__, __LINE__, l1que.nq, is);
3025  exit(1);
3026  }
3027 
3028  /* Compute queue scan limits for the Row Of Interest (ROI) in 'is' */
3029  is1 = MIN(MAX(0, is - 1), nscan - 1);
3030  is2 = MAX(MIN(nscan - 1, is + 1), 0);
3031 
3032  /* initialize output rec - only min,max,avg,cnt are actually used here*/
3033  statrec->min = -1.0 * BAD_FLT;
3034  statrec->max = BAD_FLT;
3035  statrec->avg = BAD_FLT;
3036  statrec->med = BAD_FLT;
3037  statrec->cen = BAD_FLT;
3038  statrec->sqr = BAD_FLT;
3039  statrec->sd = BAD_FLT;
3040  statrec->cnt = 0;
3041 
3042  /* compute stats for window */
3043  for (j = is1; j <= is2; j += 2) {
3044  /* don't use bad middle line - aqua detector zero channel 23 */
3045  if (j != is) {
3046  ipb = ip * nbands + ib;
3047  Bt = l1que.r[j].Bt[ipb];
3048  if (Bt > BT_LO + 0.1 && Bt < BT_HI - 0.1) {
3049  if (statrec->cnt == 0) {
3050  /* first value - start the summing for the average at zero */
3051  statrec->avg = 0.0;
3052  }
3053  statrec->avg += Bt;
3054  statrec->cnt++;
3055  }
3056  }
3057  }
3058  if (statrec->cnt > 0)
3059  statrec->avg /= statrec->cnt;
3060 
3061  return (statrec->cnt);
3062 }
3063 
3064 /* ----------------------------------------------------------------------------------- */
3065 /* sstmasked() - returns 1 if pixel was already masked (SST processing skipped) */
3066 
3067 /* ----------------------------------------------------------------------------------- */
3068 int sstmasked(int32_t *flags, int32_t ip) {
3069  if ((flags[ip] & LAND) != 0 || (flags[ip] & NAVFAIL) != 0)
3070 
3071  return (1);
3072  else
3073  return (0);
3074 }
3075 
3076 /* ------------------------------------------------------------------- */
3077 /* avhrr_ascend() - compare center lat of this line and previous line */
3078 
3079 /* ------------------------------------------------------------------- */
3080 int32_t avhrr_ascend(int32_t ny, float *diflat) {
3081  extern l1qstr l1que;
3082 
3083  int32_t nscan = l1que.nq;
3084  int32_t is, cpix;
3085  int32_t jj;
3086  //float diflat;
3087  float diflats[FILTMAX];
3088 
3089  /* make sure code is not inconsistent NQMIN in l12_parms.h */
3090  if (nscan < ny) {
3091  printf(
3092  "-E- %s line %d: L1 queue size of %d is too small for %d line box stats.\n",
3093  __FILE__, __LINE__, l1que.nq, ny);
3094  exit(-1); /* exit value changed from 1 to -1 */
3095  }
3096 
3097  /* current row is center of box */
3098  is = nscan / 2; /* is =1 for boxsize (nscan=ny=) 3 */
3099 
3100  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
3101  /* 204/409 or 1023/2048 or 1599/3200 (zero based) */
3102  /* pixel number start at 0 so first half of scan lines are:
3103  * 0..204, 0..1023, 0..1599 */
3104  cpix = (fullscanpix - 1) / 2;
3105  // if (l1que.r[is - 1].lat[cpix] < l1que.r[is].lat[cpix]) {
3106  // }
3107  // I don't know if I ever checked this for avhrr, but modis and viirs seem to wobble near the north pole (I didn't check south)
3108  // Their nadir lat's don't always increase when approaching the north pole, sometimes they decrease from one scan to the next.
3109  // so check for average difference between nadir lats to find split between ascending and descending
3110  // average didn't work. in V2016183000600, which ascends near the dateline and descends over Russia,
3111  // every 16th line (the difference between detector 15 and detector 0) is a greater change in latitude
3112  // than the difference between the lines in one scan. Hard to describe, but easy to see if plot clat[2900:3246]
3113  // Try assuming flip from ascending to descending (and reverse in the south) is only when two in a row match.
3114  *diflat = 0.0;
3115  //*diflat=*diflat+l1que.r[jj].lat[cpix] - l1que.r[jj+1].lat[cpix];
3116  //*diflat=*diflat/(nscan-1);
3117  // check the diff in clat between the lines in the box
3118  for (jj = 0; jj < nscan - 1; jj++) {
3119  diflats[jj] = l1que.r[jj].lat[cpix] - l1que.r[jj + 1].lat[cpix];
3120  }
3121  // check for no change in direction
3122  if ((diflats[is - 1] > 0.0 && diflats[is] > 0.0) ||
3123  (diflats[is - 1] < 0.0 && diflats[is] < 0.0)) {
3124  *diflat = diflats[is];
3125  } else {
3126  // there was a change. But really only if the next dif matches the same direction, otherwise it's just the flaky
3127  // change between swaths (detector 15 to detector 0)
3128  if ((diflats[is] > 0.0 && diflats[is + 1] > 0.0) ||
3129  (diflats[is] < 0.0 && diflats[is + 1] < 0.0)) {
3130  *diflat = diflats[is];
3131  } else {
3132  // not really a change, return the previous diff (or next, but not this one)
3133  *diflat = diflats[is - 1];
3134  }
3135  }
3136 
3137  //if (l1que.cscan == 156) {
3138  // printf(" in ascend, diflat avg = %f\n", *diflat);
3139  //}
3140  if (*diflat < 0.0) {
3141  /* satellite is ascending */
3142  return (1);
3143  } else {
3144  /* satellite is not ascending (usually descending) */
3145  return (0);
3146  }
3147 }
3148 
3149 /* ----------------------------------------------------------------------------------- */
3150 /* set_flags_sst() - set quality flags for long wave sea surface temperature */
3151 /* */
3152 /* B. Franz, SAIC, August 2005. */
3153 
3154 /* ----------------------------------------------------------------------------------- */
3155 void set_flags_sst(l2str *l2rec) {
3156 
3157  extern l1qstr l1que;
3158  int32_t npix = l2rec->l1rec->npix;
3159  int32_t cpix;
3160  int32_t ip, ipb, ipbir;
3161  int32_t yyyyddd;
3162  int32_t ASCEND;
3163  int32_t v5viirs;
3164  float LtRED;
3165  float LtNIR7; // change m6 to LtNIR7, as it's VIIRS 746nm;
3166  float rhoCirrus; // change m9 to rhoCirrus (VIIRS 1.38um);
3167  float Lt16; // change m10 to Lt16, as it's VIIRS 1.61 micron band;
3168  float LtNIR8;
3169  float Bt37;
3170  float Bt39;
3171  float Bt40;
3172  float Bt67;
3173  float Bt85; /* aqua trees - don't use 67, 73, 85 because of cross talk, until fixed 2016? */
3174  float Bt11;
3175  float Bt12;
3176  float dSST_ref; /* sst - ref */
3177  float dSST_SST3; /* sst - sst3 */
3178  float dSST_SST4; /* sst - sst4 */
3179  float dBt_11_12, dBt_37_11, dBt_37_12, dBt_11_37, dBt_67_11; //brightness temp diffs
3180  float dBt_40_11, dBt_85_11;
3181  float Tdeflong;
3182  float subsolar, xdoy, xrad;
3183  float diflat;
3184  statstr statrec;
3185 
3186  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
3187  cpix = (fullscanpix - 1) / 2; /* 204/409 or 1023/2048 (zero based) */
3188 
3189  get_toa_refl(l2rec, ibred, rhotRED);
3190  get_toa_refl(l2rec, ib07, rhotNIR7);
3191  get_toa_refl(l2rec, ib16, rhot16);
3192 
3193  /* check each pixel in scan */
3194  for (ip = 0; ip < npix; ip++) {
3195 
3196  /* SST not processed */
3197  if (sstmasked(l2rec->l1rec->flags, ip)) {
3198  flags_sst[ip] |= SSTF_ISMASKED;
3199  continue;
3200  }
3201 
3202  ipb = ip * l2rec->l1rec->l1file->nbands;
3203  ipbir = ip * NBANDSIR;
3204 
3205  treesum[ip] = 0.0;
3206 
3207  Bt37 = l2rec->l1rec->Bt[ipbir + ib37]; /* modis chan 20, avhrr cen3, viirs SVM12 */
3208  Bt85 = l2rec->l1rec->Bt[ipbir + ib85]; /* modis chan 29, viirs SVM14 */
3209  Bt11 = l2rec->l1rec->Bt[ipbir + ib11]; /* modis chan 31, avhrr cen4, viirs SVM15 */
3210  Bt12 = l2rec->l1rec->Bt[ipbir + ib12]; /* modis chan 32, avhrr cen5, viirs SVM16 */
3211 
3212  LtRED = l2rec->l1rec->Lt[ipb + ibred];
3213 
3214  // some bands really are sensor specific...
3215  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
3216  LtNIR8 = l2rec->l1rec->Lt[ipb + ib08];
3217  } else {
3218  rhoCirrus = l2rec->l1rec->rho_cirrus[ip];
3219  if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
3220  Bt39 = l2rec->l1rec->Bt[ipbir + ib39];
3221  /* don't use aqua Bt40 detector 0, average the previous and next scans instead */
3222  if (l2rec->l1rec->l1file->sensorID == MODISA && l2rec->l1rec->detnum == 0) {
3223  //Bt40 = Bt40_avg[ip];
3224  /* current row is center of box */
3225  int32_t is = l1que.nq / 2;
3226  if (btavg(ip, is, ib40, NBANDSIR, &statrec) > 0) {
3227  Bt40 = statrec.avg;
3228  }
3229  } else
3230  Bt40 = l2rec->l1rec->Bt[ipbir + ib40];
3231  Bt67 = l2rec->l1rec->Bt[ipbir + ib67];
3232  } else if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
3233  Bt40 = l2rec->l1rec->Bt[ipbir + ib40];
3234  LtNIR7 = l2rec->l1rec->Lt[ipb + ib07];
3235  Lt16 = l2rec->l1rec->Lt[ipb + ib16];
3236  }
3237  }
3238 
3239  /* BT could not be computed (radiance out-of-range) */
3240  if (Bt11 < BT_LO + 0.1 || Bt11 > BT_HI - 0.1 || Bt12 < BT_LO + 0.1
3241  || Bt12 > BT_HI - 0.1) {
3242  flags_sst[ip] |= SSTF_BTBAD;
3243  continue;
3244  }
3245 
3246  /* check BT range */
3247  if (Bt11 < Btmin || Bt11 > Btmax || Bt12 < Btmin || Bt12 > Btmax)
3248  flags_sst[ip] |= SSTF_BTRANGE;
3249 
3250  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
3251  if ((ASCEND = avhrr_ascend(btbox, &diflat)) == 1) {
3252  flags_sst[ip] |= SSTF_ASCEND;
3253  }
3254  /* if day and glint > thresh set bit M2B8 (shared with BT4REFDIFF) */
3255  if (l2rec->l1rec->solz[ip] < solznight) {
3256  if (l2rec->l1rec->glint_coef[ip] > glintmax)
3257  flags_sst[ip] |= SSTF_GLINT;
3258  }
3259 
3260  /* ------------------------------------------------------------------------- */
3261  /* Gross Cloud contamination test for tree test BTRANGE (M1B1) */
3262  /* Based on Kilpatrick et al: */
3263  /* Overfiew of NOAA/NASA AVHRR Pathfinder Algorithm. Figure 6. */
3264  /* ------------------------------------------------------------------------- */
3265  /* check BT range: BT must be ge -10 C and le 35 C, M1B1 */
3266  /* some noaa satellites need channel 3, and some only at night */
3267  /* if n7-15 at night or at day and not glint or
3268  * if n16-19 at night */
3269  /* Nov 2012, day only check too cold except older avhrr need too warm check also */
3270 
3271  /* if older avhrr sat and night or day and not glint check for Bt37 too warm or too cold */
3272  if ((l2rec->l1rec->l1file->subsensorID == NO07
3273  || l2rec->l1rec->l1file->subsensorID == NO09
3274  || l2rec->l1rec->l1file->subsensorID == NO10
3275  || l2rec->l1rec->l1file->subsensorID == NO11
3276  || l2rec->l1rec->l1file->subsensorID == NO12
3277  || l2rec->l1rec->l1file->subsensorID == NO14
3278  || l2rec->l1rec->l1file->subsensorID == NO15)
3279  && (l2rec->l1rec->solz[ip] >= solznight
3280  || (l2rec->l1rec->solz[ip] < solznight
3281  && l2rec->l1rec->glint_coef[ip] <= glintmax))) {
3282  if (Bt37 < Btmin || Bt37 > Btmax) {
3283  flags_sst[ip] |= SSTF_BTRANGE;
3284  }
3285  }
3286 
3287  /* if n16-19 Bt37 at night (channel 3 is a different wavelength during the day) is too warm or too cold */
3288  if ((l2rec->l1rec->l1file->subsensorID == NO16
3289  || l2rec->l1rec->l1file->subsensorID == NO17
3290  || l2rec->l1rec->l1file->subsensorID == NO18
3291  || l2rec->l1rec->l1file->subsensorID == NO19)
3292  && l2rec->l1rec->solz[ip] >= solznight) {
3293  if (Bt37 < Btmin || Bt37 > Btmax) {
3294  flags_sst[ip] |= SSTF_BTRANGE;
3295  }
3296  }
3297  } else {
3298  /* for viirs (and modis), to test l2bin dataday splitting and not mixing ascending and descending, set spare flag */
3299  if ((ASCEND = avhrr_ascend(btbox, &diflat)) == 1) {
3300 // l2rec->l1rec->flags[ip] |= ASCFLG; //NEED THIS FLAG DEFINITION FROM SOMEWHERE
3301  }
3302  //if (ip == cpix) {
3303  // printf(" pix=%d line=%d diflat=%f\n",ip,l2rec->l1rec->iscan,diflat);
3304  //}
3305  }
3306  if ((l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA)) {
3307  /* if modis Bt37 or Bt39 or Bt40 at night is too warm or too cold */
3308 
3309  if (l2rec->l1rec->solz[ip] >= solznight
3310  && (Bt37 < Btmin || Bt37 > Btmax || Bt39 < Btmin
3311  || Bt39 > Btmax || Bt40 < Btmin || Bt40 > Btmax40)) {
3312  flags_sst[ip] |= SSTF_BTRANGE;
3313  }
3314  /* if modis Bt37 or Bt39 or Bt40 is during the day too cold */
3315  if ((l2rec->l1rec->solz[ip] < solznight
3316  && l2rec->l1rec->glint_coef[ip] <= glintmax)
3317  && (Bt37 < Btmin || Bt39 < Btmin || Bt40 < Btmin)) {
3318  flags_sst[ip] |= SSTF_BTRANGE;
3319  }
3320  }
3321 
3322  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
3323  /* if viirs Bt37 or Bt40 at night is too warm or too cold */
3324  if (l2rec->l1rec->solz[ip] >= solznight
3325  && (Bt37 < Btmin || Bt37 > Btmax || Bt40 < Btmin
3326  || Bt40 > Btmax40)) {
3327  flags_sst[ip] |= SSTF_BTRANGE;
3328  }
3329 
3330  /* if viirs Bt37 or Bt40 during the day is too cold */
3331  if ((l2rec->l1rec->solz[ip] < solznight
3332  && l2rec->l1rec->glint_coef[ip] <= glintmax)
3333  && (Bt37 < Btmin || Bt40 < Btmin)) {
3334  flags_sst[ip] |= SSTF_BTRANGE;
3335  }
3336  }
3337 
3338  /* check BT diff */
3339  dBt_11_12 = Bt11 - Bt12;
3340  if (dBt_11_12 < dBtmin || dBt_11_12 > dBtmax)
3341  flags_sst[ip] |= SSTF_BTDIFF;
3342 
3343  /* check SST range, using different max limits for day and night */
3344  if ((sstq[csstbox][ip] < SSTmin)
3345  || (l2rec->l1rec->solz[ip] < solznight && sstq[csstbox][ip] > SSTmax)
3346  || (l2rec->l1rec->solz[ip] >= solznight && sstq[csstbox][ip] > SSTmaxn)) {
3347  flags_sst[ip] |= SSTF_SSTRANGE;
3348  }
3349 
3350  /* check SST difference with references */
3351  /* check for both bad, or don't bother checking anywhere?
3352  * if either are bad then SSTF_SSTRANGE will be set so
3353  * no other flags matter? */
3354  dSST_ref = sstq[csstbox][ip] - l2rec->l1rec->sstref[ip];
3355  // the "coldonly" and equatorial aerosol tests are to be run by default, but if SSTMODS is set, don't
3356  if ((evalmask & SSTMODS) == 0) {
3357  /* evaluate change to cold-test only */
3358  /* set the flag bit if sst is too much colder than reference */
3359  if (dSST_ref < -SSTdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
3360  flags_sst[ip] |= SSTF_SSTREFDIFF;
3361  if (dSST_ref < -input->sstrefdif &&
3362  l2rec->l1rec->lat[ip] >= equatorialSouth && l2rec->l1rec->lat[ip] <= equatorialNorth &&
3363  l2rec->l1rec->lon[ip] >= equatorialWest && l2rec->l1rec->lon[ip] <= equatorialEast) {
3364  /* tighter test between 10S and 30N and -105 to 105 longitude */
3365  /* equatorial aerosol test */
3366  flags_sst[ip] |= SSTF_SSTREFDIFF;
3367  }
3368  /* set the flag bit if sst is too much warmer than reference at night */
3369  if (l2rec->l1rec->solz[ip] >= solznight) {
3370  if (fabs(dSST_ref) > SSTdiff)
3371  flags_sst[ip] |= SSTF_SSTREFDIFF;
3372  }
3373  /* is sst way too much colder than reference? */
3374  if (dSST_ref < -SSTvdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
3375  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3376  /* set the flag bit if sst is way too much warmer than reference at night */
3377  if (l2rec->l1rec->solz[ip] >= solznight) {
3378  if (fabs(dSST_ref) > SSTvdiff)
3379  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3380  }
3381  } else {
3382  if (l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3383  if (fabs(dSST_ref) > SSTdiff)
3384  flags_sst[ip] |= SSTF_SSTREFDIFF;
3385  } else {
3386  if (dSST_ref < -SSTdiff || dSST_ref > (SSTdiff + 1))
3387  flags_sst[ip] |= SSTF_SSTREFDIFF;
3388  }
3389  if (fabs(dSST_ref) > SSTvdiff)
3390  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3391  } //end of to-do-or-not-to-do coldonly block
3392 
3393  /* check SST difference with 4um SST only at night */
3394  /* set BT4REFDIFF flag based on sst4 flags */
3395  if (haveSST4) {
3396  dSST_SST4 = sstq[csstbox][ip] - sst4q[csstbox][ip];
3397  if (sst4q[csstbox][ip] > sstbad + 1.0 && l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3398  if (fabs(dSST_SST4) < SST4diff1)
3399  flags_sst[ip] |= SSTF_SST4DIFF;
3400  if (fabs(dSST_SST4) < SST4diff2)
3401  flags_sst[ip] |= SSTF_SST4VDIFF;
3402  flags_sst[ip] |= (flags_sst4[ip] & SSTF_BT4REFDIFF);
3403  }
3404  }
3405 
3406  /* for viirs, check SST difference with jpss triple window SST only at night */
3407  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
3408  dSST_SST3 = sstq[csstbox][ip]-sst3q[csstbox][ip];
3409  if (sst3q[csstbox][ip] > sstbad+1.0 && l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3410  if (fabs(dSST_SST3) > SST3diff1)
3411  flags_sst[ip] |= SSTF_SST3DIFF;
3412  if (fabs(dSST_SST3) > SST3diff2)
3413  flags_sst[ip] |= SSTF_SST3VDIFF;
3414  }
3415  }
3416 
3417  /* check sensor zenith limits */
3418  if (l2rec->l1rec->senz[ip] > hisenz)
3419  flags_sst[ip] |= SSTF_HISENZ;
3420  if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
3421  if (l2rec->l1rec->senz[ip] > vhisenzv2) {
3422  /* different limit only for viirs sst (not sst3) */
3423  flags_sst[ip] |= SSTF_VHISENZ;
3424  }
3425  } else {
3426  /* not viirs */
3427  if (l2rec->l1rec->senz[ip] > vhisenz)
3428  flags_sst[ip] |= SSTF_VHISENZ;
3429  }
3430 
3431  // flag 2 edge pixels as SSTF_VHISENZ so quality gets set to 3
3432  if (l2rec->l1rec->pixnum[ip] < 2 || l2rec->l1rec->pixnum[ip] > (fullscanpix - 3))
3433  flags_sst[ip] |= SSTF_VHISENZ;
3434  // set the last 4 pixels of the scan for Terra to VHISENZ
3435  // as there is an apparent calibration issue with the BT12 for those pixels
3436  if ((l2rec->l1rec->l1file->sensorID == MODIST) && (l2rec->l1rec->pixnum[ip] > 1349))
3437  flags_sst[ip] |= SSTF_VHISENZ;
3438 
3439  /* if SST is cold (collect 5) */
3440  /* I suppose it doesn't hurt to do this for AVHRR since it isn't used to determine quality? */
3441  /* do this for all satellites */
3442  /* sstcloud checks to make sure it's day and not glint */
3443  if (sstq[csstbox][ip] > sstbad + 1.0 && l2rec->l1rec->sstref[ip] > sstbad + 1.0
3444  && sstq[csstbox][ip] - l2rec->l1rec->sstref[ip] <= -1.0)
3445  if (sstcloud(ip, cldbox, cldbox, cldthresh) == 1)
3446  flags_sst[ip] |= SSTF_REDNONUNIF;
3447 
3448  /* check homogeneity of BT */
3449  if (Bt11_maxmin[ip] > Bt11unif1)
3450  flags_sst[ip] |= SSTF_BTNONUNIF;
3451  if (Bt11_maxmin[ip] > Bt11unif2)
3452  flags_sst[ip] |= SSTF_BTVNONUNIF;
3453 
3454  if (Bt12_maxmin[ip] > Bt12unif1)
3455  flags_sst[ip] |= SSTF_BTNONUNIF;
3456  if (Bt12_maxmin[ip] > Bt12unif2)
3457  flags_sst[ip] |= SSTF_BTVNONUNIF;
3458 
3459  /* end of homogeneity checks */
3460 
3461  dBt_37_11 = Bt37 - Bt11;
3462  dBt_37_12 = Bt37 - Bt12;
3463  dBt_11_37 = Bt11 - Bt37;
3464  dBt_67_11 = Bt67 - Bt11;
3465  dBt_40_11 = Bt40 - Bt11;
3466  dBt_85_11 = Bt85 - Bt11;
3467 
3468  /* values from Kay's VIIRS trees */
3469  /* Bt's are BT_LO = -1000 or BT_HI = 1000 when bad */
3470  /* If a Bt is bad then the pixel will be bad so it doesn't matter what happens in the trees */
3471  if (Bt11 != 0.0) {
3472  Tdeflong = (Bt11 - Bt12) / Bt11;
3473  } else {
3474  Tdeflong = (Bt11 - Bt12) / (Bt11 + 0.00001); /* don't divide by zero */
3475  }
3476  double pasutime = l2rec->l1rec->scantime;
3477  int16_t year, day;
3478  double sec;
3479  unix2yds(pasutime, &year, &day, &sec);
3480 
3481  yyyyddd = year * 1000 + day;
3482 
3483  /* --------------------------------------------------------------------------*/
3484  /* decision_tree() - Tree models are based on binary recursive partitioning */
3485  /* to indicate whether data are potentially contaminated. */
3486  /* based on charts sent to NODC from Kay Kilpatrick */
3487  /* worked on by Vicky Lin, SAIC, October 2007. */
3488  /* V6 trees from Guillermo Podesta. */
3489  /* ------------------------------------------------------------------------- */
3490  switch (l2rec->l1rec->l1file->subsensorID) {
3491  case NO07:
3492  /* use the noaa-9 tree for now (Kay: Mar 2013) */
3493  // break;
3494  case NO09:
3495  if (yyyyddd < 1994001) {
3496  /* 9-Jun-97 matchups - NOAA-9 before Jan 1994 */
3497  if ((l2rec->l1rec->solz[ip] >= solznight
3498  || (l2rec->l1rec->solz[ip] < solznight
3499  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3500  && dBt_37_12 < -0.4805) {
3501  if (dBt_37_12 < -1.249) {
3502  flags_sst[ip] |= SSTF_CLOUD;
3503  } else if (dBt_11_12 >= 0.239) {
3504  flags_sst[ip] |= SSTF_CLOUD;
3505  } else if (Bt37_maxmin[ip] >= 0.7975
3506  && dBt_11_12 >= -0.1805) {
3507  flags_sst[ip] |= SSTF_CLOUD;
3508  }
3509  } else {
3510  if (dBt_11_12 < 0.307
3511  && ((l2rec->l1rec->solz[ip] >= solznight
3512  || (l2rec->l1rec->solz[ip] < solznight
3513  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3514  && Bt37_maxmin[ip] >= 1.6215)) {
3515  flags_sst[ip] |= SSTF_CLOUD;
3516  }
3517  }
3518  } else {
3519  /* 9-Jun-97 matchups - NOAA-9 after Jan 1994 */
3520  if ((l2rec->l1rec->solz[ip] >= solznight
3521  || (l2rec->l1rec->solz[ip] < solznight
3522  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3523  && dBt_37_12 < 0.62) {
3524  flags_sst[ip] |= SSTF_CLOUD;
3525  } else if (dBt_11_12 < 0.5055
3526  && ((l2rec->l1rec->solz[ip] >= solznight
3527  || (l2rec->l1rec->solz[ip] < solznight
3528  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3529  && Bt37_maxmin[ip] >= 1.033)
3530  && Bt11_maxmin[ip] >= 0.2435) {
3531  flags_sst[ip] |= SSTF_CLOUD;
3532  }
3533  }
3534  break;
3535  case NO10:
3536  break;
3537  case NO11:
3538  if ((l2rec->l1rec->solz[ip] >= solznight
3539  || (l2rec->l1rec->solz[ip] < solznight && l2rec->l1rec->glint_coef[ip] <= glintmax))
3540  && dBt_37_12 < -0.581) {
3541  flags_sst[ip] |= SSTF_CLOUD;
3542  } else {
3543  if (dBt_11_12 < 0.315) {
3544  if ((l2rec->l1rec->solz[ip] >= solznight
3545  || (l2rec->l1rec->solz[ip] < solznight
3546  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3547  && dBt_11_37 < -0.1855) {
3548  flags_sst[ip] |= SSTF_CLOUD;
3549  } else if ((l2rec->l1rec->solz[ip] >= solznight
3550  || (l2rec->l1rec->solz[ip] < solznight
3551  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3552  && Bt37_maxmin[ip] >= 0.9555) {
3553  flags_sst[ip] |= SSTF_CLOUD;
3554  }
3555  } else if ((l2rec->l1rec->solz[ip] >= solznight
3556  || (l2rec->l1rec->solz[ip] < solznight
3557  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3558  && dBt_11_37 >= 0.86 && Bt37_maxmin[ip] >= 0.5035) {
3559  flags_sst[ip] |= SSTF_CLOUD;
3560  }
3561  }
3562  break;
3563  case NO12:
3564  break;
3565  case NO14:
3566  if (yyyyddd >= 1995001 && yyyyddd < 1996001) {
3567  if ((l2rec->l1rec->solz[ip] >= solznight
3568  || (l2rec->l1rec->solz[ip] < solznight
3569  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3570  && dBt_11_37 >= 0.8015) {
3571  flags_sst[ip] |= SSTF_CLOUD;
3572  } else if ((l2rec->l1rec->solz[ip] >= solznight
3573  || (l2rec->l1rec->solz[ip] < solznight
3574  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3575  && Bt37_maxmin[ip] >= 1.216) {
3576  if (dBt_11_12 < 1.342) {
3577  flags_sst[ip] |= SSTF_CLOUD;
3578  }
3579  } else if (dBt_11_12 < 1.3225
3580  && ((l2rec->l1rec->solz[ip] >= solznight
3581  || (l2rec->l1rec->solz[ip] < solznight
3582  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3583  && dBt_11_37 < 0.994)) {
3584  flags_sst[ip] |= SSTF_CLOUD;
3585  }
3586  } else if (yyyyddd >= 1996001) {
3587  if (dBt_11_12 >= 0.755) {
3588  if ((l2rec->l1rec->solz[ip] >= solznight
3589  || (l2rec->l1rec->solz[ip] < solznight
3590  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3591  && (dBt_11_37 >= 0.7995 || LtNIR8 >= 0.7225)) {
3592  flags_sst[ip] |= SSTF_CLOUD;
3593  }
3594  } else {
3595  if ((l2rec->l1rec->solz[ip] >= solznight
3596  || (l2rec->l1rec->solz[ip] < solznight
3597  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3598  && (dBt_11_37 >= 1.1475 || LtNIR8 >= 1.2645
3599  || Bt37_maxmin[ip] >= 0.966)) {
3600  flags_sst[ip] |= SSTF_CLOUD;
3601  } else {
3602  if (((l2rec->l1rec->solz[ip] >= solznight
3603  || (l2rec->l1rec->solz[ip] < solznight
3604  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3605  && dBt_11_37 >= 0.869)
3606  || (dBt_11_12 < 0.462
3607  && ((l2rec->l1rec->solz[ip] >= solznight
3608  || (l2rec->l1rec->solz[ip] < solznight
3609  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3610  && LtNIR8 >= 0.4865))) {
3611  flags_sst[ip] |= SSTF_CLOUD;
3612  }
3613  }
3614  }
3615  }
3616  break;
3617  case NO15:
3618  break;
3619  case NO16:
3620  /* pathfinder v6 trees */
3621  /* patfinder v6 2001-2003 */
3622  if (l2rec->l1rec->solz[ip] < solznight) {
3623  /* day */
3624  if ((l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7315)
3625  || dBt_11_12 < 0.5375) {
3626  flags_sst[ip] |= SSTF_CLOUD;
3627  }
3628  } else {
3629  /* night */
3630  if (dBt_37_11 < -0.1025 || Bt37_maxmin[ip] >= 1.035
3631  || dBt_11_12 < 0.7385) {
3632  flags_sst[ip] |= SSTF_CLOUD;
3633  } else if (Bt37_maxmin[ip] < 0.7145) {
3634  if (dBt_11_12 < 1.454) {
3635  flags_sst[ip] |= SSTF_CLOUD;
3636  }
3637  } else if (dBt_11_12 < 1.271) {
3638  if (dBt_37_11 >= 0.9195) {
3639  flags_sst[ip] |= SSTF_CLOUD;
3640  }
3641  } else if (dBt_37_11 >= 2.408 && dBt_11_12 < 2.272) {
3642  flags_sst[ip] |= SSTF_CLOUD;
3643  }
3644  }
3645  break;
3646  case NO17:
3647  /* pathfinder v6 trees */
3648  /* NO17 matchups 2003-2006 - pathfinder v6 */
3649  if (l2rec->l1rec->solz[ip] < solznight) {
3650  /* day */
3651  if (l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7045) {
3652  flags_sst[ip] |= SSTF_CLOUD;
3653  } else if (l2rec->l1rec->lat[ip] >= -20.0 && dBt_11_12 < 0.4355) {
3654  flags_sst[ip] |= SSTF_CLOUD;
3655  }
3656  } else {
3657  /* night */
3658  if (dBt_37_12 < 0.0355 || Bt37_maxmin[ip] >= 0.9035
3659  || dBt_11_12 < 0.4565 || dBt_37_11 < -0.5515) {
3660  flags_sst[ip] |= SSTF_CLOUD;
3661  } else if (dBt_37_11 >= 1.635 && dBt_11_12 < 1.693) {
3662  flags_sst[ip] |= SSTF_CLOUD;
3663  }
3664  }
3665  break;
3666  case NO18:
3667  /* pathfinder v6 trees */
3668  if (l2rec->l1rec->solz[ip] < solznight) {
3669  /* day */
3670  /* NOA18 matchups 2005-2009 tree - Pathfinder version 6 */
3671  if (l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7045) {
3672  flags_sst[ip] |= SSTF_CLOUD;
3673  } else if (dBt_11_12 < 0.4995) {
3674  flags_sst[ip] |= SSTF_CLOUD;
3675  }
3676  } else {
3677  /* night */
3678  /* NOA18 matchups 2005-2009 tree - Pathfinder version 6 */
3679  if (dBt_37_12 < 0.3185) {
3680  flags_sst[ip] |= SSTF_CLOUD;
3681  } else if (Bt37_maxmin[ip] >= 0.8455) {
3682  if (dBt_11_12 < 1.477) {
3683  flags_sst[ip] |= SSTF_CLOUD;
3684  }
3685  } else if (dBt_11_12 < 0.6605) {
3686  flags_sst[ip] |= SSTF_CLOUD;
3687  } else if (dBt_37_11 < -0.5745) {
3688  flags_sst[ip] |= SSTF_CLOUD;
3689  } else if (dBt_37_11 >= 2.1015 && dBt_11_12 < 2.206) {
3690  flags_sst[ip] |= SSTF_CLOUD;
3691  }
3692  }
3693  break;
3694  case NO19:
3695  if (l2rec->l1rec->solz[ip] < solznight) {
3696  /* day */
3697  if ((l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.84648)
3698  || dBt_11_12 < 0.168655
3699  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 1.237905) {
3700  flags_sst[ip] |= SSTF_CLOUD;
3701  }
3702  } else {
3703  /* night */
3704  if (dBt_11_12 < 0.99571) {
3705  if (dBt_37_12 < -0.26932 || dBt_11_12 < 0.126555
3706  || Bt37_maxmin[ip] >= 1.02994
3707  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 0.758065) {
3708  flags_sst[ip] |= SSTF_CLOUD;
3709  }
3710  } else if (dBt_37_11 < -0.387815
3711  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 0.857825) {
3712  flags_sst[ip] |= SSTF_CLOUD;
3713  }
3714  }
3715  break;
3716  default:
3717  /* modis trees are only v6 */
3718 
3719  /* modis trees were built using gsfc extractions with rho for the vis bands */
3720  /* but l2gen puts radiance in LtRED so need to use rhotRED in the tree tests */
3721 
3722  if (l2rec->l1rec->l1file->sensorID == MODIST) {
3723  if (l2rec->l1rec->solz[ip] < solznight
3724  && l2rec->l1rec->glint_coef[ip] <= glintmax) {
3725  /* day not glint TERRA SST tree test */
3726  treesum[ip] = 0.0;
3727  /* seadas (594,662) is outside glint for comparison to glint pixels that are wrong quality */
3728  /* Bad Lt's could be -32767, 0 (class only so ignore), or 1000. */
3729  /* rho is the same sign so really could check rhotRED instead of LtRED,
3730  except where we have to check for 1000.0 */
3731  if (rhotRED[ip] < 0.206 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3732  treesum[ip] += 0.894;
3733  if (rhotRED[ip] < 0.051 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3734  treesum[ip] += 0.540;
3735  if (l2rec->l1rec->senz[ip] < 64.915) {
3736  treesum[ip] += -0.057;
3737  if (rhoCirrus_max[ip] < 0.003 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3738  treesum[ip] += 0.087;
3739  if (l2rec->l1rec->senz[ip] < 54.529) {
3740  treesum[ip] += 0.030;
3741  } else {
3742  treesum[ip] += -0.573;
3743  }
3744  } else {
3745  treesum[ip] += -0.828;
3746  }
3747  } else {
3748  treesum[ip] += -2.700;
3749  }
3750  } else {
3751  treesum[ip] += -0.490;
3752  }
3753  /* default min is 32767, but gets set to -32767 if no valid values in box */
3754  if (rhoCirrus_min[ip] < 0.004 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3755  treesum[ip] += 0.210;
3756  if (rhotRED[ip] < 0.073 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3757  treesum[ip] += 0.074;
3758  } else {
3759  treesum[ip] += -0.467;
3760  }
3761  } else {
3762  treesum[ip] += -1.008;
3763  }
3764  } else {
3765  treesum[ip] += -1.402;
3766  /* don't need LtRED_min: if Lt is -32767, 0.0 or
3767  1000.0 (which makes rho large) it will go to
3768  the negative treesum */
3769  /* default min is 32767, but gets set to -32767 if no valid values in box */
3770  if (rhotRED_min[ip] < 0.239 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
3771  treesum[ip] += 0.821;
3772  } else {
3773  treesum[ip] += -0.033;
3774  if (dBt_67_11 < -42.423) {
3775  treesum[ip] += 0.283;
3776  } else {
3777  treesum[ip] += -0.184;
3778  }
3779  }
3780  }
3781  if (dBt_37_11 < 3.493) {
3782  treesum[ip] += 0.530;
3783  if (dBt_11_12 < 0.178) {
3784  treesum[ip] += -0.589;
3785  } else {
3786  treesum[ip] += 0.174;
3787  }
3788  } else {
3789  treesum[ip] += -0.169;
3790  /* default min is 32767, but gets set to -32767 if no valid values in box */
3791  if (rhoCirrus_min[ip] < 0.008 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3792  treesum[ip] += 0.097;
3793  } else {
3794  treesum[ip] += -0.761;
3795  }
3796  }
3797  if (sst_stdev[ip] < 0.273 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3798  treesum[ip] += 0.258;
3799  } else {
3800  treesum[ip] += -0.194;
3801  }
3802  if (sst_stdev[ip] < 0.119 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3803  treesum[ip] += 0.459;
3804  } else {
3805  treesum[ip] += -0.037;
3806  if (sstq[csstbox][ip] < 27.293) {
3807  treesum[ip] += -0.054;
3808  } else {
3809  treesum[ip] += 0.411;
3810  }
3811  }
3812  } else if (l2rec->l1rec->solz[ip] < solznight) {
3813  /* day in high and low glint TERRA SST tree tests */
3814  if (LtRED < (1000.0 - 0.1)) {
3815  /* RED is not saturated */
3816  /* day in low glint TERRA SST tree test */
3817  treesum[ip] = 0.0;
3818  /* seadas (600,663) should be qual 1 not 3; terra 20141701425 */
3819  /* seadas (611,698) should be qual 3 not 1; descending so not flipped */
3820  /* don't check for > 0.0 here because we want the negative
3821  treesum if it's < 0.0 or = 1000.0 */
3822  //if (rhotRED[ip] < 0.098 || LtRED >= (1000.0-0.1)) {}
3823  /* Jun 2017, don't know why the above comment and test, but want positive tree sum only if Lt was valid */
3824  if (rhotRED[ip] < 0.098 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3825  treesum[ip] += 1.026;
3826  if (rhotRED[ip] < 0.069 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3827  treesum[ip] += 0.290;
3828  } else {
3829  treesum[ip] += -0.312;
3830  }
3831  /* rhoCirrus is -32767 when bad */
3832  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
3833  treesum[ip] += 0.154;
3834  } else {
3835  treesum[ip] += -0.574;
3836  }
3837  /* if bad then it will go negative */
3838  if (Bt73_max[ip] < -10.303) {
3839  treesum[ip] += -0.150;
3840  } else {
3841  treesum[ip] += 0.393;
3842  }
3843  } else {
3844  treesum[ip] += -0.761;
3845  if (Bt11 < 16.122) {
3846  treesum[ip] += -0.538;
3847  } else {
3848  treesum[ip] += 0.707;
3849  }
3850  }
3851  /* rhoCirrus is -32767 when bad so make sure it goes negative if bad */
3852  if (rhoCirrus_max[ip] < 0.005 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3853  treesum[ip] += 0.357;
3854  if (Tdeflong < 0.037) {
3855  treesum[ip] += -0.450;
3856  if (l2rec->l1rec->lat[ip] < 28.270) {
3857  treesum[ip] += -0.383;
3858  } else {
3859  treesum[ip] += 0.560;
3860  }
3861  } else {
3862  treesum[ip] += 0.134;
3863  }
3864  if (l2rec->l1rec->senz[ip] < 34.756) {
3865  treesum[ip] += 0.098;
3866  } else {
3867  treesum[ip] += -0.424;
3868  }
3869  } else {
3870  treesum[ip] += -1.399;
3871  /* default min is 32767, but gets set to -32767 if no valid values in box */
3872  if (rhoCirrus_min[ip] < 0.01 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3873  treesum[ip] += 0.326;
3874  } else {
3875  treesum[ip] += -1.404;
3876  }
3877  }
3878  /* sst_stdev is -32767 if bad so make sure it goes negative if bad */
3879  if (sst_stdev[ip] < 0.382 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3880  treesum[ip] += 0.130;
3881  if (l2rec->l1rec->lat[ip] < -5.005) {
3882  treesum[ip] += 0.400;
3883  } else {
3884  treesum[ip] += -0.107;
3885  }
3886  /* sst_stdev is -32767 if bad so make sure it goes negative if bad */
3887  if (sst_stdev[ip] < 0.224 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3888  treesum[ip] += 0.182;
3889  } else {
3890  treesum[ip] += -0.216;
3891  }
3892  } else {
3893  treesum[ip] += -0.632;
3894  }
3895  if (dBt_85_11 < -1.641) {
3896  treesum[ip] += 0.268;
3897  } else {
3898  treesum[ip] += -0.201;
3899  if (sstq[csstbox][ip] < 27.723) {
3900  treesum[ip] += -0.069;
3901  } else {
3902  treesum[ip] += 0.903;
3903  }
3904  }
3905  } else {
3906  /* RED is saturated */
3907  /* day high glint terra sst tree test */
3908  treesum[ip] = 0.0;
3909  /* default min is 32767, but gets set to -32767 if no valid values in box */
3910  if (Bt12_min[ip] < 15.376) {
3911  treesum[ip] += -0.908;
3912  if (dBt_11_12 < 0.27) {
3913  treesum[ip] += -0.446;
3914  } else {
3915  treesum[ip] += 0.568;
3916  }
3917  /* default min is 32767, but gets set to -32767 if no valid values in box */
3918  if (rhoCirrus_min[ip] < -0.001) {
3919  treesum[ip] += -0.981;
3920  } else {
3921  treesum[ip] += 0.137;
3922  }
3923  } else {
3924  treesum[ip] += 0.819;
3925  if (Tdeflong < 0.031) {
3926  treesum[ip] += -0.969;
3927  } else {
3928  treesum[ip] += 0.025;
3929  }
3930  /* bad is -32767 so it will go negative if bad */
3931  if (Bt73_max[ip] < -9.485) {
3932  treesum[ip] += -0.208;
3933  } else {
3934  treesum[ip] += 0.428;
3935  }
3936  }
3937  if (rhoCirrus_max[ip] < 0.005 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3938  treesum[ip] += 0.235;
3939  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
3940  treesum[ip] += 0.100;
3941  } else {
3942  treesum[ip] += -0.506;
3943  }
3944  } else {
3945  treesum[ip] += -1.176;
3946  }
3947  if (sst_stdev[ip] < 0.36 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3948  treesum[ip] += 0.222;
3949  } else {
3950  treesum[ip] += -0.657;
3951  }
3952  if (sstq[csstbox][ip] < 27.443) {
3953  treesum[ip] += -0.077;
3954  if (sst_stdev[ip] < 0.36 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3955  treesum[ip] += 0.355;
3956  } else {
3957  treesum[ip] += -0.179;
3958  }
3959  /* default min is 32767, but gets set to -32767 if no valid values in box */
3960  if (rhoCirrus_min[ip] < 0.01 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3961  treesum[ip] += 0.022;
3962  if (l2rec->l1rec->lat[ip] < 29.775) {
3963  treesum[ip] += -0.128;
3964  if (l2rec->l1rec->lat[ip] < -15.305) {
3965  treesum[ip] += 0.431;
3966  } else {
3967  treesum[ip] += -0.202;
3968  /* bad max is -32767 so it will go negative without a specific test */
3969  if (Bt11_max[ip] < 14.722) {
3970  treesum[ip] += -1.906;
3971  } else {
3972  treesum[ip] += 0.087;
3973  }
3974  }
3975  if (Tdeflong < 0.04) {
3976  treesum[ip] += -0.336;
3977  } else {
3978  treesum[ip] += 0.156;
3979  }
3980  } else {
3981  treesum[ip] += 0.379;
3982  }
3983  } else {
3984  treesum[ip] += -1.948;
3985  }
3986  } else {
3987  treesum[ip] += 0.688;
3988  }
3989  }
3990  } else {
3991  /* night TERRA SST tree test */
3992  treesum[ip] = 0.0;
3993  if (dBt_37_12 < -0.053) {
3994  treesum[ip] += -1.257;
3995  if (l2rec->l1rec->lat[ip] < 33.195) {
3996  treesum[ip] += -0.278;
3997  if (l2rec->l1rec->lat[ip] < -40.185) {
3998  treesum[ip] += 0.619;
3999  } else {
4000  treesum[ip] += -0.711;
4001  if (Bt37 < 6.477) {
4002  treesum[ip] += -3.733;
4003  } else {
4004  treesum[ip] += -0.111;
4005  }
4006  }
4007  } else {
4008  treesum[ip] += 0.333;
4009  }
4010  if (Bt37 < 9.372) {
4011  treesum[ip] += -0.292;
4012  } else {
4013  treesum[ip] += 0.764;
4014  }
4015  } else {
4016  treesum[ip] += 0.430;
4017  if (Bt11_maxmin[ip] < 0.486 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4018  treesum[ip] += 0.628;
4019  if (Bt40_stdev[ip] < 0.146 && Bt40_stdev[ip] > (BAD_FLT + 0.1)) {
4020  treesum[ip] += 0.177;
4021  } else {
4022  treesum[ip] += -0.723;
4023  }
4024  } else {
4025  treesum[ip] += -0.450;
4026  }
4027  if (dSST_SST4 < -0.878) {
4028  treesum[ip] += -1.353;
4029  if (dSST_SST4 < -1.533) {
4030  treesum[ip] += -1.439;
4031  } else {
4032  treesum[ip] += 0.346;
4033  }
4034  } else {
4035  treesum[ip] += 0.219;
4036  if (Bt37_stdev[ip] < 0.448 && Bt37_stdev[ip] > (BAD_FLT + 0.1)) {
4037  treesum[ip] += 0.290;
4038  if (dSST_SST4 < -0.422) {
4039  treesum[ip] += -0.504;
4040  } else {
4041  treesum[ip] += 0.268;
4042  }
4043  } else {
4044  treesum[ip] += -0.484;
4045  }
4046  if (Bt12 < 16.736) {
4047  treesum[ip] += -0.285;
4048  if (dBt_40_11 < -2.199) {
4049  treesum[ip] += 0.518;
4050  } else {
4051  treesum[ip] += -0.316;
4052  if (Bt12 < 11.896) {
4053  treesum[ip] += -0.527;
4054  } else {
4055  treesum[ip] += 0.400;
4056  }
4057  }
4058  } else {
4059  treesum[ip] += 0.500;
4060  }
4061  if (dSST_SST4 < 1.183) {
4062  treesum[ip] += 0.051;
4063  } else {
4064  treesum[ip] += -0.898;
4065  }
4066  }
4067  }
4068  }
4069  if (treesum[ip] <= 0.0) {
4070  flags_sst[ip] |= SSTF_CLOUD;
4071  }
4072  }/* end case for modis terra */
4073  else if (l2rec->l1rec->l1file->sensorID == MODISA) {
4074  if (l2rec->l1rec->solz[ip] < solznight && l2rec->l1rec->glint_coef[ip] <= glintmax) {
4075  /* day not glint AQUA SST tree test */
4076  treesum[ip] = 0.0;
4077  /* seadas (131,1112) should be qual 1 not 3; in l2gen it's (1353-131,2029-1112)=(1222,917) */
4078  /* seadas (1086,1375) should be qual 1 not 3; in l2gen it's (1353-1086,2029-1375)=(1264,654) */
4079  /* seadas says valid rhotRED is >= -0.1 & <= 1.0 */
4080  /* seadas says valid cirrus is >= 0.0 & <= 1.1 */
4081  if (rhotRED[ip] < 0.204 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4082  treesum[ip] += 0.982;
4083  if (rhotRED[ip] < 0.056 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4084  treesum[ip] += 0.328;
4085  /* bad cirrus is -32767 */
4086  if (rhoCirrus_max[ip] < 0.006 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4087  treesum[ip] += 0.151;
4088  if (l2rec->l1rec->senz[ip] < 59.043) {
4089  treesum[ip] += 0.029;
4090  } else {
4091  treesum[ip] += -0.618;
4092  }
4093  } else {
4094  treesum[ip] += -1.071;
4095  }
4096  } else {
4097  treesum[ip] += -0.565;
4098  }
4099  } else {
4100  treesum[ip] += -1.325;
4101  /* if Lt is 1000 then rhot will be large so this will go negative */
4102  /* default min is 32767, but gets set to -32767 if no valid values in box */
4103  if (rhotRED_min[ip] < 0.251 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4104  treesum[ip] += 0.81;
4105  } else {
4106  treesum[ip] += -0.021;
4107  if (l2rec->l1rec->lat[ip] < 48.685) {
4108  treesum[ip] += 0.605;
4109  } else {
4110  treesum[ip] += -0.088;
4111  if (l2rec->l1rec->lat[ip] < 37.755) {
4112  treesum[ip] += -0.311;
4113  } else {
4114  treesum[ip] += 0.189;
4115  }
4116  if (dBt_11_12 < 0.752) {
4117  treesum[ip] += -0.091;
4118  } else {
4119  treesum[ip] += 0.581;
4120  }
4121  }
4122  }
4123  // maybe add check to make sure Bt67 is valid? if ((Bt67 +/- 1000)
4124  if (dBt_67_11 < -38.023) {
4125  treesum[ip] += 0.221;
4126  } else {
4127  treesum[ip] += -0.201;
4128  }
4129  }
4130  /* default min is 32767, but gets set to -32767 if no valid values in box */
4131  if (rhoCirrus_min[ip] < 0.003 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4132  treesum[ip] += 0.22;
4133  if (rhoCirrus_max[ip] < 0.003 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4134  treesum[ip] += 0.069;
4135  } else {
4136  treesum[ip] += -0.422;
4137  }
4138  } else {
4139  treesum[ip] += -0.891;
4140  }
4141  if (dBt_37_11 < 3.577) {
4142  treesum[ip] += 0.719;
4143  } else {
4144  treesum[ip] += -0.227;
4145  if (sstq[csstbox][ip] < 27.668) {
4146  treesum[ip] += -0.017;
4147  } else {
4148  treesum[ip] += 0.472;
4149  }
4150  if (sst_stdev[ip] < 0.481 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4151  treesum[ip] += 0.094;
4152  } else {
4153  treesum[ip] += -0.322;
4154  }
4155  }
4156  if (sst_stdev[ip] < 0.263 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4157  treesum[ip] += 0.193;
4158  } else {
4159  treesum[ip] += -0.276;
4160  }
4161  } else if (l2rec->l1rec->solz[ip] < solznight) {
4162  /* day in high and low glint AQUA SST tree tests */
4163  if (LtRED < (1000.0 - 0.1)) {
4164  /* RED is not saturated */
4165  /* day in low glint AQUA SST tree test */
4166  /* v1 adtree */
4167  /* seadas (403,469) should be qual 1 not 3; in l2gen it's (1353-403,2029-469)=(950,1560) */
4168  treesum[ip] = 0.0;
4169  /* default min is 32767, but gets set to -32767 if no valid values in box */
4170  if (rhotRED_min[ip] < 0.087 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4171  treesum[ip] += 1.017;
4172  } else {
4173  treesum[ip] += -0.644;
4174  /* default min is 32767, but gets set to -32767 if no valid values in box */
4175  if (Bt11_min[ip] < 15.22) {
4176  treesum[ip] += -0.837;
4177  if (dBt_67_11 < -40.834) {
4178  treesum[ip] += 0.510;
4179  } else {
4180  treesum[ip] += -0.181;
4181  }
4182  } else {
4183  treesum[ip] += 0.870;
4184  if (Tdeflong < 0.036) {
4185  treesum[ip] += -0.859;
4186  } else {
4187  treesum[ip] += 0.215;
4188  }
4189  }
4190  if (l2rec->l1rec->lat[ip] < 29.675) {
4191  treesum[ip] += -0.099;
4192  } else {
4193  treesum[ip] += 0.467;
4194  }
4195  if (dBt_37_11 < 10.053) {
4196  treesum[ip] += -0.743;
4197  } else {
4198  treesum[ip] += 0.137;
4199  }
4200  }
4201  if (rhoCirrus_max[ip] < 0.004 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4202  treesum[ip] += 0.227;
4203  if (rhotRED[ip] < 0.068 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4204  treesum[ip] += 0.824;
4205  } else {
4206  treesum[ip] += -0.113;
4207  }
4208  if (rhoCirrus_maxmin[ip] < 0.001 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4209  treesum[ip] += 0.102;
4210  } else {
4211  treesum[ip] += -0.454;
4212  }
4213  } else {
4214  treesum[ip] += -1.096;
4215  if (rhoCirrus < 0.006 && rhoCirrus > (BAD_FLT + 0.1)) {
4216  treesum[ip] += 0.282;
4217  } else {
4218  treesum[ip] += -0.810;
4219  }
4220  }
4221  if (sst_stdev[ip] < 0.397 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4222  treesum[ip] += 0.235;
4223  if (l2rec->l1rec->lat[ip] < -0.365) {
4224  treesum[ip] += 0.240;
4225  } else {
4226  treesum[ip] += -0.194;
4227  if (l2rec->l1rec->lon[ip] < -67.480) {
4228  treesum[ip] += 0.241;
4229  } else {
4230  treesum[ip] += -0.229;
4231  }
4232  }
4233  if (sst_stdev[ip] < 0.188 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4234  treesum[ip] += 0.187;
4235  } else {
4236  treesum[ip] += -0.258;
4237  }
4238  } else {
4239  treesum[ip] += -0.848;
4240  }
4241  if (sstq[csstbox][ip] < 27.278) {
4242  treesum[ip] += -0.066;
4243  } else {
4244  treesum[ip] += 0.509;
4245  }
4246  // /* v2 adtree, might not be needed, try v1 first */
4247  // /* day in glint AQUA SST tree test */
4248  // treesum[ip] = 0.0;
4249  // if (rhotRED_min[ip] < 15.119 && rhotRED_min[ip] > 0.0) {
4250  // treesum[ip] += -0.759;
4251  // } else {
4252  // treesum[ip] += 0.794;
4253  // }
4254  // if (sst_stdev[ip] < 0.254) {
4255  // treesum[ip] += 0.304;
4256  // if (sst_stdev[ip] < 0.14) {
4257  // treesum[ip] += 0.237;
4258  // } else {
4259  // treesum[ip] += -0.276;
4260  // }
4261  // if (l2rec->l1rec->lat[ip] < -4.225){
4262  // treesum[ip] += 0.401;
4263  // } else {
4264  // treesum[ip] += -0.144;
4265  // }
4266  // } else {
4267  // treesum[ip] += -0.508;
4268  // }
4269  // if (dBt_11_12 < 0.303) {
4270  // treesum[ip] += -0.715;
4271  // } else {
4272  // treesum[ip] += 0.184;
4273  // if (l2rec->l1rec->lat[ip] < -23.905) {
4274  // treesum[ip] += 0.662;
4275  // } else {
4276  // treesum[ip] += -0.087;
4277  // }
4278  // if (sst_stdev[ip] < 0.494) {
4279  // treesum[ip] += 0.114;
4280  // } else {
4281  // treesum[ip] += -0.550;
4282  // }
4283  // }
4284  // if (l2rec->l1rec->lat[ip] < 32.415) {
4285  // treesum[ip] += -0.137;
4286  // if (sstq[csstbox][ip] < 26.408) {
4287  // treesum[ip] += -0.245;
4288  // } else {
4289  // treesum[ip] += 0.499;
4290  // }
4291  // if (Tdeflong < 0.038) {
4292  // treesum[ip] += -0.415;
4293  // } else {
4294  // treesum[ip] += 0.038;
4295  // }
4296  // } else {
4297  // treesum[ip] += 0.506;
4298  // }
4299  } else {
4300  /* RED is saturated */
4301  /* day high glint aqua sst tree test */
4302  treesum[ip] = 0.0;
4303  /* default min is 32767, but gets set to -32767 if no valid values in box */
4304  /* if min is bad this will go negative without a specific test */
4305  if (Bt12_min[ip] < 15.119) {
4306  treesum[ip] += -0.759;
4307  if (dBt_37_11 < 12.733) {
4308  treesum[ip] += 0.652;
4309  /* default min is 32767, but gets set to -32767 if no valid values in box */
4310  /* if min is bad this will go negative without a specific test */
4311  if (rhoCirrus_min[ip] < 0.001) {
4312  treesum[ip] += -1.096;
4313  } else {
4314  treesum[ip] += 0.130;
4315  }
4316  } else {
4317  treesum[ip] += -0.288;
4318  }
4319  if (Bt37_maxmin[ip] < 0.529 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4320  treesum[ip] += 0.926;
4321  } else {
4322  treesum[ip] += 0.026;
4323  }
4324  } else {
4325  treesum[ip] += 0.794;
4326  if (dBt_37_11 < 8.37) {
4327  treesum[ip] += 0.522;
4328  } else {
4329  treesum[ip] += -0.268;
4330  }
4331  }
4332  if (rhoCirrus_maxmin[ip] < 0.001 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4333  treesum[ip] += 0.381;
4334  if (rhoCirrus_max[ip] < 0.002 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4335  treesum[ip] += 0.006;
4336  } else {
4337  treesum[ip] += -0.584;
4338  }
4339  } else {
4340  treesum[ip] += -0.619;
4341  }
4342  if (sst_stdev[ip] < 0.397 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4343  treesum[ip] += 0.101;
4344  } else {
4345  treesum[ip] += -0.684;
4346  }
4347  if (Tdeflong < 0.036) {
4348  treesum[ip] += -0.471;
4349  if (l2rec->l1rec->lat[ip] < 32.335) {
4350  treesum[ip] += -0.481;
4351  } else {
4352  treesum[ip] += 0.684;
4353  }
4354  } else {
4355  treesum[ip] += 0.193;
4356  }
4357  /* default min is 32767, but gets set to -32767 if no valid values in box */
4358  if (rhoCirrus_min[ip] < 0.003 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4359  treesum[ip] += 0.108;
4360  if (l2rec->l1rec->lat[ip] < -25.425) {
4361  treesum[ip] += 0.517;
4362  } else {
4363  treesum[ip] += -0.082;
4364  if (l2rec->l1rec->lat[ip] < 27.215) {
4365  treesum[ip] += -0.077;
4366  if (Bt11 < 14.734) {
4367  treesum[ip] += -2.012;
4368  } else {
4369  treesum[ip] += 0.089;
4370  }
4371  } else {
4372  treesum[ip] += 0.383;
4373  }
4374  }
4375  } else {
4376  treesum[ip] += -0.913;
4377  }
4378  if (sstq[csstbox][ip] < 27.583) {
4379  treesum[ip] += -0.074;
4380  } else {
4381  treesum[ip] += 0.691;
4382  }
4383  }
4384  } else {
4385  /* night AQUA SST tree test */
4386  /* use tree that starts with 0, not -0.029 */
4387  treesum[ip] = 0.0;
4388  if (dBt_37_12 < 0.117) {
4389  treesum[ip] += -1.279;
4390  if (l2rec->l1rec->lat[ip] < 33.035) {
4391  treesum[ip] += -0.289;
4392  if (l2rec->l1rec->lat[ip] < -42.235) {
4393  treesum[ip] += 0.747;
4394  } else {
4395  treesum[ip] += -0.564;
4396  if (Bt37 < 13.117) {
4397  treesum[ip] += -0.580;
4398  } else {
4399  treesum[ip] += 0.638;
4400  }
4401  }
4402  } else {
4403  treesum[ip] += 0.355;
4404  }
4405  if (Bt37 < 9.447) {
4406  treesum[ip] += -0.307;
4407  } else {
4408  treesum[ip] += 0.747;
4409  }
4410  } else {
4411  treesum[ip] += 0.470;
4412  if (Bt11_stdev[ip] < 0.155 && Bt11_stdev[ip] > (BAD_FLT + 0.1)) {
4413  treesum[ip] += 0.690;
4414  if (Bt37_maxmin[ip] < 0.524 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4415  treesum[ip] += 0.150;
4416  } else {
4417  treesum[ip] += -0.794;
4418  }
4419  } else {
4420  treesum[ip] += -0.430;
4421  }
4422  if (dSST_SST4 < -0.787) {
4423  treesum[ip] += -1.404;
4424  if (dSST_SST4 < -1.253) {
4425  treesum[ip] += -1.086;
4426  } else {
4427  treesum[ip] += 0.388;
4428  }
4429  } else {
4430  treesum[ip] += 0.197;
4431  if (Bt37_maxmin[ip] < 1.229 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4432  treesum[ip] += 0.287;
4433  if (dSST_SST4 < -0.383) {
4434  treesum[ip] += -0.531;
4435  } else {
4436  treesum[ip] += 0.279;
4437  }
4438  } else {
4439  treesum[ip] += -0.497;
4440  }
4441  if (Bt12 < 16.353) {
4442  treesum[ip] += -0.300;
4443  if (dBt_40_11 < -2.171) {
4444  treesum[ip] += 0.516;
4445  } else {
4446  treesum[ip] += -0.395;
4447  if (sst4q[csstbox][ip] < 16.212) {
4448  treesum[ip] += -0.694;
4449  } else {
4450  treesum[ip] += 0.392;
4451  }
4452  }
4453  } else {
4454  treesum[ip] += 0.451;
4455  }
4456  }
4457  // maybe add check to make sure Bt85 is valid? if ((Bt85 +/- 1000)
4458  if (dBt_85_11 < -1.577) {
4459  treesum[ip] += 0.088;
4460  } else {
4461  treesum[ip] += -0.408;
4462  }
4463  }
4464  }
4465  if (treesum[ip] <= 0.0) {
4466  flags_sst[ip] |= SSTF_CLOUD;
4467  }
4468  }/* end case for modis */
4469 
4470  else if (l2rec->l1rec->l1file->sensorID == VIIRSN) {
4471  /* if we run this more than once I'll add a switch, but I think it's just here to make a 'previous' picture */
4472  v5viirs = 0;
4473  /* v5 viirs */
4474  if (v5viirs == 1) {
4475  if (l2rec->l1rec->solz[ip] < solznight) {
4476  /* day VIIRS SST tree test */
4477  /* v5.33 day tree */
4478 
4479  /* sd6.4 version used class files that contained values in K */
4480  /* sd7.4 version uses gsfc files with values in C */
4481  if (dBt_11_12 < 0.23435 ||
4482  Bt12 < (270.78 - CtoK)) {
4483  flags_sst[ip] |= SSTF_CLOUD;
4484  } else if (Bt37_maxmin[ip] >= 0.40919 &&
4485  (Bt37_maxmin[ip] >= 1.0765 ||
4486  Bt12_maxmin[ip] >= 0.8412 ||
4487  Bt12_maxmin[ip] < 0.3193)) {
4488  flags_sst[ip] |= SSTF_CLOUD;
4489  }
4490  } else {
4491  /* night VIIRS SST tree test */
4492  /* v5.33 night tree */
4493  if (dBt_37_11 < 0.18654) {
4494  flags_sst[ip] |= SSTF_CLOUD;
4495  } else if (dSST_SST3 < -1.021674) {
4496  if (dBt_11_12 < 2.187215 || dBt_37_11 >= 5.46243) {
4497  flags_sst[ip] |= SSTF_CLOUD;
4498  }
4499  } else if (Bt37_maxmin[ip] >= 0.48349) {
4500  if (dBt_11_12 < 0.32251 || dSST_SST3 >= 0.4414579) {
4501  flags_sst[ip] |= SSTF_CLOUD;
4502  } else if (Bt37_maxmin[ip] >= 1.012265) {
4503  if (dBt_11_12 < 0.951295 || Bt12_maxmin[ip] < 0.69827) {
4504  flags_sst[ip] |= SSTF_CLOUD;
4505  }
4506  } else {
4507  if (l2rec->l1rec->senz[ip] > 56.3) {
4508  if (dBt_37_11 < 1.101845 || dBt_11_12 >= 1.308475) {
4509  flags_sst[ip] |= SSTF_CLOUD;
4510  }
4511  }
4512  }
4513  }
4514  }
4515  } else {
4516 
4517  /* v6 viirs tree */
4518  if (l2rec->l1rec->solz[ip] >= solznight) {
4519  /* night VIIRS v6 SST tree test */
4520  treesum[ip] = 0.413;
4521  if (dBt_37_11 < 0.115) {
4522  treesum[ip] += -1.931;
4523  } else {
4524  treesum[ip] += 0.496;
4525  if (dSST_SST3 < -0.465) {
4526  treesum[ip] += -0.942;
4527  if (dSST_SST3 < -1.102) {
4528  treesum[ip] += -0.735;
4529  if (dBt_37_11 < 1.686) {
4530  treesum[ip] += 0.651;
4531  } else {
4532  treesum[ip] += -0.82;
4533  }
4534  } else {
4535  treesum[ip] += 0.207;
4536  }
4537  } else {
4538  treesum[ip] += 0.309;
4539  if (dSST_SST3 < 0.78) {
4540  treesum[ip] += 0.043;
4541  if (dSST_SST3 < -0.158) {
4542  treesum[ip] += -0.392;
4543  } else {
4544  treesum[ip] += 0.194;
4545  }
4546  } else {
4547  treesum[ip] += -1.069;
4548  }
4549  if (dBt_37_11 < 0.306) {
4550  treesum[ip] += -0.692;
4551  } else {
4552  treesum[ip] += 0.077;
4553  }
4554  }
4555  if (Bt12 < 287.407) {
4556  treesum[ip] += -0.339;
4557  } else {
4558  treesum[ip] += 0.426;
4559  }
4560  }
4561  if (sstq[csstbox][ip] < (270.414 - CtoK)) {
4562  treesum[ip] += -3.879;
4563  } else {
4564  treesum[ip] += 0.047;
4565  if (Bt37_stdev[ip] < 0.247 && Bt37_stdev[ip] > (BAD_FLT + 0.1)) {
4566  treesum[ip] += 0.176;
4567  } else {
4568  treesum[ip] += -0.187;
4569  }
4570  }
4571  if (treesum[ip] <= 0.0) {
4572  /* set bit to show failed tree test */
4573  flags_sst[ip] |= SSTF_CLOUD;
4574  }
4575  /* end case for v6 viirs night */
4576  } else if (l2rec->l1rec->glint_coef[ip] <= glintmax) {
4577  /* day not glint VIIRS v6 SST tree test */
4578  treesum[ip] = 0.0;
4579 
4580  /* add ice test here for now. May want to add it to set_qual instead of setting cloud flag */
4581  /* only for viirs, day, not glint */
4582  // check for bad rhotred? if (rhotRED[ip] > 0.5 && LtRED > (BAD_FLT+0.1) && LtRED < (1000.0-0.1) &&
4583 
4584  /* instead of lat > 40, calc sub solar point and +/- 30 degrees to get lat bounds that vary by season */
4585  /* units of x=radians=(deg/day)*day/(deg/rad)
4586  * subsolar is negative in the southern hemisphere
4587  */
4588  xdoy = day + 284.0;
4589  xrad = (360.0 / 365.0) * xdoy / RADEG;
4590  subsolar = 23.45 * sin(xrad);
4591 
4592  if (rhotRED[ip] > 0.3 &&
4593  rhot16[ip] >= 0.006 && rhot16[ip] < 0.1 &&
4594  ((l2rec->l1rec->lat[ip] > (subsolar + 30.0)) || (l2rec->l1rec->lat[ip] < (subsolar - 30.0)))) {
4595  /* set cloud bit to show ice (or wispy clouds?) */
4596  flags_sst[ip] |= SSTF_CLOUD;
4597  /* don't really need to do the tree test now, but it doesn't hurt */
4598  }
4599 
4600  /* rho's are bad when Lt's are -32767 or 0.0 (make rho <=0) or 1000 (make rho large) */
4601  // don't adjust for different f0 in Kays R code was 1.0 (from sd6.4) should have been 25.084 (what it is in sd7)
4602  /* seadas (1907,3148) should be qual 0 not 3; in l2gen it's (3199-1907,3231-3148)=(1292,83) */
4603  /* seadas (1935,3147) should be qual 0 not 0; in l2gen it's (3199-1935,3231-3147)=(1264,84) */
4604  if (rhot16[ip] < 0.16 && Lt16 > (BAD_FLT + 0.1) && Lt16 < (1000.0 - 0.1)) {
4605  treesum[ip] += 0.805;
4606  if (rhotNIR7[ip] < 0.062 && LtNIR7 > (BAD_FLT + 0.1) && LtNIR7 < (1000.0 - 0.1)) {
4607  treesum[ip] += 0.393;
4608  if (rhoCirrus < 0.004 && rhoCirrus > (BAD_FLT + 0.1)) {
4609  treesum[ip] += 0.287;
4610  if (Tdeflong < 0.002) {
4611  treesum[ip] += -0.681;
4612  } else {
4613  treesum[ip] += 0.026;
4614  if (rhotNIR7[ip] < 0.039 && LtNIR7 > (BAD_FLT + 0.1) && LtNIR7 < (1000.0 - 0.1)) {
4615  treesum[ip] += 0.364;
4616  } else {
4617  treesum[ip] += -0.21;
4618  }
4619  }
4620  } else {
4621  treesum[ip] += -1.244;
4622  }
4623  } else {
4624  treesum[ip] += -0.0572;
4625  /* default min is 32767, but gets set to -32767 if no valid values in box */
4626  if (rhot16_min[ip] < 0.032 && rhot16_min[ip] > (BAD_FLT + 0.1)) {
4627  treesum[ip] += 0.455;
4628  } else {
4629  treesum[ip] += -0.395;
4630  }
4631  }
4632  if (l2rec->l1rec->senz[ip] < 64.994) {
4633  treesum[ip] += 0.216;
4634  if (rhoCirrus < 0.007 && rhoCirrus > (BAD_FLT + 0.1)) {
4635  treesum[ip] += 0.065;
4636  } else {
4637  treesum[ip] += -1.077;
4638  }
4639  } else {
4640  treesum[ip] += -0.708;
4641  }
4642  } else {
4643  treesum[ip] += -1.755;
4644  if (rhot16[ip] < 0.266 && Lt16 > (BAD_FLT + 0.1) && Lt16 < (1000.0 - 0.1)) {
4645  treesum[ip] += 0.642;
4646  } else {
4647  treesum[ip] += -0.19;
4648  if (rhotRED_maxmin[ip] < 0.103 && rhotRED_maxmin[ip] > (BAD_FLT + 0.1)) {
4649  treesum[ip] += 0.425;
4650  } else {
4651  treesum[ip] += -0.195;
4652  }
4653  }
4654  if (dBt_11_12 < 0.235) {
4655  treesum[ip] += -0.189;
4656  } else {
4657  treesum[ip] += 0.411;
4658  }
4659  if (l2rec->l1rec->wv[ip] < 2.946) {
4660  treesum[ip] += 0.038;
4661  } else {
4662  treesum[ip] += -1.137;
4663  }
4664  }
4665  if (Bt11_maxmin[ip] < 0.762 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4666  treesum[ip] += 0.156;
4667  } else {
4668  treesum[ip] += -0.188;
4669  }
4670  if (l2rec->l1rec->wv[ip] < 1.315) {
4671  treesum[ip] += 0.327;
4672  } else {
4673  treesum[ip] += -0.054;
4674  if (sstq[csstbox][ip] < (278.171 - CtoK)) {
4675  treesum[ip] += -0.679;
4676  } else {
4677  treesum[ip] += 0.05;
4678  }
4679  }
4680  } else {
4681  /* day glint VIIRS v6 SST tree test */
4682  treesum[ip] = 0.0;
4683  /* seadas (1907,3148) should be qual 0 not 3; in l2gen it's (3199-1907,3231-3148)=(1292,83) */
4684  /* seadas (1935,3147) should be qual 0 not 0; in l2gen it's (3199-1935,3231-3147)=(1264,84) */
4685  if (rhotRED[ip] > 0.065) {
4686  /* high glint tree */
4687  treesum[ip] = 0.0;
4688  /* default min is 32767, but gets set to -32767 if no valid values in box */
4689  if (Bt85_min[ip] < (287.451 - CtoK)) {
4690  treesum[ip] += -0.812;
4691  if (l2rec->l1rec->lat[ip] < 32.315) {
4692  treesum[ip] += -0.296;
4693  if (Tdeflong < 0.001) {
4694  treesum[ip] += -1.109;
4695  if (l2rec->l1rec->lat[ip] < -30.0) {
4696  treesum[ip] += 0.525;
4697  } else {
4698  treesum[ip] += -1.827;
4699  }
4700  } else {
4701  treesum[ip] += 0.44;
4702  }
4703  if (l2rec->l1rec->wv[ip] < 2.065) {
4704  treesum[ip] += 0.49;
4705  } else {
4706  treesum[ip] += -1.104;
4707  if (Bt85 < (287.059 - CtoK)) {
4708  treesum[ip] += -1.071;
4709  } else {
4710  treesum[ip] += 0.727;
4711  }
4712  }
4713  } else {
4714  treesum[ip] += 0.669;
4715  }
4716  if (dBt_37_11 < 17.594) {
4717  treesum[ip] += 0.452;
4718  } else {
4719  treesum[ip] += -0.76;
4720  }
4721  } else {
4722  treesum[ip] += 0.858;
4723  if (l2rec->l1rec->lon[ip] < -69.475) {
4724  treesum[ip] += 0.512;
4725  } else {
4726  treesum[ip] += -0.176;
4727  if (l2rec->l1rec->lat[ip] < 1.01) {
4728  treesum[ip] += 0.64;
4729  } else {
4730  treesum[ip] += -0.176;
4731  if (tmonth < 6) {
4732  treesum[ip] += 0.69;
4733  } else {
4734  treesum[ip] += -0.305;
4735  }
4736  }
4737  }
4738  if (dBt_37_11 < 9.655) {
4739  treesum[ip] += 0.562;
4740  } else {
4741  treesum[ip] += -0.173;
4742  }
4743  if (Bt85 < (290.343 - CtoK)) {
4744  treesum[ip] += -0.39;
4745  } else {
4746  treesum[ip] += 0.243;
4747  if (Tdeflong < 0.003) {
4748  treesum[ip] += -0.817;
4749  } else {
4750  treesum[ip] += 0.267;
4751  }
4752  }
4753  if (l2rec->l1rec->lat[ip] < 22.465) {
4754  treesum[ip] += -0.206;
4755  } else {
4756  treesum[ip] += 0.523;
4757  }
4758  }
4759  if (Bt11_maxmin[ip] < 1.189 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4760  treesum[ip] += 0.059;
4761  } else {
4762  treesum[ip] += -1.22;
4763  }
4764  } else {
4765  /* low glint tree */
4766  treesum[ip] = 0.0;
4767  /* default min is 32767, but gets set to -32767 if no valid values in box */
4768  if (rhotNIR7_min[ip] < 0.104 && rhotNIR7_min[ip] > (BAD_FLT + 0.1)) {
4769  treesum[ip] += 0.91;
4770  if (rhotRED[ip] < 0.086 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4771  treesum[ip] += 0.518;
4772  /* default min is 32767, but gets set to -32767 if no valid values in box */
4773  if (rhotRED_min[ip] < 0.067 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4774  treesum[ip] += 0.558;
4775  } else {
4776  treesum[ip] += -0.263;
4777  /* don't check for > BAD_FLT here because we want the negative treesum if Lt < 0.0 (or =1000.0) */
4778  if (rhot16[ip] < 0.06 || Lt16 > (1000.0 - 0.1)) {
4779  treesum[ip] += -0.231;
4780  } else {
4781  treesum[ip] += 1.712;
4782  }
4783  /* don't check for > BAD_FLT here because we want the negative treesum if Lt < 0.0 (or =1000.0) */
4784  if (rhot16[ip] < 0.046 || Lt16 > (1000.0 - 0.1)) {
4785  treesum[ip] += -1.353;
4786  } else {
4787  treesum[ip] += 0.352;
4788  }
4789  }
4790  } else {
4791  treesum[ip] += -0.585;
4792  if (dBt_37_12 < 12.951) {
4793  treesum[ip] += -0.905;
4794  } else {
4795  treesum[ip] += 0.187;
4796  if (rhotRED[ip] < 0.098 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4797  treesum[ip] += 0.549;
4798  } else {
4799  treesum[ip] += -0.484;
4800  }
4801  }
4802  }
4803  } else {
4804  treesum[ip] += -1.819;
4805  /* default min is 32767, but gets set to -32767 if no valid values in box */
4806  if (rhotRED_min[ip] < 0.206 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4807  treesum[ip] += 0.467;
4808  } else {
4809  treesum[ip] += -1.18;
4810  if (rhotRED_maxmin[ip] < 0.037 && rhotRED_maxmin[ip] > (BAD_FLT + 0.1)) {
4811  treesum[ip] += 1.747;
4812  } else {
4813  treesum[ip] += -1.79;
4814  }
4815  }
4816  if (l2rec->l1rec->wv[ip] < 1.705) {
4817  treesum[ip] += 0.434;
4818  } else {
4819  treesum[ip] += -0.645;
4820  }
4821  }
4822  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
4823  treesum[ip] += 0.23;
4824  if (l2rec->l1rec->lat[ip] < 32.13) {
4825  treesum[ip] += -0.067;
4826  if (sstq[csstbox][ip] < (300.034 - CtoK)) {
4827  treesum[ip] += -0.146;
4828  } else {
4829  treesum[ip] += 0.824;
4830  }
4831  } else {
4832  treesum[ip] += 0.758;
4833  }
4834  } else {
4835  treesum[ip] += -1.153;
4836  /* default min is 32767, but gets set to -32767 if no valid values in box */
4837  if (rhoCirrus_min[ip] < 0.005 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4838  treesum[ip] += 0.28;
4839  } else {
4840  treesum[ip] += -0.939;
4841  }
4842  }
4843  if (l2rec->l1rec->senz[ip] < 33.204) {
4844  treesum[ip] += 0.2;
4845  } else {
4846  treesum[ip] += -0.331;
4847  }
4848  }
4849  }
4850  if (treesum[ip] <= 0.0) {
4851  flags_sst[ip] |= SSTF_CLOUD;
4852  }
4853 
4854  } /* end case for viirs v5 vs v6 */
4855 
4856  }/* end case for viirs */
4857  break;
4858  } /* end switch statement for tree tests */
4859 
4860  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
4861  /* Check stray sunlight test, M2B2 (shared bit with SST4DIFF) */
4862  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
4863 
4864  if (l2rec->l1rec->lat[ip] < 0.0
4865  && ((ASCEND == 1 && l2rec->l1rec->pixnum[ip] <= cpix)
4866  || (ASCEND == 0 && l2rec->l1rec->pixnum[ip] > cpix))
4867  && l2rec->l1rec->senz[ip] > 45.00) {
4868  flags_sst[ip] |= SSTF_SUNLIGHT;
4869  }
4870  }
4871 
4872  } /* End of pixel loop */
4873  return;
4874 }
4875 
4876 /* ----------------------------------------------------------------------------------- */
4877 /* set_flags_sst4() - set quality flags for short wave sea surface temperature */
4878 /* */
4879 /* B. Franz, SAIC, August 2005. */
4880 
4881 /* ----------------------------------------------------------------------------------- */
4882 void set_flags_sst4(l2str *l2rec) {
4883 
4884  extern l1qstr l1que;
4885  int32_t npix = l2rec->l1rec->npix;
4886  int32_t ip, ipbir;
4887  float Bt37, Bt39, Bt40, Bt85;
4888  float Bt11, Bt12;
4889  float dBt_37_12, dBt_40_11, dBt_85_11;
4890  float dBt_34;
4891  float dSST4_ref; /* sst4 - ref */
4892  float dSST_SST4; /* sst - sst4 */
4893  statstr statrec;
4894 
4895  /* check each pixel in scan */
4896  for (ip = 0; ip < npix; ip++) {
4897 
4898  d3940ref[ip] = sstbad;
4899 
4900  /* SST not processed */
4901  if (sstmasked(l2rec->l1rec->flags, ip)) {
4902  flags_sst4[ip] |= SSTF_ISMASKED;
4903  continue;
4904  }
4905 
4906  ipbir = ip * NBANDSIR;
4907 
4908  Bt37 = l2rec->l1rec->Bt[ipbir + ib37]; /* modis chan 20, avhrr cen3, viirs m12 */
4909  Bt39 = l2rec->l1rec->Bt[ipbir + ib39]; /* modis chan 22 */
4910  Bt40 = l2rec->l1rec->Bt[ipbir + ib40]; /* modis chan 23, viirs m13 */
4911  Bt85 = l2rec->l1rec->Bt[ipbir + ib85]; /* modis chan 29, viirs m14 */
4912  Bt11 = l2rec->l1rec->Bt[ipbir + ib11]; /* modis chan 31, avhrr cen4, viirs m15 */
4913  Bt12 = l2rec->l1rec->Bt[ipbir + ib12]; /* modis chan 32, avhrr cen5, viirs m16 */
4914 
4915  /* aqua detector 0 channel 23 has a problem so average detectors 9 and 1 instead */
4916  if ((l2rec->l1rec->l1file->sensorID == MODISA)
4917  && l2rec->l1rec->detnum == 0) {
4918  /* current row is center of box */
4919  int32_t is = l1que.nq / 2;
4920  if (btavg(ip, is, ib40, NBANDSIR, &statrec) > 0) {
4921  Bt40 = statrec.avg;
4922  }
4923  }
4924  dBt_37_12 = Bt37 - Bt12;
4925  dBt_40_11 = Bt40 - Bt11;
4926  dBt_85_11 = Bt85 - Bt11;
4927 
4928  /* if BT could not be computed (radiance out-of-range) */
4929  if (Bt39 < BT_LO + 0.1 || Bt39 > BT_HI - 0.1 || Bt40 < BT_LO + 0.1
4930  || Bt40 > BT_HI - 0.1) {
4931  flags_sst4[ip] |= SSTF_BTBAD;
4932  continue;
4933  }
4934 
4935  /* check BT range (don't care about Bt37 for sst4) */
4936  if (Bt39 < Btmin || Bt39 > Btmax || Bt40 < Btmin || Bt40 > Btmax40)
4937  flags_sst4[ip] |= SSTF_BTRANGE;
4938 
4939  /* check BT diff */
4940  dBt_34 = Bt39 - Bt40;
4941  if (dBt_34 < dBt4min || dBt_34 > dBt4max)
4942  flags_sst4[ip] |= SSTF_BTDIFF;
4943 
4944  /* check BT diff against BT reference */
4945 
4946  /* pathfinder v6 */
4947  d3940ref[ip] = btrefdiffv6(ip, Bt39, Bt40, l2rec);
4948 
4949  if (d3940ref[ip] < dBtrefmin || d3940ref[ip] > dBtrefmax)
4950  flags_sst4[ip] |= SSTF_BT4REFDIFF;
4951 
4952  /* check SST range */
4953  if (sst4q[csstbox][ip] < SSTmin || sst4q[csstbox][ip] > SSTmaxn)
4954  flags_sst4[ip] |= SSTF_SSTRANGE;
4955 
4956  /* check SST difference with reference */
4957  if (sst4q[csstbox][ip] < sstbad + 1.0 || l2rec->l1rec->sstref[ip] < sstbad + 1.0) {
4958  dSST4_ref = sstbad;
4959  } else {
4960  dSST4_ref = sst4q[csstbox][ip] - l2rec->l1rec->sstref[ip];
4961  }
4962  // the "coldonly" and equatorial aerosol tests are to be run by default, but if SSTMODS is set, don't
4963  if ((evalmask & SSTMODS) == 0) {
4964  /* evaluate change to cold-test only */
4965  /* set the flag bit if sst4 is too much colder than reference */
4966  if (dSST4_ref < -SSTdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
4967  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4968  if (dSST4_ref < -input->sstrefdif &&
4969  l2rec->l1rec->lat[ip] >= equatorialSouth && l2rec->l1rec->lat[ip] <= equatorialNorth &&
4970  l2rec->l1rec->lon[ip] >= equatorialWest && l2rec->l1rec->lon[ip] <= equatorialEast) {
4971  /* tighter test for avhrr between 10S and 30N and -105 to 105 longitude */
4972  /* equatorial aerosol test */
4973  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4974  }
4975  /* set the flag bit if sst4 is too much warmer than reference at night */
4976  if (l2rec->l1rec->solz[ip] >= solznight) {
4977  if (fabs(dSST4_ref) > SSTdiff)
4978  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4979  if (fabs(dSST4_ref) > SSTvdiff)
4980  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4981  }
4982  /* is sst4 way too much colder than reference? */
4983  if (dSST4_ref < -SSTvdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
4984  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4985  } else {
4986  if (l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
4987  if (fabs(dSST4_ref) > SSTdiff)
4988  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4989  } else {
4990  if (dSST4_ref < -SSTdiff || dSST4_ref > (SSTdiff + 1))
4991  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4992  }
4993  if (fabs(dSST4_ref) > SSTvdiff)
4994  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4995  }
4996  /* check SST4 difference with 11um SST */
4997  dSST_SST4 = sstq[csstbox][ip] - sst4q[csstbox][ip];
4998  if (sstq[csstbox][ip] > sstbad + 1.0) {
4999  if (fabs(dSST_SST4) < SST4diff1)
5000  flags_sst4[ip] |= SSTF_SST4DIFF;
5001  if (fabs(dSST_SST4) < SST4diff2)
5002  flags_sst4[ip] |= SSTF_SST4VDIFF;
5003  }
5004 
5005  /* check sensor zenith limits */
5006  if (l2rec->l1rec->senz[ip] > hisenz)
5007  flags_sst4[ip] |= SSTF_HISENZ;
5008  if (l2rec->l1rec->senz[ip] > vhisenz)
5009  flags_sst4[ip] |= SSTF_VHISENZ;
5010  // flag 2 edge pixels as SSTF_VHISENZ so quality gets set to 3
5011  if (l2rec->l1rec->pixnum[ip] < 2 || l2rec->l1rec->pixnum[ip] > (fullscanpix - 3))
5012  flags_sst4[ip] |= SSTF_VHISENZ;
5013  // set the last 4 pixels of the scan for Terra to VHISENZ
5014  // as there is an apparent calibration issue with the BT12 for those pixels
5015  if ((l2rec->l1rec->l1file->sensorID == MODIST) && (l2rec->l1rec->pixnum[ip] > 1349))
5016  flags_sst4[ip] |= SSTF_VHISENZ;
5017 
5018  // /* if 11um SST is cold, check for clouds (collect 5) */
5019  // /* sstcloud checks to make sure it's day and not glint */
5020  // if (sstq[csstbox][ip] > sstbad+1.0 && l2rec->l1rec->sstref[ip] > sstbad+1.0 &&
5021  // sstq[csstbox][ip]-l2rec->l1rec->sstref[ip] <= -1.0)
5022  // if (sstcloud(ip,cldbox,cldbox,cldthresh) == 1)
5023  // flags_sst4[ip] |= SSTF_REDNONUNIF;
5024 
5025  /* check homogeneity of BT */
5026 
5027  if (Bt39_maxmin[ip] > Bt39unif1)
5028  flags_sst4[ip] |= SSTF_BTNONUNIF;
5029  if (Bt39_maxmin[ip] > Bt39unif2)
5030  flags_sst4[ip] |= SSTF_BTVNONUNIF;
5031 
5032  if (Bt40_maxmin[ip] > Bt40unif1)
5033  flags_sst4[ip] |= SSTF_BTNONUNIF;
5034  if (Bt40_maxmin[ip] > Bt40unif2)
5035  flags_sst4[ip] |= SSTF_BTVNONUNIF;
5036  /* end of homogeneity checks */
5037 
5038  if (l2rec->l1rec->solz[ip] >= solznight) {
5039  if (l2rec->l1rec->l1file->sensorID == MODIST) {
5040  /* night TERRA SST tree test */
5041  treesum[ip] = 0.0;
5042  if (dBt_37_12 < -0.053) {
5043  treesum[ip] += -1.257;
5044  if (l2rec->l1rec->lat[ip] < 33.195) {
5045  treesum[ip] += -0.278;
5046  if (l2rec->l1rec->lat[ip] < -40.185) {
5047  treesum[ip] += 0.619;
5048  } else {
5049  treesum[ip] += -0.711;
5050  if (Bt37 < 6.477) {
5051  treesum[ip] += -3.733;
5052  } else {
5053  treesum[ip] += -0.111;
5054  }
5055  }
5056  } else {
5057  treesum[ip] += 0.333;
5058  }
5059  if (Bt37 < 9.372) {
5060  treesum[ip] += -0.292;
5061  } else {
5062  treesum[ip] += 0.764;
5063  }
5064  } else {
5065  treesum[ip] += 0.430;
5066  if (Bt11_maxmin[ip] < 0.486 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
5067  treesum[ip] += 0.628;
5068  if (Bt40_stdev[ip] < 0.146) {
5069  treesum[ip] += 0.177;
5070  } else {
5071  treesum[ip] += -0.723;
5072  }
5073  } else {
5074  treesum[ip] += -0.450;
5075  }
5076  if (dSST_SST4 < -0.878) {
5077  treesum[ip] += -1.353;
5078  if (dSST_SST4 < -1.533) {
5079  treesum[ip] += -1.439;
5080  } else {
5081  treesum[ip] += 0.346;
5082  }
5083  } else {
5084  treesum[ip] += 0.219;
5085  if (Bt37_stdev[ip] < 0.448) {
5086  treesum[ip] += 0.290;
5087  if (dSST_SST4 < -0.422) {
5088  treesum[ip] += -0.504;
5089  } else {
5090  treesum[ip] += 0.268;
5091  }
5092  } else {
5093  treesum[ip] += -0.484;
5094  }
5095  if (Bt12 < 16.736) {
5096  treesum[ip] += -0.285;
5097  if (dBt_40_11 < -2.199) {
5098  treesum[ip] += 0.518;
5099  } else {
5100  treesum[ip] += -0.316;
5101  if (Bt12 < 11.896) {
5102  treesum[ip] += -0.527;
5103  } else {
5104  treesum[ip] += 0.400;
5105  }
5106  }
5107  } else {
5108  treesum[ip] += 0.500;
5109  }
5110  if (dSST_SST4 < 1.183) {
5111  treesum[ip] += 0.051;
5112  } else {
5113  treesum[ip] += -0.898;
5114  }
5115  }
5116  }
5117 
5118  } else if (l2rec->l1rec->l1file->sensorID == MODISA) {
5119  /* MODIS AQUA SST4 Night decision Tree test */
5120  /* night AQUA SST tree test */
5121  /* use tree that starts with 0, not -0.029 */
5122  treesum[ip] += 0.0;
5123  if (dBt_37_12 < 0.117) {
5124  treesum[ip] += -1.279;
5125  if (l2rec->l1rec->lat[ip] < 33.035) {
5126  treesum[ip] += -0.289;
5127  if (l2rec->l1rec->lat[ip] < -42.235) {
5128  treesum[ip] += 0.747;
5129  } else {
5130  treesum[ip] += -0.564;
5131  if (Bt37 < 13.117) {
5132  treesum[ip] += -0.580;
5133  } else {
5134  treesum[ip] += 0.638;
5135  }
5136  }
5137  } else {
5138  treesum[ip] += 0.355;
5139  }
5140  if (Bt37 < 9.447) {
5141  treesum[ip] += -0.307;
5142  } else {
5143  treesum[ip] += 0.747;
5144  }
5145  } else {
5146  treesum[ip] += 0.470;
5147  if (Bt11_stdev[ip] < 0.155) {
5148  treesum[ip] += 0.690;
5149  if (Bt37_maxmin[ip] < 0.524 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
5150  treesum[ip] += 0.150;
5151  } else {
5152  treesum[ip] += -0.794;
5153  }
5154  } else {
5155  treesum[ip] += -0.430;
5156  }
5157  if (dSST_SST4 < -0.787) {
5158  treesum[ip] += -1.404;
5159  if (dSST_SST4 < -1.253) {
5160  treesum[ip] += -1.086;
5161  } else {
5162  treesum[ip] += 0.388;
5163  }
5164  } else {
5165  treesum[ip] += 0.197;
5166  if (Bt37_maxmin[ip] < 1.229 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
5167  treesum[ip] += 0.287;
5168  if (dSST_SST4 < -0.383) {
5169  treesum[ip] += -0.531;
5170  } else {
5171  treesum[ip] += 0.279;
5172  }
5173  } else {
5174  treesum[ip] += -0.497;
5175  }
5176  if (Bt12 < 16.353) {
5177  treesum[ip] += -0.300;
5178  if (dBt_40_11 < -2.171) {
5179  treesum[ip] += 0.516;
5180  } else {
5181  treesum[ip] += -0.395;
5182  if (sst4q[csstbox][ip] < 16.212) {
5183  treesum[ip] += -0.694;
5184  } else {
5185  treesum[ip] += 0.392;
5186  }
5187  }
5188  } else {
5189  treesum[ip] += 0.451;
5190  }
5191  }
5192  // maybe add check to make sure Bt85 is valid? if ((Bt85 +/- 1000)
5193  if (dBt_85_11 < -1.577) {
5194  treesum[ip] += 0.088;
5195  } else {
5196  treesum[ip] += -0.408;
5197  }
5198  }
5199  }
5200  if (treesum[ip] <= 0.0) {
5201  flags_sst4[ip] |= SSTF_CLOUD;
5202  }
5203  } /* end case for sst4 */
5204  } /* end pixel loop */
5205  return;
5206 }
5207 
5208 /* ----------------------------------------------------------------------------------- */
5209 /* set_qual_sst() - set quality levels for long wave sea surface temperature */
5210 /* */
5211 /* B. Franz, SAIC, August 2005. */
5212 
5213 /* ----------------------------------------------------------------------------------- */
5214 void set_qual_sst(l2str *l2rec) {
5215  int32_t ip, ipb;
5216  float rhot;
5217 
5218  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
5219 
5220  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
5221  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5222  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5223 
5224  qual_sst[ip] = 8;
5225 
5226  } else if ((flags_sst[ip] & SSTF_BTRANGE) == 0
5227  && (flags_sst[ip] & SSTF_BTVNONUNIF) == 0
5228  && (flags_sst[ip] & SSTF_VHISENZ) == 0
5229  && (flags_sst[ip] & SSTF_SSTRANGE) == 0
5230  && (flags_sst[ip] & SSTF_SUNLIGHT) == 0) {
5231 
5232  if ((flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5233  && (flags_sst[ip] & SSTF_BTNONUNIF) == 0
5234  && (flags_sst[ip] & SSTF_HISENZ) == 0
5235  && (flags_sst[ip] & SSTF_GLINT) == 0
5236  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5237  qual_sst[ip] = 0;
5238  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0
5239  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5240  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5241  qual_sst[ip] = 1;
5242  } else if ((flags_sst[ip] & SSTF_HISENZ) == 0
5243  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5244  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5245  qual_sst[ip] = 2;
5246  } else if ((flags_sst[ip] & SSTF_CLOUD) == 0
5247  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0) {
5248  qual_sst[ip] = 3;
5249  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0
5250  && (flags_sst[ip] & SSTF_HISENZ) == 0) {
5251  qual_sst[ip] = 4;
5252  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0) {
5253  qual_sst[ip] = 5;
5254  } else if ((flags_sst[ip] & SSTF_HISENZ) == 0) {
5255  qual_sst[ip] = 6;
5256  } else {
5257  qual_sst[ip] = 7;
5258  }
5259  if (l2rec->l1rec->ssttype[ip] == 0 && qual_sst[ip] < 4) {
5260  qual_sst[ip] = 4;
5261  }
5262  } else {
5263  qual_sst[ip] = 7;
5264  }
5265  } else {
5266 
5267  if (l2rec->l1rec->solz[ip] < solznight) {
5268 
5269  /* daytime 11um SST quality level */
5270 
5271  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5272  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5273 
5274  qual_sst[ip] = 4;
5275 
5276  } else if ((flags_sst[ip] & SSTF_VHISENZ) > 0
5277  || (flags_sst[ip] & SSTF_BTRANGE) > 0
5278  || (flags_sst[ip] & SSTF_SSTRANGE) > 0
5279  || (flags_sst[ip] & SSTF_BTVNONUNIF) > 0
5280  || (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0
5281  || (flags_sst[ip] & SSTF_CLOUD) > 0
5282  || l2rec->l1rec->ssttype[ip] == 0) {
5283 
5284  qual_sst[ip] = 3;
5285 
5286  } else if (((input->viirsnv7 >= 1)
5287  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5288  || (flags_sst[ip] & SSTF_SSTREFDIFF) > 0
5289  || (flags_sst[ip] & SSTF_REDNONUNIF) > 0) {
5290 
5291  qual_sst[ip] = 2;
5292 
5293  } else if (((input->viirsnv7 <= 0)
5294  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5295  || (l2rec->l1rec->glint_coef[ip] > glintmax)
5296  || (flags_sst[ip] & SSTF_HISENZ) > 0) {
5297  /* if this changes then change comp_sses_sstv6mv */
5298  /* because we use qual 0 sses stats for qual 1 due to glint only */
5299 
5300  qual_sst[ip] = 1;
5301 
5302  } else {
5303 
5304  qual_sst[ip] = 0;
5305 
5306  }
5307 
5308  // Kay thinks this was supposed to be taken out for modis also, but definitely don't want it for viirs?
5309  //if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
5310  /* Reduce quality if red-band reflectance is high and sst is cold (collect 5) */
5311  /* only in non glint areas */
5312  if (haveRed && l2rec-> l1rec->glint_coef[ip] <= glintmax) {
5313  ipb = ip * nbvis + ibred;
5314  /* Sue: should it divide by csolz here as it does everywhere else?:
5315  * rhot = PI * l2rec->l1rec->Lt[ipb] / l2rec->Fo[ibred] / l2rec->csolz[ip];
5316  * but then the 0.05 limit needs to change?
5317  */
5318  rhot = PI * l2rec->l1rec->Lt[ipb] / l2rec->l1rec->Fo[ibred];
5319  if (qual_sst[ip] < 3 && sstq[csstbox][ip] - l2rec->l1rec->sstref[ip] < -1.0
5320  && rhot > 0.05) {
5321  qual_sst[ip]++;
5322  if ((input->viirsnv7 >= 1)
5323  && (qual_sst[ip] == 1)) {
5324  /* viirsnv7 test want only hisenz to be qual 1 so degrade q0 to q2 */
5325  qual_sst[ip]++;
5326  }
5327 
5328  }
5329  }
5330  //}
5331  } else {
5332  /* nighttime 11um SST quality level */
5333 
5334  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5335  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5336 
5337  qual_sst[ip] = 4;
5338 
5339  } else if ((flags_sst[ip] & SSTF_BTRANGE) > 0
5340  || (flags_sst[ip] & SSTF_SSTRANGE) > 0
5341  || (flags_sst[ip] & SSTF_BT4REFDIFF) > 0
5342  || (flags_sst[ip] & SSTF_BTVNONUNIF) > 0
5343  || (flags_sst[ip] & SSTF_CLOUD) > 0
5344  || (flags_sst[ip] & SSTF_VHISENZ) > 0
5345  || (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0) {
5346  /* BTVNONUNIF should be qual 3 as it is for day */
5347  /* don't need to check for haveSST4 because BT4REFDIF is only set if haveSST4 */
5348 
5349  qual_sst[ip] = 3;
5350 
5351  } else if ((flags_sst[ip] & SSTF_SSTREFDIFF) > 0
5352  || (flags_sst[ip] & SSTF_SST4VDIFF) > 0
5353  || ((input->viirsnv7 >= 1)
5354  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5355  || ((input->viirsnv7 >= 1)
5356  && ((flags_sst[ip] & SSTF_SST4DIFF) > 0))
5357  || ((input->viirsnv7 >= 1)
5358  && ((flags_sst[ip] & SSTF_SST3DIFF) > 0))
5359  || (flags_sst[ip] & SSTF_SST3VDIFF) > 0) {
5360  /* for now sst4vdiff and sst3vdiff are the same bit */
5361  /* BTNONUNIF should be qual 2 as it is for day */
5362 
5363  qual_sst[ip] = 2;
5364 
5365  } else if (((input->viirsnv7 <= 0)
5366  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5367  || ((input->viirsnv7 <= 0)
5368  && ((flags_sst[ip] & SSTF_SST4DIFF) > 0))
5369  || ((input->viirsnv7 <= 0)
5370  && ((flags_sst[ip] & SSTF_SST3DIFF) > 0))
5371  || (flags_sst[ip] & SSTF_HISENZ) > 0) {
5372  /* for now sst4diff and sst3diff are the same bit */
5373 
5374  qual_sst[ip] = 1;
5375 
5376  } else {
5377 
5378  qual_sst[ip] = 0;
5379  }
5380 
5381  /* reduce quality if 4um BTs show non-uniformity (collect 5) */
5382  if (haveSST4 && sst4q[csstbox][ip] > sstbad + 1.0) {
5383  if (qual_sst[ip] < 3
5384  && (flags_sst4[ip] & SSTF_BTNONUNIF) > 0)
5385  qual_sst[ip]++;
5386  if ((input->viirsnv7 >= 1) && (qual_sst[ip] == 1)) {
5387  /* viirsnv7 test want only hisenz to be qual 1 so degrade q0 to q2 */
5388  qual_sst[ip]++;
5389  }
5390  }
5391 
5392  } /* end night sst quality */
5393 
5394  } /* end modis and viirs sst quality */
5395  } /* end pixel loop */
5396  return;
5397 }
5398 
5399 /* ----------------------------------------------------------------------------------- */
5400 /* set_qual_sst4() - set quality levels for short wave sea surface temperature */
5401 /* */
5402 /* B. Franz, SAIC, August 2005. */
5403 
5404 /* ----------------------------------------------------------------------------------- */
5405 void set_qual_sst4(l2str *l2rec) {
5406  int32_t ip;
5407 
5408  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
5409 
5410  if (l2rec->l1rec->solz[ip] < solznight) {
5411 
5412  /* daytime 4um SST quality level */
5413 
5414  if ((flags_sst4[ip] & SSTF_ISMASKED) > 0
5415  || (flags_sst4[ip] & SSTF_BTBAD) > 0) {
5416 
5417  qual_sst4[ip] = 4;
5418 
5419  } else {
5420 
5421  qual_sst4[ip] = 3;
5422 
5423  }
5424 
5425  } else {
5426 
5427  /* nighttime 4um SST quality level */
5428 
5429  if ((flags_sst4[ip] & SSTF_ISMASKED) > 0
5430  || (flags_sst4[ip] & SSTF_BTBAD) > 0) {
5431 
5432  qual_sst4[ip] = 4;
5433 
5434  } else if ((flags_sst4[ip] & SSTF_BTRANGE) > 0
5435  || (flags_sst4[ip] & SSTF_SSTRANGE) > 0
5436  || (flags_sst4[ip] & SSTF_BT4REFDIFF) > 0
5437  || (flags_sst4[ip] & SSTF_BTVNONUNIF) > 0
5438  || (flags_sst4[ip] & SSTF_CLOUD) > 0
5439  || (flags_sst4[ip] & SSTF_VHISENZ) > 0
5440  || (flags_sst4[ip] & SSTF_SSTREFVDIFF) > 0) {
5441  /* BTVNONUNIF and VHISENZ are very bad */
5442 
5443  qual_sst4[ip] = 3;
5444 
5445  } else if ((flags_sst4[ip] & SSTF_SSTREFDIFF) > 0
5446  || (flags_sst4[ip] & SSTF_SST4VDIFF) > 0) {
5447  /* I'm pretty sure this SSTF_SSTREFDIFF should be in flags_sst4, not flags_sst */
5448 
5449  qual_sst4[ip] = 2;
5450 
5451  } else if ((flags_sst4[ip] & SSTF_BTNONUNIF) > 0
5452  || (flags_sst4[ip] & SSTF_SST4DIFF) > 0
5453  || (flags_sst4[ip] & SSTF_HISENZ) > 0) {
5454 
5455  qual_sst4[ip] = 1;
5456 
5457  } else {
5458 
5459  qual_sst4[ip] = 0;
5460  }
5461  }
5462 
5463  }
5464  return;
5465 }
5466 
5467 /* ----------------------------------------------------------------------------------- */
5468 /* nlsst() - nonlinear sst, long wave sea surface temperature */
5469 /* */
5470 /* B. Franz, SAIC, August 2005. */
5471 
5472 /* ----------------------------------------------------------------------------------- */
5473 float nlsst(float Bt11, float Bt12, int32_t is, float sstref,
5474  float **bounds, float **coef, float sstrefoffday,
5475  float sstrefoffnight, int32_t ip, int32_t xsatid,
5476  size_t nlatbands, size_t ncoefficients) {
5477  extern l1qstr l1que;
5478  static float dBtlo = 0.5;
5479  static float dBthi = 0.9;
5480 
5481  float lat = l1que.r[is].lat[ip];
5482  float mu = l1que.r[is].csenz[ip];
5483 
5484  int32_t sensorID = l1que.r[is].l1file->sensorID;
5485  int mside = l1que.r[is].mside;
5486 
5487  float dBt = Bt11 - Bt12; /* aka "thing" */
5488  float sstlo, ssthi, lsst;
5489  float satzdir;
5490  int i, j;
5491 
5492  /* it doesn't hurt, but V5 and viirs osi-saf do not have latband coeffs so they don't use 'ic' */
5493  /* bounds go from -90 to 90 so we want the first set */
5494  /* that has max lat greater than pixel lat */
5495  int ic = nlatbands - 1;
5496  while ((ic > 0) && (lat <= bounds[ic - 1][1] + latwin)) {
5497  ic--;
5498  }
5499 
5500  /* convert sstref to K or whatever offset is required for the coeffs - usually 0.0 except for VIIRS */
5501  if (l1que.r[is].solz[ip] >= solznight)
5502  sstref = sstref + sstrefoffnight;
5503  else
5504  sstref = sstref + sstrefoffday;
5505 
5506  /* sstref is supposed to warm the retrieved temperature, not cool it */
5507  /* this is the scaling factor for temp deficit term */
5508  /* not a problem for deg F or K, but deg C can go negative and we don't want that */
5509  if (sstref < 0.0) {
5510  sstref = 0.0;
5511  }
5512  // VIIRS coefficients made for brightness temperatures are in Kelvin...so add the offset
5513  if (sensorID == VIIRSN) {
5514  Bt11 += CtoK
5515  ;
5516  }
5517  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
5518  satzdir = (l1que.r[is].pixnum[ip] < fullscanpix / 2) ? -1.0 : 1.0;
5519 
5520  if (isV5) {
5521  /* Noaa-16 uses different bounds than all the other AVHRR satellites */
5522  if (xsatid == NO16) {
5523  dBtlo = 1.0;
5524  dBthi = 1.4;
5525  }
5526  if(nlatbands < 2) {
5527  printf("ERROR - need nlatbands >= 2, found %ld\n", nlatbands);
5528  exit(EXIT_FAILURE);
5529  }
5530  if(ncoefficients < 4) {
5531  printf("ERROR - need ncoefficients >= 4, found %ld\n", ncoefficients);
5532  exit(EXIT_FAILURE);
5533  }
5534  if (dBt <= dBtlo)
5535  lsst = coef[0][0] + coef[0][1] * Bt11 + coef[0][2] * dBt * sstref + coef[0][3] * dBt * ((1.0 / mu) - 1.0);
5536  else if (dBt >= dBthi)
5537  lsst = coef[1][0] + coef[1][1] * Bt11 + coef[1][2] * dBt * sstref + coef[1][3] * dBt * ((1.0 / mu) - 1.0);
5538  else {
5539  sstlo = coef[0][0] + coef[0][1] * Bt11 + coef[0][2] * dBt * sstref + coef[0][3] * dBt * ((1.0 / mu) - 1.0);
5540  ssthi = coef[1][0] + coef[1][1] * Bt11 + coef[1][2] * dBt * sstref + coef[1][3] * dBt * ((1.0 / mu) - 1.0);
5541  lsst = sstlo + ((dBt - dBtlo) / (dBthi - dBtlo))*(ssthi - sstlo);
5542  }
5543 
5544  } else if (input->viirsnosisaf == 1) {
5545  if (l1que.r[is].solz[ip] >= solznight)
5546  ic = 2;
5547  else
5548  ic = 0;
5549  /* OSI-SAF only has one set of coeffs for day and one for night, not two */
5550  /* OSI-SAF equation: Ts = b0 + (b1 + bt2 Stheta) T11 + [b3 + b4 (Ts0 - 273.15) + b5 Stheta] dT11-12 + b6 Stheta */
5551  /* Ts0 is in degrees Kelvin. Our sstref is in Deg C so don't need to subract 273.15 */
5552  if(ncoefficients < 7) {
5553  printf("ERROR - need ncoefficients >= 7, found %ld\n", ncoefficients);
5554  exit(EXIT_FAILURE);
5555  }
5556  lsst = coef[ic][0]
5557  + (coef[ic][1] + coef[ic][2] * ((1.0 / mu) - 1.0)) * Bt11
5558  + (coef[ic][3] + coef[ic][4] * sstref + coef[ic][5] * ((1.0 / mu) - 1.0)) * dBt
5559  + coef[ic][6] * ((1.0 / mu) - 1.0);
5560 
5561  } else if (input->viirsnv7 >= 1) {
5562  /* v7 high senz latband coeffs */
5563  /* choose coeffs by latband */
5564 
5565  if (ncoefficients < 7) {
5566  printf("ERROR - need ncoefficients >= 7, found %ld\n", ncoefficients);
5567  exit(EXIT_FAILURE);
5568  }
5569  if (lat < bounds[ic][1] - latwin || ic == (nlatbands - 1)) {
5570  lsst = coef[ic][0] + coef[ic][1] * Bt11 + coef[ic][2] * dBt * sstref
5571  + coef[ic][3] * dBt * ((1.0 / mu) - 1.0)
5572  + coef[ic][4] * ((1.0 / mu) - 1.0)
5573  + coef[ic][5] * pow(((1.0 / mu) - 1.0), coef[ic][6]);
5574 
5575  } else {
5576  dBtlo = bounds[ic][1] - latwin;
5577  dBthi = bounds[ic][1] + latwin;
5578  sstlo = coef[ic][0] + coef[ic][1] * Bt11 + coef[ic][2] * dBt * sstref
5579  + coef[ic][3] * dBt * ((1.0 / mu) - 1.0)
5580  + coef[ic][4] * ((1.0 / mu) - 1.0)
5581  + coef[ic][5] * pow(((1.0 / mu) - 1.0), coef[ic][6]);
5582  ssthi = coef[ic + 1][0] + coef[ic + 1][1] * Bt11
5583  + coef[ic + 1][2] * dBt * sstref
5584  + coef[ic + 1][3] * dBt * ((1.0 / mu) - 1.0)
5585  + coef[ic + 1][4] * ((1.0 / mu) - 1.0)
5586  + coef[ic + 1][5] * pow(((1.0 / mu) - 1.0), coef[ic + 1][6]);
5587 
5588  lsst = sstlo + ((lat - dBtlo) / (dBthi - dBtlo)) * (ssthi - sstlo);
5589 
5590  }
5591 
5592  } else {
5593  /* v6 latband coeffs */
5594 
5595  /* Only coeff 4 for terra and aqua, they're zeroed for others */
5596  /* Only coeffs 5-6 for terra, aqua and viirs; they're zeroed for avhrr */
5597  static float* C0;
5598  static float* C1;
5599  static int last_ic = -1;
5600 
5601  if (C0 == NULL) {
5602  C0 = (float *) calloc(7, sizeof(float));
5603  C1 = (float *) calloc(7, sizeof(float));
5604  }
5605  if (last_ic != ic) {
5606  for (i = 0; i < ncoefficients; i++){
5607  j = i;
5608  if (((sensorID == VIIRSN) || (sensorID == AVHRR)) && (i > 3)){
5609  j += 1;
5610  }
5611  C0[j] = coef[ic][i];
5612  if (ic < nlatbands -1)
5613  C1[j] = coef[ic + 1][i];
5614  }
5615  }
5616  if (lat < bounds[ic][1] - latwin || ic == (nlatbands - 1)) {
5617  lsst = C0[0]
5618  + C0[1] * Bt11
5619  + C0[2] * dBt * sstref
5620  + C0[3] * dBt * ((1.0 / mu) - 1.0)
5621  + C0[4] * mside
5622  + C0[5] * (l1que.r[is].senz[ip] * satzdir)
5623  + C0[6] * pow((l1que.r[is].senz[ip] * satzdir), 2);
5624  } else {
5625  dBtlo = bounds[ic][1] - latwin;
5626  dBthi = bounds[ic][1] + latwin;
5627  sstlo = C0[0]
5628  + C0[1] * Bt11
5629  + C0[2] * dBt * sstref
5630  + C0[3] * dBt * ((1.0 / mu) - 1.0)
5631  + C0[4] * mside
5632  + C0[5] * (l1que.r[is].senz[ip] * satzdir)
5633  + C0[6] * pow((l1que.r[is].senz[ip] * satzdir), 2);
5634  ssthi = C1[0]
5635  + C1[1] * Bt11
5636  + C1[2] * dBt * sstref
5637  + C1[3] * dBt * ((1.0 / mu) - 1.0)
5638  + C1[4] * mside
5639  + C1[5] * (l1que.r[is].senz[ip] * satzdir)
5640  + C1[6] * pow((l1que.r[is].senz[ip] * satzdir), 2);
5641 
5642  lsst = sstlo + ((lat - dBtlo) / (dBthi - dBtlo)) * (ssthi - sstlo);
5643 
5644  }
5645  }
5646  if (sensorID == VIIRSN) {
5647  lsst = lsst - CtoK; /* internally, want all sst's in Deg C */
5648  }
5649  return (lsst);
5650 }
5651 
5652 /* ----------------------------------------------------------------------------------- */
5653 /* nlsst4() - nonlinear sst4, 4um sea surface temperature (3.9 & 4.0 um) */
5654 /* */
5655 /* B. Franz, SAIC, August 2005. */
5656 
5657 /* ----------------------------------------------------------------------------------- */
5