ocssw V2022
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 = l1_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 || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
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 (l1_input->resolution == 250) {
1585  fullscanpix = 5416;
1586  } else if (l1_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 || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
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 || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
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  case VIIRSJ1:
1934  strcpy(mission, "VIIR");
1935  break;
1936  default:
1937  strcpy(mission, "AQUA");
1938  break;
1939  }
1940 
1941  /* Open the file */
1942  if ((fp = fopen(input->sstcoeffile, "r")) == NULL) {
1943  fprintf(stderr,
1944  "-E- %s line %d: unable to open sst coef file %s for reading\n",
1945  __FILE__, __LINE__, input->sstcoeffile);
1946  exit(1);
1947  }
1948  /* Find month */
1949  leap = (isleap(year) == TRUE ? 1 : 0);
1950  if (tmonth < 0) {
1951  for (tmonth = 11; tmonth >= 0; tmonth--) {
1952  /* day is one based, StartOfMonth is zero based */
1953  if (day > StartOfMonth[leap][tmonth]) {
1954  break;
1955  }
1956  }
1957  }
1958 
1959  if (sensorID == VIIRSN || sensorID == VIIRSJ1) {
1960  /* VIIRS latband coeffs are this format: */
1961  /* sat, start yyyyddd, start hhmmss, end yyyyddd, end hhmmss, min lat, max lat, c0, c1, c2, c3 */
1962  /* VIIR 2012001 000000 2012069 235900 -90 -40 -4.202194 1.0196764 0.002352195 1.206611 */
1963 
1964  /* Form date string */
1965 
1966  /* msec is the scan line start time */
1967  /* should calculate the actual pixel time (extrapolate from msec) for field 57? */
1968  /* modis, viirs, and avhhr: year, day are per scan line, not just time of first scan */
1969  // pasutime = yds2unix(year, day, ((double) (msec)) / 1000.0);
1970  ztime = ydhmsf(pasutime, 'G');
1971  /* ztime is yyyydddhhmmssfff */
1972  strncpy(odatel, ztime, 13);
1973  odatel[13] = '\0';
1974 
1975  /* Loop through to find bounding times */
1976 
1977  indx = 0;
1978  while (fgets(line, 200, fp)) {
1979  if (line[0] == '#') {
1980  /* look for lines with: # variable = value */
1981  if (!(p = strchr(line, '=')))
1982  continue;
1983  p1 = line + 1; /* look for white space after # and before variable name */
1984  while (isspace(*p1))
1985  p1++;
1986  p2 = p - 1; /* look for white space before = and after variable name */
1987  while (isspace(*p2))
1988  p2--;
1989  /* get variable name from input line */
1990  strncpy(name, p1, p2 - p1 + 1);
1991  name[p2 - p1 + 1] = '\0';
1992 
1993  /*
1994  * Parse parameter value string
1995  */
1996  /* start at character after = and ignore white space before value */
1997  p1 = p + 1;
1998  while (isspace(*p1))
1999  p1++;
2000  p2 = p1;
2001  /* look for white space to determine end of value */
2002  while (!isspace(*p2))
2003  p2++;
2004  /* get value from input line */
2005  strncpy(value, p1, p2 - p1);
2006  value[p2 - p1] = '\0';
2007  /*
2008  * Copy value to appropriate variable
2009  */
2010  if (strcmp(name, "sstref_day_offset") == 0) {
2011  *sstrefoffday = (float) atof(value);
2012  gotsstrefoffday = 1;
2013  }
2014  if (strcmp(name, "sstref_night_offset") == 0) {
2015  *sstrefoffnight = (float) atof(value);
2016  gotsstrefoffnight = 1;
2017  }
2018 
2019  } else {
2020  /* read viirs sst latband coeffs */
2021  if (strncmp(line, mission, 4) == 0) {
2022  if (input->viirsnosisaf == 1) {
2023  sscanf(line, "%4s %7s %7s %1s %f %f %f %f %f %f %f",
2024  mission2, sdates, edates, dorn,
2025  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2026  &coef[indx][3], &coef[indx][4], &coef[indx][5],
2027  &coef[indx][6]);
2028  } else if (input->viirsnv7 >= 0) {
2029  /* all latband versions except v6.4.1 */
2030  sscanf(line,
2031  "%4s %7s %6s %7s %6s %f %f %f %f %f %f %f %f %f",
2032  mission2, sdates, stime, edates, etime,
2033  &bounds[indx][0], &bounds[indx][1],
2034  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2035  &coef[indx][3], &coef[indx][4], &coef[indx][5],
2036  &coef[indx][6]);
2037  sprintf(sdatel, "%s%s", sdates, stime);
2038  sprintf(edatel, "%s%s", edates, etime);
2039  if (strcmp(odatel, sdatel) >= 0
2040  && (strcmp(odatel, edatel) <= 0
2041  || strcmp(edates, "0000000") == 0
2042  || strcmp(edates, "") == 0)) {
2043  indx++;
2044  if (indx == 7) {
2045  found = 1;
2046  break;
2047  }
2048  }
2049  } else {
2050  /* viirs nlsst v6.4.1 has extra satz terms, but not mirror side */
2051  /* and is by month, not start and end dates */
2052  sscanf(line, "%4s %d %f %f %f %f %f %f %f %f",
2053  mission2, &month, &bounds[indx][0], &bounds[indx][1],
2054  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2055  &coef[indx][3], &coef[indx][5], &coef[indx][6]);
2056  /* no mirror side for viirs */
2057  coef[indx][4] = 0.0;
2058  if (month == tmonth + 1) {
2059  indx++;
2060  if (indx == 7) {
2061  found = 1;
2062  break;
2063  }
2064  }
2065  }
2066 
2067  }
2068  }
2069  }
2070 
2071  } else {
2072  /* get latband coeffs for AVHRR or MODIS */
2073 
2074  /* Loop through to find 6 sets of coeffs for this month */
2075  indx = 0;
2076  while (fgets(line, 200, fp)) {
2077  if (line[0] == '#') {
2078  /* look for lines with: # variable = value */
2079  if (!(p = strchr(line, '=')))
2080  continue;
2081  p1 = line + 1; /* look for white space after # and before variable name */
2082  while (isspace(*p1))
2083  p1++;
2084  p2 = p - 1; /* look for white space before = and after variable name */
2085  while (isspace(*p2))
2086  p2--;
2087  /* get variable name from input line */
2088  strncpy(name, p1, p2 - p1 + 1);
2089  name[p2 - p1 + 1] = '\0';
2090 
2091  /*
2092  * Parse parameter value string
2093  */
2094  /* start at character after = and ignore white space before value */
2095  p1 = p + 1;
2096  while (isspace(*p1))
2097  p1++;
2098  p2 = p1;
2099  /* look for white space to determine end of value */
2100  while (!isspace(*p2))
2101  p2++;
2102  /* get value from input line */
2103  strncpy(value, p1, p2 - p1);
2104  value[p2 - p1] = '\0';
2105  /*
2106  * Copy value to appropriate variable
2107  */
2108  if (strcmp(name, "sstref_day_offset") == 0) {
2109  *sstrefoffday = (float) atof(value);
2110  gotsstrefoffday = 1;
2111  }
2112  if (strcmp(name, "sstref_night_offset") == 0) {
2113  *sstrefoffnight = (float) atof(value);
2114  gotsstrefoffnight = 1;
2115  }
2116 
2117  } else {
2118  /* read sst latband coeffs for AVHRR or MODIS */
2119  if (strncmp(line, mission, 4) == 0) {
2120  if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
2121  /* terra and aqua nlsst has extra mirror side and satz terms */
2122  sscanf(line,
2123  "%4s %d %f %f %f %f %f %f %f %f %f %d %f %f",
2124  mission2, &month, &bounds[indx][0],
2125  &bounds[indx][1], &coef[indx][0],
2126  &coef[indx][1], &coef[indx][2], &coef[indx][3],
2127  &coef[indx][4], &coef[indx][5], &coef[indx][6],
2128  &tmp1, &tmp2, &tmp3);
2129  } else {
2130  /* don't have extra terms for avhrr nlsst */
2131  sscanf(line, "%4s %d %f %f %f %f %f %f", mission2,
2132  &month, &bounds[indx][0], &bounds[indx][1],
2133  &coef[indx][0], &coef[indx][1], &coef[indx][2],
2134  &coef[indx][3]);
2135  coef[indx][4] = 0.0;
2136  coef[indx][5] = 0.0;
2137  coef[indx][6] = 0.0;
2138  }
2139  if (month == tmonth + 1) {
2140  indx++;
2141  if (indx == 7){
2142  found = 1;
2143  break;
2144  }
2145  }
2146  }
2147  }
2148  }
2149  }
2150 
2151  fclose(fp);
2152 
2153  if (found == 1) {
2154 
2155  if ((sensorID == VIIRSN || sensorID == VIIRSJ1)
2156  && (gotsstrefoffday == 0 || gotsstrefoffnight == 0)) {
2157  fprintf(stderr,
2158  "-E- %s line %d: Day and night sst reference offsets not found in %s\n",
2159  __FILE__, __LINE__, input->sstcoeffile);
2160  exit(1);
2161  }
2162 
2163  printf("Loading SST lat band coefficients from %s:\n",
2164  input->sstcoeffile);
2165 
2166  if ((sensorID == VIIRSN || sensorID == VIIRSJ1) && input->viirsnosisaf == 1) {
2167  /* viirs osi-saf coefs are not latband */
2168  for (indx = 0; indx < 4; indx++) {
2169  printf("%s %s %s %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2170  coeflabel[indx], sdates, edates, coef[indx][0], coef[indx][1], coef[indx][2],
2171  coef[indx][3], coef[indx][4], coef[indx][5], coef[indx][6]);
2172  }
2173  } else {
2174  for (indx = 0; indx < 7; indx++) {
2175  if ((sensorID == VIIRSN || sensorID == VIIRSJ1) && input->viirsnv7 >= 1) {
2176  /* v7 viirs have start and end dates, not months */
2177  printf(
2178  "%s %s %s %s %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2179  sdates, stime, edates, etime, bounds[indx][0], bounds[indx][1],
2180  coef[indx][0], coef[indx][1], coef[indx][2], coef[indx][3],
2181  coef[indx][4], coef[indx][5], coef[indx][6]);
2182  } else {
2183  /* most latband coeffs are by month */
2184  printf(
2185  "%d %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2186  month, bounds[indx][0], bounds[indx][1], coef[indx][0],
2187  coef[indx][1], coef[indx][2], coef[indx][3],
2188  coef[indx][4], coef[indx][5], coef[indx][6]);
2189  }
2190  }
2191  }
2192  printf(" sst reference day offset = %f\n", *sstrefoffday);
2193  printf(" sst reference night offset = %f\n", *sstrefoffnight);
2194 
2195  } else {
2196  fprintf(stderr,
2197  "-E- %s line %d: unable to locate valid SST coefficients for %s in %s\n",
2198  __FILE__, __LINE__, odates, input->sstcoeffile);
2199  exit(1);
2200  }
2201 
2202  return;
2203 }
2204 
2205 /* ----------------------------------------------------------------------------------- */
2206 /* read_sst4_coeff() - reads 4um sst coefficients for the specified date */
2207 /* */
2208 /* B. Franz, SAIC, 11 August 2005 */
2209 
2210 /* ----------------------------------------------------------------------------------- */
2211 void read_sst4_coeff(int32 sensorID, char *filename, double scantime,
2212  float **bounds, float **coef) {
2213  char mission[5];
2214  char mission2[5];
2215 
2216  FILE *fp;
2217  char *line;
2218  line = (char *) calloc(200, sizeof (char));
2219  char odate[8];
2220  int found = 0;
2221  // int tmonth;
2222  int16_t year, day;
2223  double sec;
2224  unix2yds(scantime, &year, &day, &sec);
2225 
2226  int month;
2227  int indx;
2228  // char tmp1[200] = "\0"; /* random length to slurp up whatevers after the coefs */
2229  int32 leap;
2230 
2231  switch (sensorID) {
2232  case MODISA:
2233  strcpy(mission, "AQUA");
2234  break;
2235  case MODIST:
2236  strcpy(mission, "TERR");
2237  break;
2238  case VIIRSN:
2239  case VIIRSJ1:
2240  strcpy(mission, "VIIR");
2241  break;
2242  default:
2243  strcpy(mission, "AQUA");
2244  break;
2245  }
2246 
2247  /* Open the file */
2248  if ((fp = fopen(filename, "r")) == NULL) {
2249  fprintf(stderr, "-E- %s line %d: unable to open %s for reading\n",
2250  __FILE__, __LINE__, filename);
2251  exit(1);
2252  }
2253 
2254  sprintf(odate, "%4d%03d", year, day);
2255 
2256  /* Find month */
2257  leap = (isleap(year) == TRUE ? 1 : 0);
2258  if (tmonth < 0) {
2259  for (tmonth = 11; tmonth >= 0; tmonth--) {
2260  if (day > StartOfMonth[leap][tmonth]) {
2261  break;
2262  }
2263  }
2264  }
2265  /* Loop through to find 6 sets of coeffs for this month */
2266  indx = 0;
2267  fprintf(stderr, " looking for month %d mission %s\n", tmonth + 1, mission);
2268  while (fgets(line, 200, fp)) {
2269  /* read sst4 latband coeffs */
2270  if (strncmp(line, mission, 4) == 0) {
2271  sscanf(line, "%4s %d %f %f %f %f %f %f %f %f %f", &mission2[0],
2272  &month, &bounds[indx][0], &bounds[indx][1], &coef[indx][0],
2273  &coef[indx][1], &coef[indx][2], &coef[indx][3],
2274  &coef[indx][4], &coef[indx][5], &coef[indx][6]);
2275  if (month == tmonth + 1) {
2276  indx++;
2277  if (indx == 7) {
2278  found = 1;
2279  break;
2280  }
2281  }
2282  }
2283  }
2284 
2285  fclose(fp);
2286 
2287  if (found == 1) {
2288 
2289  printf("Loading SST4 lat band coefficients from %s:\n", filename);
2290  for (indx = 0; indx < 7; indx++) {
2291  printf("%d %6.1f %6.1f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f\n",
2292  month, bounds[indx][0], bounds[indx][1], coef[indx][0], coef[indx][1], coef[indx][2],
2293  coef[indx][3], coef[indx][4], coef[indx][5], coef[indx][6]);
2294  }
2295 
2296  } else {
2297  fprintf(stderr,
2298  "-E- %s line %d: unable to locate valid SST4 coefficients for %s in %s\n",
2299  __FILE__, __LINE__, odate, filename);
2300  exit(1);
2301  }
2302  free(line);
2303 
2304  return;
2305 }
2306 
2307 /* ----------------------------------------------------------------------------------- */
2308 /* sstcloud() - uses red band to test homogeneity in nx x ny box (RSMAS) */
2309 /* */
2310 /* B. Franz, SAIC, August 2005. */
2311 
2312 /* ----------------------------------------------------------------------------------- */
2313 int sstcloud(int32_t ip, int32_t nx, int32_t ny, float thresh) {
2314  extern l1qstr l1que;
2315 
2316  int32_t nscan = l1que.nq;
2317  int32_t npix = l1que.r[0].npix;
2318  int32_t is = nscan / 2;
2319  int32_t ip1, ip2;
2320  int32_t is1, is2;
2321  int32_t i, j;
2322  int32_t ipb;
2323  float rhot;
2324  int flag = 0;
2325  int cnt = 0;
2326  float maxv = BAD_FLT;
2327  float minv = -1.0 * BAD_FLT; /* initial minimum has to be large positive number */
2328 
2329  if (!haveRed)
2330  return (flag);
2331 
2332  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2333  if (nscan < ny) {
2334  printf(
2335  "-E- %s line %d: L1 queue size of %d is too small for requested homogeneity test of %d x %d.\n",
2336  __FILE__, __LINE__, l1que.nq, nx, ny);
2337  exit(1);
2338  }
2339 
2340  /* algorithm is only valid for daytime and in non glint area */
2341  if (l1que.r[is].solz[ip] >= solznight
2342  || l1que.r[is].glint_coef[ip] > glintmax)
2343  return (flag);
2344 
2345  /* compute pix/scan limits for the Row Of Interest (ROI) */
2346  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2347  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2348  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2349  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2350  // filt_dist = -nx / 2;
2351  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2352  // ipmod -= l1que.r[is].spix;
2353  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2354  //
2355  // filt_dist = nx / 2;
2356  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2357  // ipmod -= l1que.r[is].spix;
2358  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2359  // } else {
2360  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2361  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2362  // }
2363 
2364  /* compute max and min normalized reflectance in window */
2365  /* l1_hmodis_hdf.c converts visible bands, i.e. red, to radiance */
2366  /* so we need to convert it to reflectance */
2367  /* atmospheric transmittance and path radiance due to scattering and attenuation by
2368  * atmospheric gases and aerosols is a problem for ocean color, but not for sst
2369  * so, while we probably don't need the extra terms here, we definitely don't need
2370  * them for the cloud decision trees */
2371  for (j = is1; j <= is2; j++)
2372  for (i = ip1; i <= ip2; i++) {
2373  ipb = i * nbvis + ibred;
2374  if (l1que.r[j].Lt[ipb] > 0.0 && l1que.r[j].Lt[ipb] < satred) {
2375  rhot = PI * l1que.r[j].Lt[ipb] / l1que.r[j].Fo[ibred]
2376  / l1que.r[j].tg_sol[ipb] / l1que.r[j].tg_sen[ipb]
2377  / l1que.r[j].t_sol[ipb] / l1que.r[j].t_sen[ipb]
2378  / cos(l1que.r[j].solz[ip] / RADEG);
2379  maxv = MAX(maxv, rhot);
2380  minv = MIN(minv, rhot);
2381  cnt++;
2382  }
2383  }
2384 
2385  /* added for collect 5 (if all saturated but 1, then cloud) */
2386  if (cnt < 2 || (maxv - minv) > thresh) {
2387  flag = 1;
2388  }
2389 
2390  return (flag);
2391 }
2392 
2393 /* ----------------------------------------------------------------------------------- */
2394 /* rhoCboxstats() - test homogeneity in nx x ny box of rho_cirrus (RSMAS) */
2395 
2396 /* ----------------------------------------------------------------------------------- */
2397 int32_t rhoCboxstats(int32_t ip, int32_t nx, int32_t ny, statstr *statrec) {
2398  extern l1qstr l1que;
2399 
2400  static int32_t len = 0;
2401  int32_t nscan = l1que.nq;
2402  int32_t npix = l1que.r[0].npix;
2403  int32_t ip1, ip2;
2404  int32_t is1, is2;
2405  int32_t is, i, j, ix;
2406  float Bt;
2407  static float *x = NULL;
2408 
2409 
2410  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2411  if (nscan < ny) {
2412  printf("-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2413  __FILE__, __LINE__, l1que.nq, nx, ny);
2414  exit(1);
2415  }
2416 
2417  /* allocate sufficient workspace for the filter size */
2418  if (x == NULL || len < nx * ny) {
2419  len = nx*ny;
2420  if (x != NULL) free(x);
2421  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2422  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2423  __FILE__, __LINE__);
2424  exit(1);
2425  }
2426  }
2427 
2428  /* Compute queue scan limits for the Row Of Interest (ROI) */
2429  is = nscan / 2;
2430  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2431  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2432 
2433  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2434  /* for aggregated VIIRS, limits are aggregation zone dependent */
2435  // if( ( l1que.r[is].l1file->sensorID == VIIRS ) && ( l1que.r[is].scn_fmt == 0 ) ) {
2436  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2437  // filt_dist = -nx / 2;
2438  // viirs_pxcvt_agdel( ip_scan, filt_dist, &ipmod );
2439  // ipmod -= l1que.r[is].spix;
2440  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2441  //
2442  // filt_dist = nx / 2;
2443  // viirs_pxcvt_agdel( ip_scan, filt_dist, &ipmod );
2444  // ipmod -= l1que.r[is].spix;
2445  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2446  // } else {
2447  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2448  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2449  // }
2450 
2451  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2452  statrec->min = -1.0 * BAD_FLT;
2453  statrec->max = BAD_FLT;
2454  statrec->avg = BAD_FLT;
2455  statrec->med = BAD_FLT;
2456  statrec->cen = BAD_FLT;
2457  statrec->sqr = BAD_FLT;
2458  statrec->sd = BAD_FLT;
2459  statrec->cnt = 0;
2460  for (ix = 0; ix < len; ix++) {
2461  x[ix] = 0.0;
2462  }
2463 
2464  /* compute stats for window */
2465  for (j = is1; j <= is2; j++) for (i = ip1; i <= ip2; i++) {
2466  Bt = l1que.r[j].rho_cirrus[i];
2467  if (Bt > (BAD_FLT + 0.1)) {
2468  statrec->max = MAX(statrec->max, Bt);
2469  statrec->min = MIN(statrec->min, Bt);
2470  if (statrec->cnt == 0) {
2471  /* first value - start the summing for the average at zero */
2472  statrec->avg = 0.0;
2473  }
2474  statrec->avg += Bt;
2475  statrec->sqr += Bt * Bt;
2476  x[statrec->cnt] = Bt;
2477  statrec->cnt++;
2478  }
2479  }
2480  if (statrec->cnt > 2) {
2481  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2482  statrec->max = x[statrec->cnt - 1];
2483  statrec->min = x[0];
2484  if (statrec->cnt > 1) {
2485  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2486  if (statrec->sd > 0.0) {
2487  statrec->sd = sqrt(statrec->sd);
2488  } else {
2489  statrec->sd = BAD_FLT;
2490  }
2491  }
2492 
2493  statrec->med = x[statrec->cnt / 2];
2494  }
2495  if (statrec->cnt > 0)
2496  statrec->avg /= statrec->cnt;
2497 
2498  return (statrec->cnt);
2499 }
2500 
2501 /* ----------------------------------------------------------------------------------- */
2502 /* rhotboxstats() - test homogeneity in nx x ny box of rhot's (RSMAS) */
2503 /* */
2504 
2505 /* ----------------------------------------------------------------------------------- */
2506 int32_t rhotboxstats(int32_t ip, int ib, int nbands, int32_t nx, int32_t ny,
2507  statstr *statrec) {
2508  extern l1qstr l1que;
2509 
2510  static int32_t len = 0;
2511  int32_t nscan = l1que.nq;
2512  int32_t npix = l1que.r[0].npix;
2513  int32_t ip1, ip2;
2514  int32_t is1, is2;
2515  int32_t is, i, j, ix;
2516  int32_t ipb;
2517  // int32_t ip_scan, filt_dist, ipmod;
2518  float rhot;
2519  static float *x = NULL;
2520 
2521  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2522  if (nscan < ny) {
2523  printf(
2524  "-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2525  __FILE__, __LINE__, l1que.nq, nx, ny);
2526  exit(1);
2527  }
2528 
2529  /* allocate sufficient workspace for the filter size */
2530  if (x == NULL || len < nx * ny) {
2531  len = nx*ny;
2532  if (x != NULL) free(x);
2533  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2534  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2535  __FILE__, __LINE__);
2536  exit(1);
2537  }
2538  }
2539 
2540  /* Compute queue scan limits for the Row Of Interest (ROI) */
2541  is = nscan / 2;
2542  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2543  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2544 
2545  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2546  /* for aggregated VIIRS, limits are aggregation zone dependent */
2547  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2548  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2549  // filt_dist = -nx / 2;
2550  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2551  // ipmod -= l1que.r[is].spix;
2552  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2553  //
2554  // filt_dist = nx / 2;
2555  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2556  // ipmod -= l1que.r[is].spix;
2557  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2558  // } else {
2559  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2560  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2561  // }
2562 
2563  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2564  statrec->min = -1.0 * BAD_FLT;
2565  statrec->max = BAD_FLT;
2566  statrec->avg = BAD_FLT;
2567  statrec->med = BAD_FLT;
2568  statrec->cen = BAD_FLT;
2569  statrec->sqr = BAD_FLT;
2570  statrec->sd = BAD_FLT;
2571  statrec->cnt = 0;
2572  for (ix = 0; ix < len; ix++) {
2573  x[ix] = 0.0;
2574  }
2575 
2576  /* compute stats for window */
2577  for (j = is1; j <= is2; j++)
2578  for (i = ip1; i <= ip2; i++) {
2579  ipb = i * nbands + ib;
2580  /* saturated Lt is 1000.0 or sstbad=-32767 (class bad could be 0.0, but don't care?) */
2581  if (l1que.r[j].Lt[ipb] > sstbad + 1.0 && l1que.r[j].Lt[ipb] < 1000.0 - 1.0) {
2582  rhot = M_PI * l1que.r[j].Lt[ipb] / l1que.r[j].Fo[ib] / l1que.r[j].csolz[i];
2583  statrec->max = MAX(statrec->max, rhot);
2584  statrec->min = MIN(statrec->min, rhot);
2585  if (statrec->cnt == 0) {
2586  /* first value - start the summing for the average at zero */
2587  statrec->avg = 0.0;
2588  }
2589  statrec->avg += rhot;
2590  statrec->sqr += rhot * rhot;
2591  x[statrec->cnt] = rhot;
2592  statrec->cnt++;
2593  }
2594  }
2595  if (statrec->cnt > 2) {
2596  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2597  statrec->max = x[statrec->cnt - 1];
2598  statrec->min = x[0];
2599  if (statrec->cnt > 1) {
2600  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2601  if (statrec->sd > 0.0) {
2602  statrec->sd = sqrt(statrec->sd);
2603  } else {
2604  statrec->sd = BAD_FLT;
2605  }
2606  }
2607 
2608  statrec->med = x[statrec->cnt / 2];
2609  }
2610  if (statrec->cnt > 0)
2611  statrec->avg /= statrec->cnt;
2612 
2613  return (statrec->cnt);
2614 }
2615 
2616 /* ----------------------------------------------------------------------------------- */
2617 /* btboxstats() - test homogeneity in nx x ny box of Bt's (RSMAS) */
2618 /* */
2619 /* B. Franz, SAIC, August 2005. */
2620 
2621 /* ----------------------------------------------------------------------------------- */
2622 int32_t btboxstats(int32_t ip, int ib, int nbands, int32_t nx, int32_t ny,
2623  statstr *statrec) {
2624  extern l1qstr l1que;
2625 
2626  static int32_t len = 0;
2627  int32_t nscan = l1que.nq;
2628  int32_t npix = l1que.r[0].npix;
2629  int32_t ip1, ip2;
2630  int32_t is1, is2;
2631  int32_t is, i, j, ix;
2632  int32_t ipb;
2633  float Bt;
2634  static float *x = NULL;
2635 
2636  /* make sure code is not inconsistent NQMIN in l12_parms.h */
2637  if (nscan < ny) {
2638  printf(
2639  "-E- %s line %d: L1 queue size of %d is too small for %d x %d box stats.\n",
2640  __FILE__, __LINE__, l1que.nq, nx, ny);
2641  exit(1);
2642  }
2643 
2644  /* allocate sufficient workspace for the filter size */
2645  if (x == NULL || len < nx * ny) {
2646  len = nx*ny;
2647  if (x != NULL) free(x);
2648  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2649  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2650  __FILE__, __LINE__);
2651  exit(1);
2652  }
2653  }
2654 
2655  /* Compute queue scan limits for the Row Of Interest (ROI) */
2656  is = nscan / 2;
2657  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2658  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2659 
2660  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2661  /* for aggregated VIIRS, limits are aggregation zone dependent */
2662  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2663  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2664  // filt_dist = -nx / 2;
2665  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2666  // ipmod -= l1que.r[is].spix;
2667  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2668  //
2669  // filt_dist = nx / 2;
2670  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2671  // ipmod -= l1que.r[is].spix;
2672  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2673  // } else {
2674  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2675  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2676  // }
2677 
2678  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2679  statrec->min = -1.0 * BAD_FLT;
2680  statrec->max = BAD_FLT;
2681  statrec->avg = BAD_FLT;
2682  statrec->med = BAD_FLT;
2683  statrec->cen = BAD_FLT;
2684  statrec->sqr = BAD_FLT;
2685  statrec->sd = BAD_FLT;
2686  statrec->cnt = 0;
2687  for (ix = 0; ix < len; ix++) {
2688  x[ix] = 0.0;
2689  }
2690 
2691  /* compute stats for window */
2692  for (j = is1; j <= is2; j++) {
2693  /* don't use values from aqua detector zero Bt40 */
2694  /* Peter and Kay say to not use average, just skip detector zero */
2695  if ((l1que.r[j].l1file->sensorID != MODISA)
2696  || (ib != ib40)
2697  || (l1que.r[j].detnum != 0)) {
2698 
2699  for (i = ip1; i <= ip2; i++) {
2700  ipb = i * nbands + ib;
2701  Bt = l1que.r[j].Bt[ipb];
2702  if (Bt > BT_LO + 0.1 && Bt < BT_HI - 0.1) {
2703  statrec->max = MAX(statrec->max, Bt);
2704  statrec->min = MIN(statrec->min, Bt);
2705  if (statrec->cnt == 0) {
2706  /* first value - start the summing for the average at zero */
2707  statrec->avg = 0.0;
2708  statrec->sqr = 0.0;
2709  }
2710  statrec->avg += Bt;
2711  statrec->sqr += Bt * Bt;
2712  x[statrec->cnt] = Bt;
2713  statrec->cnt++;
2714  }
2715  }
2716  }
2717  }
2718  if (statrec->cnt > 2) {
2719  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2720  statrec->max = x[statrec->cnt - 1];
2721  statrec->min = x[0];
2722  if (statrec->cnt > 1) {
2723  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2724  if (statrec->sd > 0.0) {
2725  statrec->sd = sqrt(statrec->sd);
2726  } else {
2727  statrec->sd = BAD_FLT;
2728  }
2729  }
2730 
2731  statrec->med = x[statrec->cnt / 2];
2732  }
2733  if (statrec->cnt > 0) {
2734  statrec->avg /= statrec->cnt;
2735  }
2736 
2737  return (statrec->cnt);
2738 }
2739 
2740 /* ----------------------------------------------------------------------------------- */
2741 /* sstboxstats() - test homogeneity in nx x ny box of Bt's (RSMAS) */
2742 /* */
2743 /* B. Franz, SAIC, August 2005. */
2744 
2745 /* ----------------------------------------------------------------------------------- */
2746 int32_t sstboxstats(int32_t ip, int32_t nx, int32_t ny, statstr *statrec) {
2747  extern l1qstr l1que;
2748 
2749  static int32_t len = 0;
2750  int32_t nscan = l1que.nq;
2751  int32_t npix = l1que.r[0].npix;
2752  int32_t ip1, ip2;
2753  int32_t is1, is2;
2754  int32_t is, i, j, ix;
2755  // int32_t ip_scan, filt_dist, ipmod;
2756  float sst;
2757  static float *x = NULL;
2758 
2759  /* make sure code is not inconsistent */
2760  if (nscan < ny) {
2761  printf(
2762  "-E- %s line %d: sst queue size of %d is too small for %d x %d box stats.\n",
2763  __FILE__, __LINE__, nscan, nx, ny);
2764  exit(1);
2765  }
2766 
2767  /* allocate sufficient workspace for the filter size */
2768  if (x == NULL || len < nx * ny) {
2769  len = nx*ny;
2770  if (x != NULL) free(x);
2771  if ((x = (float *) malloc(len * sizeof (float))) == NULL) {
2772  printf("-E- %s line %d: Unable to allocate workspace for median and stdev\n",
2773  __FILE__, __LINE__);
2774  exit(1);
2775  }
2776  }
2777 
2778  /* Compute queue scan limits for the Row Of Interest (ROI) */
2779  is = nscan / 2;
2780  is1 = MIN(MAX(0, is - ny / 2), nscan - 1);
2781  is2 = MAX(MIN(nscan - 1, is + ny / 2), 0);
2782 
2783  /* compute pixel neighbor limits for the Row Of Interest (ROI) */
2784  /* for aggregated VIIRS, limits are aggregation zone dependent */
2785  // if ((l1que.r[is].l1file->sensorID == VIIRS) && (l1que.r[is].scn_fmt == 0)) {
2786  // ip_scan = ip + l1que.r[is].spix; /* scan pixel */
2787  // filt_dist = -nx / 2;
2788  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2789  // ipmod -= l1que.r[is].spix;
2790  // ip1 = MIN( MAX( 0, ipmod ), npix-1 );
2791  //
2792  // filt_dist = nx / 2;
2793  // viirs_pxcvt_agdel(ip_scan, filt_dist, &ipmod);
2794  // ipmod -= l1que.r[is].spix;
2795  // ip2 = MAX( MIN( npix-1, ipmod ), 0 );
2796  // } else {
2797  ip1 = MIN(MAX(0, ip - nx / 2), npix - 1);
2798  ip2 = MAX(MIN(npix - 1, ip + nx / 2), 0);
2799  // }
2800 
2801  /* initialize output rec - only min,max,avg,cnt are actually used here*/
2802  statrec->min = -1.0 * BAD_FLT;
2803  statrec->max = BAD_FLT;
2804  statrec->avg = BAD_FLT;
2805  statrec->med = BAD_FLT;
2806  statrec->cen = BAD_FLT;
2807  statrec->sqr = BAD_FLT;
2808  statrec->sd = BAD_FLT;
2809  statrec->cnt = 0;
2810  for (ix = 0; ix < len; ix++) {
2811  x[ix] = 0.0;
2812  }
2813 
2814  /* compute stats for window */
2815  for (j = is1; j <= is2; j++) {
2816  for (i = ip1; i <= ip2; i++) {
2817  sst = sstq[j][i];
2818  if (sst >= SSTmin && sst <= SSTmax) {
2819  statrec->max = MAX(statrec->max, sst);
2820  statrec->min = MIN(statrec->min, sst);
2821  if (statrec->cnt == 0) {
2822  /* first value - start the summing for the average at zero */
2823  statrec->avg = 0.0;
2824  statrec->sqr = 0.0;
2825  }
2826  statrec->avg += sst;
2827  statrec->sqr += sst * sst;
2828  x[statrec->cnt] = sst;
2829  statrec->cnt++;
2830  }
2831  }
2832  }
2833  if (statrec->cnt > 2) {
2834  qsort(x, statrec->cnt, sizeof (float), (int (*)(const void *, const void *)) compfloat);
2835  statrec->max = x[statrec->cnt - 1];
2836  statrec->min = x[0];
2837  if (statrec->cnt > 1) {
2838  statrec->sd = (statrec->sqr - statrec->avg * statrec->avg / statrec->cnt) / (statrec->cnt - 1);
2839  if (statrec->sd > 0.0) {
2840  statrec->sd = sqrt(statrec->sd);
2841  } else {
2842  statrec->sd = BAD_FLT;
2843  }
2844  }
2845 
2846  statrec->med = x[statrec->cnt / 2];
2847  }
2848  if (statrec->cnt > 0) {
2849  statrec->avg /= statrec->cnt;
2850  }
2851 
2852  return (statrec->cnt);
2853 }
2854 /* ----------------------------------------------------------------------------------- */
2855 /* run_sstboxstdev() - calculate the sst standard deviation */
2856 
2857 /* ----------------------------------------------------------------------------------- */
2858 void run_sstboxstdev(int npix, float *sst_stdev) {
2859 
2860  int ip;
2861  statstr statrec;
2862  for (ip = 0; ip < npix; ip++) {
2863  sst_stdev[ip] = sstbad;
2864  if (sstboxstats(ip, btbox, btbox, &statrec) > 0) {
2865  sst_stdev[ip] = statrec.sd;
2866  }
2867  }
2868 }
2869 
2870 /* ----------------------------------------------------------------------------------- */
2871 /* run_rhotboxmin() - initialize the static min arrays */
2872 
2873 /* ----------------------------------------------------------------------------------- */
2874 void run_rhotboxmin(int npix, int btidx, int nbands, float *minarr) {
2875 
2876  int ip;
2877  statstr statrec;
2878  for (ip = 0; ip < npix; ip++) {
2879  minarr[ip] = sstbad;
2880  if (rhotboxstats(ip, btidx, nbands, btbox, btbox, &statrec) > 0) {
2881  minarr[ip] = statrec.min;
2882  }
2883  }
2884 }
2885 
2886 /* ----------------------------------------------------------------------------------- */
2887 /* run_rhotboxmaxmin() - initialize the static maxmin arrays */
2888 
2889 /* ----------------------------------------------------------------------------------- */
2890 void run_rhotboxmaxmin(int npix, int btidx, int nbands, float *maxminarr, float *maxarr) {
2891 
2892  int ip;
2893  statstr statrec;
2894  for (ip = 0; ip < npix; ip++) {
2895  maxminarr[ip] = sstbad;
2896  maxarr[ip] = sstbad;
2897  if (rhotboxstats(ip, btidx, nbands, btbox, btbox, &statrec) > 0) {
2898  maxminarr[ip] = statrec.max - statrec.min;
2899  maxarr[ip] = statrec.max;
2900  }
2901  }
2902 }
2903 
2904 /* ----------------------------------------------------------------------------------- */
2905 /* run_btboxavg() - initialize the static avg arrays */
2906 
2907 /* ----------------------------------------------------------------------------------- */
2908 void run_btboxavg(int npix, int btidx, float *avgarr) {
2909 
2910  int ip;
2911  statstr statrec;
2912  for (ip = 0; ip < npix; ip++) {
2913  avgarr[ip] = sstbad;
2914  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2915  avgarr[ip] = statrec.avg;
2916  }
2917  }
2918 }
2919 
2920 /* ----------------------------------------------------------------------------------- */
2921 /* run_btboxmaxmin() - initialize the static maxmin arrays */
2922 
2923 /* ----------------------------------------------------------------------------------- */
2924 void run_btboxmaxmin(int npix, int btidx, float *maxminarr) {
2925 
2926  int ip;
2927  statstr statrec;
2928  for (ip = 0; ip < npix; ip++) {
2929  maxminarr[ip] = sstbad;
2930  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2931  maxminarr[ip] = statrec.max - statrec.min;
2932  }
2933  }
2934 }
2935 /* ----------------------------------------------------------------------------------- */
2936 /* run_btboxmin() - initialize the static min arrays */
2937 
2938 /* ----------------------------------------------------------------------------------- */
2939 void run_btboxmin(int npix, int btidx, float *minarr) {
2940 
2941  int ip;
2942  statstr statrec;
2943  for (ip = 0; ip < npix; ip++) {
2944  minarr[ip] = sstbad;
2945  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2946  minarr[ip] = statrec.min;
2947  }
2948  }
2949 }
2950 /* ----------------------------------------------------------------------------------- */
2951 /* run_btboxmax() - initialize the static max arrays */
2952 
2953 /* ----------------------------------------------------------------------------------- */
2954 void run_btboxmax(int npix, int btidx, float *maxarr) {
2955 
2956  int ip;
2957  statstr statrec;
2958  for (ip = 0; ip < npix; ip++) {
2959  maxarr[ip] = sstbad;
2960  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2961  maxarr[ip] = statrec.max;
2962  }
2963  }
2964 }
2965 
2966 void run_rhoCboxmaxmin(int npix, float *maxminarr, float *minarr, float *maxarr) {
2967 
2968  int ip;
2969  statstr statrec;
2970  for (ip = 0; ip < npix; ip++) {
2971  maxminarr[ip] = sstbad;
2972  minarr[ip] = sstbad;
2973  maxarr[ip] = sstbad;
2974  if (rhoCboxstats(ip, btbox, btbox, &statrec) > 0) {
2975  maxminarr[ip] = statrec.max - statrec.min;
2976  minarr[ip] = statrec.min;
2977  maxarr[ip] = statrec.max;
2978  }
2979  }
2980 }
2981 
2982 void run_btboxstdev(int npix, int btidx, float *stdevarr) {
2983  int ip;
2984  statstr statrec;
2985  for (ip = 0; ip < npix; ip++) {
2986  stdevarr[ip] = sstbad;
2987  if (btboxstats(ip, btidx, NBANDSIR, btbox, btbox, &statrec) > 0) {
2988  stdevarr[ip] = statrec.sd;
2989  }
2990  }
2991 }
2992 /* ----------------------------------------------------------------------------------- */
2993 /* btavg() - Average detector 9 with detector 1 to replace aqua channel 23 detector 0 */
2994 /* */
2995 /* S. Walsh, RSMAS, August 2013. */
2996 
2997 /* ----------------------------------------------------------------------------------- */
2998 int32_t btavg(int32_t ip, int is, int ib, int nbands, statstr *statrec) {
2999  extern l1qstr l1que;
3000 
3001  int32_t nscan = l1que.nq;
3002  int32_t is1, is2;
3003  int32_t j;
3004  int32_t ipb;
3005  float Bt;
3006 
3007  /* this routine should only be used to fake Bt40's for aqua detector zero */
3008  if ((l1que.r[is].l1file->sensorID != MODISA)
3009  || (ib != ib40)
3010  || (l1que.r[is].detnum != 0)) {
3011  printf(
3012  "-E- %s line %d: btavg should only be used for AQUA detector zero Bt40.\n",
3013  __FILE__, __LINE__);
3014  exit(1);
3015  }
3016 
3017  /* make sure code is not inconsistent NQMIN in l12_parms.h */
3018  /* Right now modisa uses 3x3 boxes for aqua.
3019  * When we switch to 5x5 then modisa/msl12_filter.dat will have to change to keep a 7x7 box
3020  * so that we have the line before and after the 5x5 box to calculate the Bt40 average for
3021  * detector 0 in whichever line of the box is being calculated
3022  */
3023  if ((nscan < (btbox + 2)) && (is == 0 || is == btbox)) {
3024  printf(
3025  "-E- %s line %d: L1 queue size of %d is too small for 3 line Bt40 average around line %d.\n",
3026  __FILE__, __LINE__, l1que.nq, is);
3027  exit(1);
3028  }
3029 
3030  /* Compute queue scan limits for the Row Of Interest (ROI) in 'is' */
3031  is1 = MIN(MAX(0, is - 1), nscan - 1);
3032  is2 = MAX(MIN(nscan - 1, is + 1), 0);
3033 
3034  /* initialize output rec - only min,max,avg,cnt are actually used here*/
3035  statrec->min = -1.0 * BAD_FLT;
3036  statrec->max = BAD_FLT;
3037  statrec->avg = BAD_FLT;
3038  statrec->med = BAD_FLT;
3039  statrec->cen = BAD_FLT;
3040  statrec->sqr = BAD_FLT;
3041  statrec->sd = BAD_FLT;
3042  statrec->cnt = 0;
3043 
3044  /* compute stats for window */
3045  for (j = is1; j <= is2; j += 2) {
3046  /* don't use bad middle line - aqua detector zero channel 23 */
3047  if (j != is) {
3048  ipb = ip * nbands + ib;
3049  Bt = l1que.r[j].Bt[ipb];
3050  if (Bt > BT_LO + 0.1 && Bt < BT_HI - 0.1) {
3051  if (statrec->cnt == 0) {
3052  /* first value - start the summing for the average at zero */
3053  statrec->avg = 0.0;
3054  }
3055  statrec->avg += Bt;
3056  statrec->cnt++;
3057  }
3058  }
3059  }
3060  if (statrec->cnt > 0)
3061  statrec->avg /= statrec->cnt;
3062 
3063  return (statrec->cnt);
3064 }
3065 
3066 /* ----------------------------------------------------------------------------------- */
3067 /* sstmasked() - returns 1 if pixel was already masked (SST processing skipped) */
3068 
3069 /* ----------------------------------------------------------------------------------- */
3070 int sstmasked(int32_t *flags, int32_t ip) {
3071  if ((flags[ip] & LAND) != 0 || (flags[ip] & NAVFAIL) != 0)
3072 
3073  return (1);
3074  else
3075  return (0);
3076 }
3077 
3078 /* ------------------------------------------------------------------- */
3079 /* avhrr_ascend() - compare center lat of this line and previous line */
3080 
3081 /* ------------------------------------------------------------------- */
3082 int32_t avhrr_ascend(int32_t ny, float *diflat) {
3083  extern l1qstr l1que;
3084 
3085  int32_t nscan = l1que.nq;
3086  int32_t is, cpix;
3087  int32_t jj;
3088  //float diflat;
3089  float diflats[FILTMAX];
3090 
3091  /* make sure code is not inconsistent NQMIN in l12_parms.h */
3092  if (nscan < ny) {
3093  printf(
3094  "-E- %s line %d: L1 queue size of %d is too small for %d line box stats.\n",
3095  __FILE__, __LINE__, l1que.nq, ny);
3096  exit(-1); /* exit value changed from 1 to -1 */
3097  }
3098 
3099  /* current row is center of box */
3100  is = nscan / 2; /* is =1 for boxsize (nscan=ny=) 3 */
3101 
3102  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
3103  /* 204/409 or 1023/2048 or 1599/3200 (zero based) */
3104  /* pixel number start at 0 so first half of scan lines are:
3105  * 0..204, 0..1023, 0..1599 */
3106  cpix = (fullscanpix - 1) / 2;
3107  // if (l1que.r[is - 1].lat[cpix] < l1que.r[is].lat[cpix]) {
3108  // }
3109  // 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)
3110  // Their nadir lat's don't always increase when approaching the north pole, sometimes they decrease from one scan to the next.
3111  // so check for average difference between nadir lats to find split between ascending and descending
3112  // average didn't work. in V2016183000600, which ascends near the dateline and descends over Russia,
3113  // every 16th line (the difference between detector 15 and detector 0) is a greater change in latitude
3114  // than the difference between the lines in one scan. Hard to describe, but easy to see if plot clat[2900:3246]
3115  // Try assuming flip from ascending to descending (and reverse in the south) is only when two in a row match.
3116  *diflat = 0.0;
3117  //*diflat=*diflat+l1que.r[jj].lat[cpix] - l1que.r[jj+1].lat[cpix];
3118  //*diflat=*diflat/(nscan-1);
3119  // check the diff in clat between the lines in the box
3120  for (jj = 0; jj < nscan - 1; jj++) {
3121  diflats[jj] = l1que.r[jj].lat[cpix] - l1que.r[jj + 1].lat[cpix];
3122  }
3123  // check for no change in direction
3124  if ((diflats[is - 1] > 0.0 && diflats[is] > 0.0) ||
3125  (diflats[is - 1] < 0.0 && diflats[is] < 0.0)) {
3126  *diflat = diflats[is];
3127  } else {
3128  // there was a change. But really only if the next dif matches the same direction, otherwise it's just the flaky
3129  // change between swaths (detector 15 to detector 0)
3130  if ((diflats[is] > 0.0 && diflats[is + 1] > 0.0) ||
3131  (diflats[is] < 0.0 && diflats[is + 1] < 0.0)) {
3132  *diflat = diflats[is];
3133  } else {
3134  // not really a change, return the previous diff (or next, but not this one)
3135  *diflat = diflats[is - 1];
3136  }
3137  }
3138 
3139  //if (l1que.cscan == 156) {
3140  // printf(" in ascend, diflat avg = %f\n", *diflat);
3141  //}
3142  if (*diflat < 0.0) {
3143  /* satellite is ascending */
3144  return (1);
3145  } else {
3146  /* satellite is not ascending (usually descending) */
3147  return (0);
3148  }
3149 }
3150 
3151 /* ----------------------------------------------------------------------------------- */
3152 /* set_flags_sst() - set quality flags for long wave sea surface temperature */
3153 /* */
3154 /* B. Franz, SAIC, August 2005. */
3155 
3156 /* ----------------------------------------------------------------------------------- */
3157 void set_flags_sst(l2str *l2rec) {
3158 
3159  extern l1qstr l1que;
3160  int32_t npix = l2rec->l1rec->npix;
3161  int32_t cpix;
3162  int32_t ip, ipb, ipbir;
3163  int32_t yyyyddd;
3164  int32_t ASCEND;
3165  int32_t v5viirs;
3166  float LtRED;
3167  float LtNIR7; // change m6 to LtNIR7, as it's VIIRS 746nm;
3168  float rhoCirrus; // change m9 to rhoCirrus (VIIRS 1.38um);
3169  float Lt16; // change m10 to Lt16, as it's VIIRS 1.61 micron band;
3170  float LtNIR8;
3171  float Bt37;
3172  float Bt39;
3173  float Bt40;
3174  float Bt67;
3175  float Bt85; /* aqua trees - don't use 67, 73, 85 because of cross talk, until fixed 2016? */
3176  float Bt11;
3177  float Bt12;
3178  float dSST_ref; /* sst - ref */
3179  float dSST_SST3; /* sst - sst3 */
3180  float dSST_SST4; /* sst - sst4 */
3181  float dBt_11_12, dBt_37_11, dBt_37_12, dBt_11_37, dBt_67_11; //brightness temp diffs
3182  float dBt_40_11, dBt_85_11;
3183  float Tdeflong;
3184  float subsolar, xdoy, xrad;
3185  float diflat;
3186  statstr statrec;
3187 
3188  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
3189  cpix = (fullscanpix - 1) / 2; /* 204/409 or 1023/2048 (zero based) */
3190 
3191  get_toa_refl(l2rec, ibred, rhotRED);
3192  get_toa_refl(l2rec, ib07, rhotNIR7);
3193  get_toa_refl(l2rec, ib16, rhot16);
3194 
3195  /* check each pixel in scan */
3196  for (ip = 0; ip < npix; ip++) {
3197 
3198  /* SST not processed */
3199  if (sstmasked(l2rec->l1rec->flags, ip)) {
3200  flags_sst[ip] |= SSTF_ISMASKED;
3201  continue;
3202  }
3203 
3204  ipb = ip * l2rec->l1rec->l1file->nbands;
3205  ipbir = ip * NBANDSIR;
3206 
3207  treesum[ip] = 0.0;
3208 
3209  Bt37 = l2rec->l1rec->Bt[ipbir + ib37]; /* modis chan 20, avhrr cen3, viirs SVM12 */
3210  Bt85 = l2rec->l1rec->Bt[ipbir + ib85]; /* modis chan 29, viirs SVM14 */
3211  Bt11 = l2rec->l1rec->Bt[ipbir + ib11]; /* modis chan 31, avhrr cen4, viirs SVM15 */
3212  Bt12 = l2rec->l1rec->Bt[ipbir + ib12]; /* modis chan 32, avhrr cen5, viirs SVM16 */
3213 
3214  LtRED = l2rec->l1rec->Lt[ipb + ibred];
3215 
3216  // some bands really are sensor specific...
3217  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
3218  LtNIR8 = l2rec->l1rec->Lt[ipb + ib08];
3219  } else {
3220  rhoCirrus = l2rec->l1rec->rho_cirrus[ip];
3221  if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
3222  Bt39 = l2rec->l1rec->Bt[ipbir + ib39];
3223  /* don't use aqua Bt40 detector 0, average the previous and next scans instead */
3224  if (l2rec->l1rec->l1file->sensorID == MODISA && l2rec->l1rec->detnum == 0) {
3225  //Bt40 = Bt40_avg[ip];
3226  /* current row is center of box */
3227  int32_t is = l1que.nq / 2;
3228  if (btavg(ip, is, ib40, NBANDSIR, &statrec) > 0) {
3229  Bt40 = statrec.avg;
3230  }
3231  } else
3232  Bt40 = l2rec->l1rec->Bt[ipbir + ib40];
3233  Bt67 = l2rec->l1rec->Bt[ipbir + ib67];
3234  } else if (l2rec->l1rec->l1file->sensorID == VIIRSN || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
3235  Bt40 = l2rec->l1rec->Bt[ipbir + ib40];
3236  LtNIR7 = l2rec->l1rec->Lt[ipb + ib07];
3237  Lt16 = l2rec->l1rec->Lt[ipb + ib16];
3238  }
3239  }
3240 
3241  /* BT could not be computed (radiance out-of-range) */
3242  if (Bt11 < BT_LO + 0.1 || Bt11 > BT_HI - 0.1 || Bt12 < BT_LO + 0.1
3243  || Bt12 > BT_HI - 0.1) {
3244  flags_sst[ip] |= SSTF_BTBAD;
3245  continue;
3246  }
3247 
3248  /* check BT range */
3249  if (Bt11 < Btmin || Bt11 > Btmax || Bt12 < Btmin || Bt12 > Btmax)
3250  flags_sst[ip] |= SSTF_BTRANGE;
3251 
3252  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
3253  if ((ASCEND = avhrr_ascend(btbox, &diflat)) == 1) {
3254  flags_sst[ip] |= SSTF_ASCEND;
3255  }
3256  /* if day and glint > thresh set bit M2B8 (shared with BT4REFDIFF) */
3257  if (l2rec->l1rec->solz[ip] < solznight) {
3258  if (l2rec->l1rec->glint_coef[ip] > glintmax)
3259  flags_sst[ip] |= SSTF_GLINT;
3260  }
3261 
3262  /* ------------------------------------------------------------------------- */
3263  /* Gross Cloud contamination test for tree test BTRANGE (M1B1) */
3264  /* Based on Kilpatrick et al: */
3265  /* Overfiew of NOAA/NASA AVHRR Pathfinder Algorithm. Figure 6. */
3266  /* ------------------------------------------------------------------------- */
3267  /* check BT range: BT must be ge -10 C and le 35 C, M1B1 */
3268  /* some noaa satellites need channel 3, and some only at night */
3269  /* if n7-15 at night or at day and not glint or
3270  * if n16-19 at night */
3271  /* Nov 2012, day only check too cold except older avhrr need too warm check also */
3272 
3273  /* if older avhrr sat and night or day and not glint check for Bt37 too warm or too cold */
3274  if ((l2rec->l1rec->l1file->subsensorID == NO07
3275  || l2rec->l1rec->l1file->subsensorID == NO09
3276  || l2rec->l1rec->l1file->subsensorID == NO10
3277  || l2rec->l1rec->l1file->subsensorID == NO11
3278  || l2rec->l1rec->l1file->subsensorID == NO12
3279  || l2rec->l1rec->l1file->subsensorID == NO14
3280  || l2rec->l1rec->l1file->subsensorID == NO15)
3281  && (l2rec->l1rec->solz[ip] >= solznight
3282  || (l2rec->l1rec->solz[ip] < solznight
3283  && l2rec->l1rec->glint_coef[ip] <= glintmax))) {
3284  if (Bt37 < Btmin || Bt37 > Btmax) {
3285  flags_sst[ip] |= SSTF_BTRANGE;
3286  }
3287  }
3288 
3289  /* if n16-19 Bt37 at night (channel 3 is a different wavelength during the day) is too warm or too cold */
3290  if ((l2rec->l1rec->l1file->subsensorID == NO16
3291  || l2rec->l1rec->l1file->subsensorID == NO17
3292  || l2rec->l1rec->l1file->subsensorID == NO18
3293  || l2rec->l1rec->l1file->subsensorID == NO19)
3294  && l2rec->l1rec->solz[ip] >= solznight) {
3295  if (Bt37 < Btmin || Bt37 > Btmax) {
3296  flags_sst[ip] |= SSTF_BTRANGE;
3297  }
3298  }
3299  } else {
3300  /* for viirs (and modis), to test l2bin dataday splitting and not mixing ascending and descending, set spare flag */
3301  if ((ASCEND = avhrr_ascend(btbox, &diflat)) == 1) {
3302 // l2rec->l1rec->flags[ip] |= ASCFLG; //NEED THIS FLAG DEFINITION FROM SOMEWHERE
3303  }
3304  //if (ip == cpix) {
3305  // printf(" pix=%d line=%d diflat=%f\n",ip,l2rec->l1rec->iscan,diflat);
3306  //}
3307  }
3308  if ((l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA)) {
3309  /* if modis Bt37 or Bt39 or Bt40 at night is too warm or too cold */
3310 
3311  if (l2rec->l1rec->solz[ip] >= solznight
3312  && (Bt37 < Btmin || Bt37 > Btmax || Bt39 < Btmin
3313  || Bt39 > Btmax || Bt40 < Btmin || Bt40 > Btmax40)) {
3314  flags_sst[ip] |= SSTF_BTRANGE;
3315  }
3316  /* if modis Bt37 or Bt39 or Bt40 is during the day too cold */
3317  if ((l2rec->l1rec->solz[ip] < solznight
3318  && l2rec->l1rec->glint_coef[ip] <= glintmax)
3319  && (Bt37 < Btmin || Bt39 < Btmin || Bt40 < Btmin)) {
3320  flags_sst[ip] |= SSTF_BTRANGE;
3321  }
3322  }
3323 
3324  if (l2rec->l1rec->l1file->sensorID == VIIRSN || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
3325  /* if viirs Bt37 or Bt40 at night is too warm or too cold */
3326  if (l2rec->l1rec->solz[ip] >= solznight
3327  && (Bt37 < Btmin || Bt37 > Btmax || Bt40 < Btmin
3328  || Bt40 > Btmax40)) {
3329  flags_sst[ip] |= SSTF_BTRANGE;
3330  }
3331 
3332  /* if viirs Bt37 or Bt40 during the day is too cold */
3333  if ((l2rec->l1rec->solz[ip] < solznight
3334  && l2rec->l1rec->glint_coef[ip] <= glintmax)
3335  && (Bt37 < Btmin || Bt40 < Btmin)) {
3336  flags_sst[ip] |= SSTF_BTRANGE;
3337  }
3338  }
3339 
3340  /* check BT diff */
3341  dBt_11_12 = Bt11 - Bt12;
3342  if (dBt_11_12 < dBtmin || dBt_11_12 > dBtmax)
3343  flags_sst[ip] |= SSTF_BTDIFF;
3344 
3345  /* check SST range, using different max limits for day and night */
3346  if ((sstq[csstbox][ip] < SSTmin)
3347  || (l2rec->l1rec->solz[ip] < solznight && sstq[csstbox][ip] > SSTmax)
3348  || (l2rec->l1rec->solz[ip] >= solznight && sstq[csstbox][ip] > SSTmaxn)) {
3349  flags_sst[ip] |= SSTF_SSTRANGE;
3350  }
3351 
3352  /* check SST difference with references */
3353  /* check for both bad, or don't bother checking anywhere?
3354  * if either are bad then SSTF_SSTRANGE will be set so
3355  * no other flags matter? */
3356  dSST_ref = sstq[csstbox][ip] - l2rec->l1rec->sstref[ip];
3357  // the "coldonly" and equatorial aerosol tests are to be run by default, but if SSTMODS is set, don't
3358  if ((evalmask & SSTMODS) == 0) {
3359  /* evaluate change to cold-test only */
3360  /* set the flag bit if sst is too much colder than reference */
3361  if (dSST_ref < -SSTdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
3362  flags_sst[ip] |= SSTF_SSTREFDIFF;
3363  if (dSST_ref < -input->sstrefdif &&
3364  l2rec->l1rec->lat[ip] >= equatorialSouth && l2rec->l1rec->lat[ip] <= equatorialNorth &&
3365  l2rec->l1rec->lon[ip] >= equatorialWest && l2rec->l1rec->lon[ip] <= equatorialEast) {
3366  /* tighter test between 10S and 30N and -105 to 105 longitude */
3367  /* equatorial aerosol test */
3368  flags_sst[ip] |= SSTF_SSTREFDIFF;
3369  }
3370  /* set the flag bit if sst is too much warmer than reference at night */
3371  if (l2rec->l1rec->solz[ip] >= solznight) {
3372  if (fabs(dSST_ref) > SSTdiff)
3373  flags_sst[ip] |= SSTF_SSTREFDIFF;
3374  }
3375  /* is sst way too much colder than reference? */
3376  if (dSST_ref < -SSTvdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
3377  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3378  /* set the flag bit if sst is way too much warmer than reference at night */
3379  if (l2rec->l1rec->solz[ip] >= solznight) {
3380  if (fabs(dSST_ref) > SSTvdiff)
3381  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3382  }
3383  } else {
3384  if (l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3385  if (fabs(dSST_ref) > SSTdiff)
3386  flags_sst[ip] |= SSTF_SSTREFDIFF;
3387  } else {
3388  if (dSST_ref < -SSTdiff || dSST_ref > (SSTdiff + 1))
3389  flags_sst[ip] |= SSTF_SSTREFDIFF;
3390  }
3391  if (fabs(dSST_ref) > SSTvdiff)
3392  flags_sst[ip] |= SSTF_SSTREFVDIFF;
3393  } //end of to-do-or-not-to-do coldonly block
3394 
3395  /* check SST difference with 4um SST only at night */
3396  /* set BT4REFDIFF flag based on sst4 flags */
3397  if (haveSST4) {
3398  dSST_SST4 = sstq[csstbox][ip] - sst4q[csstbox][ip];
3399  if (sst4q[csstbox][ip] > sstbad + 1.0 && l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3400  if (fabs(dSST_SST4) < SST4diff1)
3401  flags_sst[ip] |= SSTF_SST4DIFF;
3402  if (fabs(dSST_SST4) < SST4diff2)
3403  flags_sst[ip] |= SSTF_SST4VDIFF;
3404  flags_sst[ip] |= (flags_sst4[ip] & SSTF_BT4REFDIFF);
3405  }
3406  }
3407 
3408  /* for viirs, check SST difference with jpss triple window SST only at night */
3409  if (l2rec->l1rec->l1file->sensorID == VIIRSN || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
3410  dSST_SST3 = sstq[csstbox][ip]-sst3q[csstbox][ip];
3411  if (sst3q[csstbox][ip] > sstbad+1.0 && l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
3412  if (fabs(dSST_SST3) > SST3diff1)
3413  flags_sst[ip] |= SSTF_SST3DIFF;
3414  if (fabs(dSST_SST3) > SST3diff2)
3415  flags_sst[ip] |= SSTF_SST3VDIFF;
3416  }
3417  }
3418 
3419  /* check sensor zenith limits */
3420  if (l2rec->l1rec->senz[ip] > hisenz)
3421  flags_sst[ip] |= SSTF_HISENZ;
3422  if (l2rec->l1rec->l1file->sensorID == VIIRSN || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
3423  if (l2rec->l1rec->senz[ip] > vhisenzv2) {
3424  /* different limit only for viirs sst (not sst3) */
3425  flags_sst[ip] |= SSTF_VHISENZ;
3426  }
3427  } else {
3428  /* not viirs */
3429  if (l2rec->l1rec->senz[ip] > vhisenz)
3430  flags_sst[ip] |= SSTF_VHISENZ;
3431  }
3432 
3433  // flag 2 edge pixels as SSTF_VHISENZ so quality gets set to 3
3434  if (l2rec->l1rec->pixnum[ip] < 2 || l2rec->l1rec->pixnum[ip] > (fullscanpix - 3))
3435  flags_sst[ip] |= SSTF_VHISENZ;
3436  // set the last 4 pixels of the scan for Terra to VHISENZ
3437  // as there is an apparent calibration issue with the BT12 for those pixels
3438  if ((l2rec->l1rec->l1file->sensorID == MODIST) && (l2rec->l1rec->pixnum[ip] > 1349))
3439  flags_sst[ip] |= SSTF_VHISENZ;
3440 
3441  /* if SST is cold (collect 5) */
3442  /* I suppose it doesn't hurt to do this for AVHRR since it isn't used to determine quality? */
3443  /* do this for all satellites */
3444  /* sstcloud checks to make sure it's day and not glint */
3445  if (sstq[csstbox][ip] > sstbad + 1.0 && l2rec->l1rec->sstref[ip] > sstbad + 1.0
3446  && sstq[csstbox][ip] - l2rec->l1rec->sstref[ip] <= -1.0)
3447  if (sstcloud(ip, cldbox, cldbox, cldthresh) == 1)
3448  flags_sst[ip] |= SSTF_REDNONUNIF;
3449 
3450  /* check homogeneity of BT */
3451  if (Bt11_maxmin[ip] > Bt11unif1)
3452  flags_sst[ip] |= SSTF_BTNONUNIF;
3453  if (Bt11_maxmin[ip] > Bt11unif2)
3454  flags_sst[ip] |= SSTF_BTVNONUNIF;
3455 
3456  if (Bt12_maxmin[ip] > Bt12unif1)
3457  flags_sst[ip] |= SSTF_BTNONUNIF;
3458  if (Bt12_maxmin[ip] > Bt12unif2)
3459  flags_sst[ip] |= SSTF_BTVNONUNIF;
3460 
3461  /* end of homogeneity checks */
3462 
3463  dBt_37_11 = Bt37 - Bt11;
3464  dBt_37_12 = Bt37 - Bt12;
3465  dBt_11_37 = Bt11 - Bt37;
3466  dBt_67_11 = Bt67 - Bt11;
3467  dBt_40_11 = Bt40 - Bt11;
3468  dBt_85_11 = Bt85 - Bt11;
3469 
3470  /* values from Kay's VIIRS trees */
3471  /* Bt's are BT_LO = -1000 or BT_HI = 1000 when bad */
3472  /* If a Bt is bad then the pixel will be bad so it doesn't matter what happens in the trees */
3473  if (Bt11 != 0.0) {
3474  Tdeflong = (Bt11 - Bt12) / Bt11;
3475  } else {
3476  Tdeflong = (Bt11 - Bt12) / (Bt11 + 0.00001); /* don't divide by zero */
3477  }
3478  double pasutime = l2rec->l1rec->scantime;
3479  int16_t year, day;
3480  double sec;
3481  unix2yds(pasutime, &year, &day, &sec);
3482 
3483  yyyyddd = year * 1000 + day;
3484 
3485  /* --------------------------------------------------------------------------*/
3486  /* decision_tree() - Tree models are based on binary recursive partitioning */
3487  /* to indicate whether data are potentially contaminated. */
3488  /* based on charts sent to NODC from Kay Kilpatrick */
3489  /* worked on by Vicky Lin, SAIC, October 2007. */
3490  /* V6 trees from Guillermo Podesta. */
3491  /* ------------------------------------------------------------------------- */
3492  switch (l2rec->l1rec->l1file->subsensorID) {
3493  case NO07:
3494  /* use the noaa-9 tree for now (Kay: Mar 2013) */
3495  // break;
3496  case NO09:
3497  if (yyyyddd < 1994001) {
3498  /* 9-Jun-97 matchups - NOAA-9 before Jan 1994 */
3499  if ((l2rec->l1rec->solz[ip] >= solznight
3500  || (l2rec->l1rec->solz[ip] < solznight
3501  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3502  && dBt_37_12 < -0.4805) {
3503  if (dBt_37_12 < -1.249) {
3504  flags_sst[ip] |= SSTF_CLOUD;
3505  } else if (dBt_11_12 >= 0.239) {
3506  flags_sst[ip] |= SSTF_CLOUD;
3507  } else if (Bt37_maxmin[ip] >= 0.7975
3508  && dBt_11_12 >= -0.1805) {
3509  flags_sst[ip] |= SSTF_CLOUD;
3510  }
3511  } else {
3512  if (dBt_11_12 < 0.307
3513  && ((l2rec->l1rec->solz[ip] >= solznight
3514  || (l2rec->l1rec->solz[ip] < solznight
3515  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3516  && Bt37_maxmin[ip] >= 1.6215)) {
3517  flags_sst[ip] |= SSTF_CLOUD;
3518  }
3519  }
3520  } else {
3521  /* 9-Jun-97 matchups - NOAA-9 after Jan 1994 */
3522  if ((l2rec->l1rec->solz[ip] >= solznight
3523  || (l2rec->l1rec->solz[ip] < solznight
3524  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3525  && dBt_37_12 < 0.62) {
3526  flags_sst[ip] |= SSTF_CLOUD;
3527  } else if (dBt_11_12 < 0.5055
3528  && ((l2rec->l1rec->solz[ip] >= solznight
3529  || (l2rec->l1rec->solz[ip] < solznight
3530  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3531  && Bt37_maxmin[ip] >= 1.033)
3532  && Bt11_maxmin[ip] >= 0.2435) {
3533  flags_sst[ip] |= SSTF_CLOUD;
3534  }
3535  }
3536  break;
3537  case NO10:
3538  break;
3539  case NO11:
3540  if ((l2rec->l1rec->solz[ip] >= solznight
3541  || (l2rec->l1rec->solz[ip] < solznight && l2rec->l1rec->glint_coef[ip] <= glintmax))
3542  && dBt_37_12 < -0.581) {
3543  flags_sst[ip] |= SSTF_CLOUD;
3544  } else {
3545  if (dBt_11_12 < 0.315) {
3546  if ((l2rec->l1rec->solz[ip] >= solznight
3547  || (l2rec->l1rec->solz[ip] < solznight
3548  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3549  && dBt_11_37 < -0.1855) {
3550  flags_sst[ip] |= SSTF_CLOUD;
3551  } else if ((l2rec->l1rec->solz[ip] >= solznight
3552  || (l2rec->l1rec->solz[ip] < solznight
3553  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3554  && Bt37_maxmin[ip] >= 0.9555) {
3555  flags_sst[ip] |= SSTF_CLOUD;
3556  }
3557  } else if ((l2rec->l1rec->solz[ip] >= solznight
3558  || (l2rec->l1rec->solz[ip] < solznight
3559  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3560  && dBt_11_37 >= 0.86 && Bt37_maxmin[ip] >= 0.5035) {
3561  flags_sst[ip] |= SSTF_CLOUD;
3562  }
3563  }
3564  break;
3565  case NO12:
3566  break;
3567  case NO14:
3568  if (yyyyddd >= 1995001 && yyyyddd < 1996001) {
3569  if ((l2rec->l1rec->solz[ip] >= solznight
3570  || (l2rec->l1rec->solz[ip] < solznight
3571  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3572  && dBt_11_37 >= 0.8015) {
3573  flags_sst[ip] |= SSTF_CLOUD;
3574  } else if ((l2rec->l1rec->solz[ip] >= solznight
3575  || (l2rec->l1rec->solz[ip] < solznight
3576  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3577  && Bt37_maxmin[ip] >= 1.216) {
3578  if (dBt_11_12 < 1.342) {
3579  flags_sst[ip] |= SSTF_CLOUD;
3580  }
3581  } else if (dBt_11_12 < 1.3225
3582  && ((l2rec->l1rec->solz[ip] >= solznight
3583  || (l2rec->l1rec->solz[ip] < solznight
3584  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3585  && dBt_11_37 < 0.994)) {
3586  flags_sst[ip] |= SSTF_CLOUD;
3587  }
3588  } else if (yyyyddd >= 1996001) {
3589  if (dBt_11_12 >= 0.755) {
3590  if ((l2rec->l1rec->solz[ip] >= solznight
3591  || (l2rec->l1rec->solz[ip] < solznight
3592  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3593  && (dBt_11_37 >= 0.7995 || LtNIR8 >= 0.7225)) {
3594  flags_sst[ip] |= SSTF_CLOUD;
3595  }
3596  } else {
3597  if ((l2rec->l1rec->solz[ip] >= solznight
3598  || (l2rec->l1rec->solz[ip] < solznight
3599  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3600  && (dBt_11_37 >= 1.1475 || LtNIR8 >= 1.2645
3601  || Bt37_maxmin[ip] >= 0.966)) {
3602  flags_sst[ip] |= SSTF_CLOUD;
3603  } else {
3604  if (((l2rec->l1rec->solz[ip] >= solznight
3605  || (l2rec->l1rec->solz[ip] < solznight
3606  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3607  && dBt_11_37 >= 0.869)
3608  || (dBt_11_12 < 0.462
3609  && ((l2rec->l1rec->solz[ip] >= solznight
3610  || (l2rec->l1rec->solz[ip] < solznight
3611  && l2rec->l1rec->glint_coef[ip] <= glintmax))
3612  && LtNIR8 >= 0.4865))) {
3613  flags_sst[ip] |= SSTF_CLOUD;
3614  }
3615  }
3616  }
3617  }
3618  break;
3619  case NO15:
3620  break;
3621  case NO16:
3622  /* pathfinder v6 trees */
3623  /* patfinder v6 2001-2003 */
3624  if (l2rec->l1rec->solz[ip] < solznight) {
3625  /* day */
3626  if ((l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7315)
3627  || dBt_11_12 < 0.5375) {
3628  flags_sst[ip] |= SSTF_CLOUD;
3629  }
3630  } else {
3631  /* night */
3632  if (dBt_37_11 < -0.1025 || Bt37_maxmin[ip] >= 1.035
3633  || dBt_11_12 < 0.7385) {
3634  flags_sst[ip] |= SSTF_CLOUD;
3635  } else if (Bt37_maxmin[ip] < 0.7145) {
3636  if (dBt_11_12 < 1.454) {
3637  flags_sst[ip] |= SSTF_CLOUD;
3638  }
3639  } else if (dBt_11_12 < 1.271) {
3640  if (dBt_37_11 >= 0.9195) {
3641  flags_sst[ip] |= SSTF_CLOUD;
3642  }
3643  } else if (dBt_37_11 >= 2.408 && dBt_11_12 < 2.272) {
3644  flags_sst[ip] |= SSTF_CLOUD;
3645  }
3646  }
3647  break;
3648  case NO17:
3649  /* pathfinder v6 trees */
3650  /* NO17 matchups 2003-2006 - pathfinder v6 */
3651  if (l2rec->l1rec->solz[ip] < solznight) {
3652  /* day */
3653  if (l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7045) {
3654  flags_sst[ip] |= SSTF_CLOUD;
3655  } else if (l2rec->l1rec->lat[ip] >= -20.0 && dBt_11_12 < 0.4355) {
3656  flags_sst[ip] |= SSTF_CLOUD;
3657  }
3658  } else {
3659  /* night */
3660  if (dBt_37_12 < 0.0355 || Bt37_maxmin[ip] >= 0.9035
3661  || dBt_11_12 < 0.4565 || dBt_37_11 < -0.5515) {
3662  flags_sst[ip] |= SSTF_CLOUD;
3663  } else if (dBt_37_11 >= 1.635 && dBt_11_12 < 1.693) {
3664  flags_sst[ip] |= SSTF_CLOUD;
3665  }
3666  }
3667  break;
3668  case NO18:
3669  /* pathfinder v6 trees */
3670  if (l2rec->l1rec->solz[ip] < solznight) {
3671  /* day */
3672  /* NOA18 matchups 2005-2009 tree - Pathfinder version 6 */
3673  if (l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.7045) {
3674  flags_sst[ip] |= SSTF_CLOUD;
3675  } else if (dBt_11_12 < 0.4995) {
3676  flags_sst[ip] |= SSTF_CLOUD;
3677  }
3678  } else {
3679  /* night */
3680  /* NOA18 matchups 2005-2009 tree - Pathfinder version 6 */
3681  if (dBt_37_12 < 0.3185) {
3682  flags_sst[ip] |= SSTF_CLOUD;
3683  } else if (Bt37_maxmin[ip] >= 0.8455) {
3684  if (dBt_11_12 < 1.477) {
3685  flags_sst[ip] |= SSTF_CLOUD;
3686  }
3687  } else if (dBt_11_12 < 0.6605) {
3688  flags_sst[ip] |= SSTF_CLOUD;
3689  } else if (dBt_37_11 < -0.5745) {
3690  flags_sst[ip] |= SSTF_CLOUD;
3691  } else if (dBt_37_11 >= 2.1015 && dBt_11_12 < 2.206) {
3692  flags_sst[ip] |= SSTF_CLOUD;
3693  }
3694  }
3695  break;
3696  case NO19:
3697  if (l2rec->l1rec->solz[ip] < solznight) {
3698  /* day */
3699  if ((l2rec->l1rec->glint_coef[ip] <= glintmax && LtRED >= 0.84648)
3700  || dBt_11_12 < 0.168655
3701  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 1.237905) {
3702  flags_sst[ip] |= SSTF_CLOUD;
3703  }
3704  } else {
3705  /* night */
3706  if (dBt_11_12 < 0.99571) {
3707  if (dBt_37_12 < -0.26932 || dBt_11_12 < 0.126555
3708  || Bt37_maxmin[ip] >= 1.02994
3709  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 0.758065) {
3710  flags_sst[ip] |= SSTF_CLOUD;
3711  }
3712  } else if (dBt_37_11 < -0.387815
3713  || fabsf(sstq[csstbox][ip] - l2rec->l1rec->sstref[ip]) >= 0.857825) {
3714  flags_sst[ip] |= SSTF_CLOUD;
3715  }
3716  }
3717  break;
3718  default:
3719  /* modis trees are only v6 */
3720 
3721  /* modis trees were built using gsfc extractions with rho for the vis bands */
3722  /* but l2gen puts radiance in LtRED so need to use rhotRED in the tree tests */
3723 
3724  if (l2rec->l1rec->l1file->sensorID == MODIST) {
3725  if (l2rec->l1rec->solz[ip] < solznight
3726  && l2rec->l1rec->glint_coef[ip] <= glintmax) {
3727  /* day not glint TERRA SST tree test */
3728  treesum[ip] = 0.0;
3729  /* seadas (594,662) is outside glint for comparison to glint pixels that are wrong quality */
3730  /* Bad Lt's could be -32767, 0 (class only so ignore), or 1000. */
3731  /* rho is the same sign so really could check rhotRED instead of LtRED,
3732  except where we have to check for 1000.0 */
3733  if (rhotRED[ip] < 0.206 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3734  treesum[ip] += 0.894;
3735  if (rhotRED[ip] < 0.051 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3736  treesum[ip] += 0.540;
3737  if (l2rec->l1rec->senz[ip] < 64.915) {
3738  treesum[ip] += -0.057;
3739  if (rhoCirrus_max[ip] < 0.003 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3740  treesum[ip] += 0.087;
3741  if (l2rec->l1rec->senz[ip] < 54.529) {
3742  treesum[ip] += 0.030;
3743  } else {
3744  treesum[ip] += -0.573;
3745  }
3746  } else {
3747  treesum[ip] += -0.828;
3748  }
3749  } else {
3750  treesum[ip] += -2.700;
3751  }
3752  } else {
3753  treesum[ip] += -0.490;
3754  }
3755  /* default min is 32767, but gets set to -32767 if no valid values in box */
3756  if (rhoCirrus_min[ip] < 0.004 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3757  treesum[ip] += 0.210;
3758  if (rhotRED[ip] < 0.073 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3759  treesum[ip] += 0.074;
3760  } else {
3761  treesum[ip] += -0.467;
3762  }
3763  } else {
3764  treesum[ip] += -1.008;
3765  }
3766  } else {
3767  treesum[ip] += -1.402;
3768  /* don't need LtRED_min: if Lt is -32767, 0.0 or
3769  1000.0 (which makes rho large) it will go to
3770  the negative treesum */
3771  /* default min is 32767, but gets set to -32767 if no valid values in box */
3772  if (rhotRED_min[ip] < 0.239 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
3773  treesum[ip] += 0.821;
3774  } else {
3775  treesum[ip] += -0.033;
3776  if (dBt_67_11 < -42.423) {
3777  treesum[ip] += 0.283;
3778  } else {
3779  treesum[ip] += -0.184;
3780  }
3781  }
3782  }
3783  if (dBt_37_11 < 3.493) {
3784  treesum[ip] += 0.530;
3785  if (dBt_11_12 < 0.178) {
3786  treesum[ip] += -0.589;
3787  } else {
3788  treesum[ip] += 0.174;
3789  }
3790  } else {
3791  treesum[ip] += -0.169;
3792  /* default min is 32767, but gets set to -32767 if no valid values in box */
3793  if (rhoCirrus_min[ip] < 0.008 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3794  treesum[ip] += 0.097;
3795  } else {
3796  treesum[ip] += -0.761;
3797  }
3798  }
3799  if (sst_stdev[ip] < 0.273 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3800  treesum[ip] += 0.258;
3801  } else {
3802  treesum[ip] += -0.194;
3803  }
3804  if (sst_stdev[ip] < 0.119 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3805  treesum[ip] += 0.459;
3806  } else {
3807  treesum[ip] += -0.037;
3808  if (sstq[csstbox][ip] < 27.293) {
3809  treesum[ip] += -0.054;
3810  } else {
3811  treesum[ip] += 0.411;
3812  }
3813  }
3814  } else if (l2rec->l1rec->solz[ip] < solznight) {
3815  /* day in high and low glint TERRA SST tree tests */
3816  if (LtRED < (1000.0 - 0.1)) {
3817  /* RED is not saturated */
3818  /* day in low glint TERRA SST tree test */
3819  treesum[ip] = 0.0;
3820  /* seadas (600,663) should be qual 1 not 3; terra 20141701425 */
3821  /* seadas (611,698) should be qual 3 not 1; descending so not flipped */
3822  /* don't check for > 0.0 here because we want the negative
3823  treesum if it's < 0.0 or = 1000.0 */
3824  //if (rhotRED[ip] < 0.098 || LtRED >= (1000.0-0.1)) {}
3825  /* Jun 2017, don't know why the above comment and test, but want positive tree sum only if Lt was valid */
3826  if (rhotRED[ip] < 0.098 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3827  treesum[ip] += 1.026;
3828  if (rhotRED[ip] < 0.069 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
3829  treesum[ip] += 0.290;
3830  } else {
3831  treesum[ip] += -0.312;
3832  }
3833  /* rhoCirrus is -32767 when bad */
3834  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
3835  treesum[ip] += 0.154;
3836  } else {
3837  treesum[ip] += -0.574;
3838  }
3839  /* if bad then it will go negative */
3840  if (Bt73_max[ip] < -10.303) {
3841  treesum[ip] += -0.150;
3842  } else {
3843  treesum[ip] += 0.393;
3844  }
3845  } else {
3846  treesum[ip] += -0.761;
3847  if (Bt11 < 16.122) {
3848  treesum[ip] += -0.538;
3849  } else {
3850  treesum[ip] += 0.707;
3851  }
3852  }
3853  /* rhoCirrus is -32767 when bad so make sure it goes negative if bad */
3854  if (rhoCirrus_max[ip] < 0.005 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3855  treesum[ip] += 0.357;
3856  if (Tdeflong < 0.037) {
3857  treesum[ip] += -0.450;
3858  if (l2rec->l1rec->lat[ip] < 28.270) {
3859  treesum[ip] += -0.383;
3860  } else {
3861  treesum[ip] += 0.560;
3862  }
3863  } else {
3864  treesum[ip] += 0.134;
3865  }
3866  if (l2rec->l1rec->senz[ip] < 34.756) {
3867  treesum[ip] += 0.098;
3868  } else {
3869  treesum[ip] += -0.424;
3870  }
3871  } else {
3872  treesum[ip] += -1.399;
3873  /* default min is 32767, but gets set to -32767 if no valid values in box */
3874  if (rhoCirrus_min[ip] < 0.01 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3875  treesum[ip] += 0.326;
3876  } else {
3877  treesum[ip] += -1.404;
3878  }
3879  }
3880  /* sst_stdev is -32767 if bad so make sure it goes negative if bad */
3881  if (sst_stdev[ip] < 0.382 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3882  treesum[ip] += 0.130;
3883  if (l2rec->l1rec->lat[ip] < -5.005) {
3884  treesum[ip] += 0.400;
3885  } else {
3886  treesum[ip] += -0.107;
3887  }
3888  /* sst_stdev is -32767 if bad so make sure it goes negative if bad */
3889  if (sst_stdev[ip] < 0.224 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3890  treesum[ip] += 0.182;
3891  } else {
3892  treesum[ip] += -0.216;
3893  }
3894  } else {
3895  treesum[ip] += -0.632;
3896  }
3897  if (dBt_85_11 < -1.641) {
3898  treesum[ip] += 0.268;
3899  } else {
3900  treesum[ip] += -0.201;
3901  if (sstq[csstbox][ip] < 27.723) {
3902  treesum[ip] += -0.069;
3903  } else {
3904  treesum[ip] += 0.903;
3905  }
3906  }
3907  } else {
3908  /* RED is saturated */
3909  /* day high glint terra sst tree test */
3910  treesum[ip] = 0.0;
3911  /* default min is 32767, but gets set to -32767 if no valid values in box */
3912  if (Bt12_min[ip] < 15.376) {
3913  treesum[ip] += -0.908;
3914  if (dBt_11_12 < 0.27) {
3915  treesum[ip] += -0.446;
3916  } else {
3917  treesum[ip] += 0.568;
3918  }
3919  /* default min is 32767, but gets set to -32767 if no valid values in box */
3920  if (rhoCirrus_min[ip] < -0.001) {
3921  treesum[ip] += -0.981;
3922  } else {
3923  treesum[ip] += 0.137;
3924  }
3925  } else {
3926  treesum[ip] += 0.819;
3927  if (Tdeflong < 0.031) {
3928  treesum[ip] += -0.969;
3929  } else {
3930  treesum[ip] += 0.025;
3931  }
3932  /* bad is -32767 so it will go negative if bad */
3933  if (Bt73_max[ip] < -9.485) {
3934  treesum[ip] += -0.208;
3935  } else {
3936  treesum[ip] += 0.428;
3937  }
3938  }
3939  if (rhoCirrus_max[ip] < 0.005 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
3940  treesum[ip] += 0.235;
3941  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
3942  treesum[ip] += 0.100;
3943  } else {
3944  treesum[ip] += -0.506;
3945  }
3946  } else {
3947  treesum[ip] += -1.176;
3948  }
3949  if (sst_stdev[ip] < 0.36 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3950  treesum[ip] += 0.222;
3951  } else {
3952  treesum[ip] += -0.657;
3953  }
3954  if (sstq[csstbox][ip] < 27.443) {
3955  treesum[ip] += -0.077;
3956  if (sst_stdev[ip] < 0.36 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
3957  treesum[ip] += 0.355;
3958  } else {
3959  treesum[ip] += -0.179;
3960  }
3961  /* default min is 32767, but gets set to -32767 if no valid values in box */
3962  if (rhoCirrus_min[ip] < 0.01 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
3963  treesum[ip] += 0.022;
3964  if (l2rec->l1rec->lat[ip] < 29.775) {
3965  treesum[ip] += -0.128;
3966  if (l2rec->l1rec->lat[ip] < -15.305) {
3967  treesum[ip] += 0.431;
3968  } else {
3969  treesum[ip] += -0.202;
3970  /* bad max is -32767 so it will go negative without a specific test */
3971  if (Bt11_max[ip] < 14.722) {
3972  treesum[ip] += -1.906;
3973  } else {
3974  treesum[ip] += 0.087;
3975  }
3976  }
3977  if (Tdeflong < 0.04) {
3978  treesum[ip] += -0.336;
3979  } else {
3980  treesum[ip] += 0.156;
3981  }
3982  } else {
3983  treesum[ip] += 0.379;
3984  }
3985  } else {
3986  treesum[ip] += -1.948;
3987  }
3988  } else {
3989  treesum[ip] += 0.688;
3990  }
3991  }
3992  } else {
3993  /* night TERRA SST tree test */
3994  treesum[ip] = 0.0;
3995  if (dBt_37_12 < -0.053) {
3996  treesum[ip] += -1.257;
3997  if (l2rec->l1rec->lat[ip] < 33.195) {
3998  treesum[ip] += -0.278;
3999  if (l2rec->l1rec->lat[ip] < -40.185) {
4000  treesum[ip] += 0.619;
4001  } else {
4002  treesum[ip] += -0.711;
4003  if (Bt37 < 6.477) {
4004  treesum[ip] += -3.733;
4005  } else {
4006  treesum[ip] += -0.111;
4007  }
4008  }
4009  } else {
4010  treesum[ip] += 0.333;
4011  }
4012  if (Bt37 < 9.372) {
4013  treesum[ip] += -0.292;
4014  } else {
4015  treesum[ip] += 0.764;
4016  }
4017  } else {
4018  treesum[ip] += 0.430;
4019  if (Bt11_maxmin[ip] < 0.486 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4020  treesum[ip] += 0.628;
4021  if (Bt40_stdev[ip] < 0.146 && Bt40_stdev[ip] > (BAD_FLT + 0.1)) {
4022  treesum[ip] += 0.177;
4023  } else {
4024  treesum[ip] += -0.723;
4025  }
4026  } else {
4027  treesum[ip] += -0.450;
4028  }
4029  if (dSST_SST4 < -0.878) {
4030  treesum[ip] += -1.353;
4031  if (dSST_SST4 < -1.533) {
4032  treesum[ip] += -1.439;
4033  } else {
4034  treesum[ip] += 0.346;
4035  }
4036  } else {
4037  treesum[ip] += 0.219;
4038  if (Bt37_stdev[ip] < 0.448 && Bt37_stdev[ip] > (BAD_FLT + 0.1)) {
4039  treesum[ip] += 0.290;
4040  if (dSST_SST4 < -0.422) {
4041  treesum[ip] += -0.504;
4042  } else {
4043  treesum[ip] += 0.268;
4044  }
4045  } else {
4046  treesum[ip] += -0.484;
4047  }
4048  if (Bt12 < 16.736) {
4049  treesum[ip] += -0.285;
4050  if (dBt_40_11 < -2.199) {
4051  treesum[ip] += 0.518;
4052  } else {
4053  treesum[ip] += -0.316;
4054  if (Bt12 < 11.896) {
4055  treesum[ip] += -0.527;
4056  } else {
4057  treesum[ip] += 0.400;
4058  }
4059  }
4060  } else {
4061  treesum[ip] += 0.500;
4062  }
4063  if (dSST_SST4 < 1.183) {
4064  treesum[ip] += 0.051;
4065  } else {
4066  treesum[ip] += -0.898;
4067  }
4068  }
4069  }
4070  }
4071  if (treesum[ip] <= 0.0) {
4072  flags_sst[ip] |= SSTF_CLOUD;
4073  }
4074  }/* end case for modis terra */
4075  else if (l2rec->l1rec->l1file->sensorID == MODISA) {
4076  if (l2rec->l1rec->solz[ip] < solznight && l2rec->l1rec->glint_coef[ip] <= glintmax) {
4077  /* day not glint AQUA SST tree test */
4078  treesum[ip] = 0.0;
4079  /* seadas (131,1112) should be qual 1 not 3; in l2gen it's (1353-131,2029-1112)=(1222,917) */
4080  /* seadas (1086,1375) should be qual 1 not 3; in l2gen it's (1353-1086,2029-1375)=(1264,654) */
4081  /* seadas says valid rhotRED is >= -0.1 & <= 1.0 */
4082  /* seadas says valid cirrus is >= 0.0 & <= 1.1 */
4083  if (rhotRED[ip] < 0.204 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4084  treesum[ip] += 0.982;
4085  if (rhotRED[ip] < 0.056 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4086  treesum[ip] += 0.328;
4087  /* bad cirrus is -32767 */
4088  if (rhoCirrus_max[ip] < 0.006 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4089  treesum[ip] += 0.151;
4090  if (l2rec->l1rec->senz[ip] < 59.043) {
4091  treesum[ip] += 0.029;
4092  } else {
4093  treesum[ip] += -0.618;
4094  }
4095  } else {
4096  treesum[ip] += -1.071;
4097  }
4098  } else {
4099  treesum[ip] += -0.565;
4100  }
4101  } else {
4102  treesum[ip] += -1.325;
4103  /* if Lt is 1000 then rhot will be large so this will go negative */
4104  /* default min is 32767, but gets set to -32767 if no valid values in box */
4105  if (rhotRED_min[ip] < 0.251 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4106  treesum[ip] += 0.81;
4107  } else {
4108  treesum[ip] += -0.021;
4109  if (l2rec->l1rec->lat[ip] < 48.685) {
4110  treesum[ip] += 0.605;
4111  } else {
4112  treesum[ip] += -0.088;
4113  if (l2rec->l1rec->lat[ip] < 37.755) {
4114  treesum[ip] += -0.311;
4115  } else {
4116  treesum[ip] += 0.189;
4117  }
4118  if (dBt_11_12 < 0.752) {
4119  treesum[ip] += -0.091;
4120  } else {
4121  treesum[ip] += 0.581;
4122  }
4123  }
4124  }
4125  // maybe add check to make sure Bt67 is valid? if ((Bt67 +/- 1000)
4126  if (dBt_67_11 < -38.023) {
4127  treesum[ip] += 0.221;
4128  } else {
4129  treesum[ip] += -0.201;
4130  }
4131  }
4132  /* default min is 32767, but gets set to -32767 if no valid values in box */
4133  if (rhoCirrus_min[ip] < 0.003 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4134  treesum[ip] += 0.22;
4135  if (rhoCirrus_max[ip] < 0.003 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4136  treesum[ip] += 0.069;
4137  } else {
4138  treesum[ip] += -0.422;
4139  }
4140  } else {
4141  treesum[ip] += -0.891;
4142  }
4143  if (dBt_37_11 < 3.577) {
4144  treesum[ip] += 0.719;
4145  } else {
4146  treesum[ip] += -0.227;
4147  if (sstq[csstbox][ip] < 27.668) {
4148  treesum[ip] += -0.017;
4149  } else {
4150  treesum[ip] += 0.472;
4151  }
4152  if (sst_stdev[ip] < 0.481 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4153  treesum[ip] += 0.094;
4154  } else {
4155  treesum[ip] += -0.322;
4156  }
4157  }
4158  if (sst_stdev[ip] < 0.263 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4159  treesum[ip] += 0.193;
4160  } else {
4161  treesum[ip] += -0.276;
4162  }
4163  } else if (l2rec->l1rec->solz[ip] < solznight) {
4164  /* day in high and low glint AQUA SST tree tests */
4165  if (LtRED < (1000.0 - 0.1)) {
4166  /* RED is not saturated */
4167  /* day in low glint AQUA SST tree test */
4168  /* v1 adtree */
4169  /* seadas (403,469) should be qual 1 not 3; in l2gen it's (1353-403,2029-469)=(950,1560) */
4170  treesum[ip] = 0.0;
4171  /* default min is 32767, but gets set to -32767 if no valid values in box */
4172  if (rhotRED_min[ip] < 0.087 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4173  treesum[ip] += 1.017;
4174  } else {
4175  treesum[ip] += -0.644;
4176  /* default min is 32767, but gets set to -32767 if no valid values in box */
4177  if (Bt11_min[ip] < 15.22) {
4178  treesum[ip] += -0.837;
4179  if (dBt_67_11 < -40.834) {
4180  treesum[ip] += 0.510;
4181  } else {
4182  treesum[ip] += -0.181;
4183  }
4184  } else {
4185  treesum[ip] += 0.870;
4186  if (Tdeflong < 0.036) {
4187  treesum[ip] += -0.859;
4188  } else {
4189  treesum[ip] += 0.215;
4190  }
4191  }
4192  if (l2rec->l1rec->lat[ip] < 29.675) {
4193  treesum[ip] += -0.099;
4194  } else {
4195  treesum[ip] += 0.467;
4196  }
4197  if (dBt_37_11 < 10.053) {
4198  treesum[ip] += -0.743;
4199  } else {
4200  treesum[ip] += 0.137;
4201  }
4202  }
4203  if (rhoCirrus_max[ip] < 0.004 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4204  treesum[ip] += 0.227;
4205  if (rhotRED[ip] < 0.068 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4206  treesum[ip] += 0.824;
4207  } else {
4208  treesum[ip] += -0.113;
4209  }
4210  if (rhoCirrus_maxmin[ip] < 0.001 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4211  treesum[ip] += 0.102;
4212  } else {
4213  treesum[ip] += -0.454;
4214  }
4215  } else {
4216  treesum[ip] += -1.096;
4217  if (rhoCirrus < 0.006 && rhoCirrus > (BAD_FLT + 0.1)) {
4218  treesum[ip] += 0.282;
4219  } else {
4220  treesum[ip] += -0.810;
4221  }
4222  }
4223  if (sst_stdev[ip] < 0.397 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4224  treesum[ip] += 0.235;
4225  if (l2rec->l1rec->lat[ip] < -0.365) {
4226  treesum[ip] += 0.240;
4227  } else {
4228  treesum[ip] += -0.194;
4229  if (l2rec->l1rec->lon[ip] < -67.480) {
4230  treesum[ip] += 0.241;
4231  } else {
4232  treesum[ip] += -0.229;
4233  }
4234  }
4235  if (sst_stdev[ip] < 0.188 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4236  treesum[ip] += 0.187;
4237  } else {
4238  treesum[ip] += -0.258;
4239  }
4240  } else {
4241  treesum[ip] += -0.848;
4242  }
4243  if (sstq[csstbox][ip] < 27.278) {
4244  treesum[ip] += -0.066;
4245  } else {
4246  treesum[ip] += 0.509;
4247  }
4248  // /* v2 adtree, might not be needed, try v1 first */
4249  // /* day in glint AQUA SST tree test */
4250  // treesum[ip] = 0.0;
4251  // if (rhotRED_min[ip] < 15.119 && rhotRED_min[ip] > 0.0) {
4252  // treesum[ip] += -0.759;
4253  // } else {
4254  // treesum[ip] += 0.794;
4255  // }
4256  // if (sst_stdev[ip] < 0.254) {
4257  // treesum[ip] += 0.304;
4258  // if (sst_stdev[ip] < 0.14) {
4259  // treesum[ip] += 0.237;
4260  // } else {
4261  // treesum[ip] += -0.276;
4262  // }
4263  // if (l2rec->l1rec->lat[ip] < -4.225){
4264  // treesum[ip] += 0.401;
4265  // } else {
4266  // treesum[ip] += -0.144;
4267  // }
4268  // } else {
4269  // treesum[ip] += -0.508;
4270  // }
4271  // if (dBt_11_12 < 0.303) {
4272  // treesum[ip] += -0.715;
4273  // } else {
4274  // treesum[ip] += 0.184;
4275  // if (l2rec->l1rec->lat[ip] < -23.905) {
4276  // treesum[ip] += 0.662;
4277  // } else {
4278  // treesum[ip] += -0.087;
4279  // }
4280  // if (sst_stdev[ip] < 0.494) {
4281  // treesum[ip] += 0.114;
4282  // } else {
4283  // treesum[ip] += -0.550;
4284  // }
4285  // }
4286  // if (l2rec->l1rec->lat[ip] < 32.415) {
4287  // treesum[ip] += -0.137;
4288  // if (sstq[csstbox][ip] < 26.408) {
4289  // treesum[ip] += -0.245;
4290  // } else {
4291  // treesum[ip] += 0.499;
4292  // }
4293  // if (Tdeflong < 0.038) {
4294  // treesum[ip] += -0.415;
4295  // } else {
4296  // treesum[ip] += 0.038;
4297  // }
4298  // } else {
4299  // treesum[ip] += 0.506;
4300  // }
4301  } else {
4302  /* RED is saturated */
4303  /* day high glint aqua sst tree test */
4304  treesum[ip] = 0.0;
4305  /* default min is 32767, but gets set to -32767 if no valid values in box */
4306  /* if min is bad this will go negative without a specific test */
4307  if (Bt12_min[ip] < 15.119) {
4308  treesum[ip] += -0.759;
4309  if (dBt_37_11 < 12.733) {
4310  treesum[ip] += 0.652;
4311  /* default min is 32767, but gets set to -32767 if no valid values in box */
4312  /* if min is bad this will go negative without a specific test */
4313  if (rhoCirrus_min[ip] < 0.001) {
4314  treesum[ip] += -1.096;
4315  } else {
4316  treesum[ip] += 0.130;
4317  }
4318  } else {
4319  treesum[ip] += -0.288;
4320  }
4321  if (Bt37_maxmin[ip] < 0.529 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4322  treesum[ip] += 0.926;
4323  } else {
4324  treesum[ip] += 0.026;
4325  }
4326  } else {
4327  treesum[ip] += 0.794;
4328  if (dBt_37_11 < 8.37) {
4329  treesum[ip] += 0.522;
4330  } else {
4331  treesum[ip] += -0.268;
4332  }
4333  }
4334  if (rhoCirrus_maxmin[ip] < 0.001 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4335  treesum[ip] += 0.381;
4336  if (rhoCirrus_max[ip] < 0.002 && rhoCirrus_max[ip] > (BAD_FLT + 0.1)) {
4337  treesum[ip] += 0.006;
4338  } else {
4339  treesum[ip] += -0.584;
4340  }
4341  } else {
4342  treesum[ip] += -0.619;
4343  }
4344  if (sst_stdev[ip] < 0.397 && sst_stdev[ip] > (BAD_FLT + 0.1)) {
4345  treesum[ip] += 0.101;
4346  } else {
4347  treesum[ip] += -0.684;
4348  }
4349  if (Tdeflong < 0.036) {
4350  treesum[ip] += -0.471;
4351  if (l2rec->l1rec->lat[ip] < 32.335) {
4352  treesum[ip] += -0.481;
4353  } else {
4354  treesum[ip] += 0.684;
4355  }
4356  } else {
4357  treesum[ip] += 0.193;
4358  }
4359  /* default min is 32767, but gets set to -32767 if no valid values in box */
4360  if (rhoCirrus_min[ip] < 0.003 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4361  treesum[ip] += 0.108;
4362  if (l2rec->l1rec->lat[ip] < -25.425) {
4363  treesum[ip] += 0.517;
4364  } else {
4365  treesum[ip] += -0.082;
4366  if (l2rec->l1rec->lat[ip] < 27.215) {
4367  treesum[ip] += -0.077;
4368  if (Bt11 < 14.734) {
4369  treesum[ip] += -2.012;
4370  } else {
4371  treesum[ip] += 0.089;
4372  }
4373  } else {
4374  treesum[ip] += 0.383;
4375  }
4376  }
4377  } else {
4378  treesum[ip] += -0.913;
4379  }
4380  if (sstq[csstbox][ip] < 27.583) {
4381  treesum[ip] += -0.074;
4382  } else {
4383  treesum[ip] += 0.691;
4384  }
4385  }
4386  } else {
4387  /* night AQUA SST tree test */
4388  /* use tree that starts with 0, not -0.029 */
4389  treesum[ip] = 0.0;
4390  if (dBt_37_12 < 0.117) {
4391  treesum[ip] += -1.279;
4392  if (l2rec->l1rec->lat[ip] < 33.035) {
4393  treesum[ip] += -0.289;
4394  if (l2rec->l1rec->lat[ip] < -42.235) {
4395  treesum[ip] += 0.747;
4396  } else {
4397  treesum[ip] += -0.564;
4398  if (Bt37 < 13.117) {
4399  treesum[ip] += -0.580;
4400  } else {
4401  treesum[ip] += 0.638;
4402  }
4403  }
4404  } else {
4405  treesum[ip] += 0.355;
4406  }
4407  if (Bt37 < 9.447) {
4408  treesum[ip] += -0.307;
4409  } else {
4410  treesum[ip] += 0.747;
4411  }
4412  } else {
4413  treesum[ip] += 0.470;
4414  if (Bt11_stdev[ip] < 0.155 && Bt11_stdev[ip] > (BAD_FLT + 0.1)) {
4415  treesum[ip] += 0.690;
4416  if (Bt37_maxmin[ip] < 0.524 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4417  treesum[ip] += 0.150;
4418  } else {
4419  treesum[ip] += -0.794;
4420  }
4421  } else {
4422  treesum[ip] += -0.430;
4423  }
4424  if (dSST_SST4 < -0.787) {
4425  treesum[ip] += -1.404;
4426  if (dSST_SST4 < -1.253) {
4427  treesum[ip] += -1.086;
4428  } else {
4429  treesum[ip] += 0.388;
4430  }
4431  } else {
4432  treesum[ip] += 0.197;
4433  if (Bt37_maxmin[ip] < 1.229 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
4434  treesum[ip] += 0.287;
4435  if (dSST_SST4 < -0.383) {
4436  treesum[ip] += -0.531;
4437  } else {
4438  treesum[ip] += 0.279;
4439  }
4440  } else {
4441  treesum[ip] += -0.497;
4442  }
4443  if (Bt12 < 16.353) {
4444  treesum[ip] += -0.300;
4445  if (dBt_40_11 < -2.171) {
4446  treesum[ip] += 0.516;
4447  } else {
4448  treesum[ip] += -0.395;
4449  if (sst4q[csstbox][ip] < 16.212) {
4450  treesum[ip] += -0.694;
4451  } else {
4452  treesum[ip] += 0.392;
4453  }
4454  }
4455  } else {
4456  treesum[ip] += 0.451;
4457  }
4458  }
4459  // maybe add check to make sure Bt85 is valid? if ((Bt85 +/- 1000)
4460  if (dBt_85_11 < -1.577) {
4461  treesum[ip] += 0.088;
4462  } else {
4463  treesum[ip] += -0.408;
4464  }
4465  }
4466  }
4467  if (treesum[ip] <= 0.0) {
4468  flags_sst[ip] |= SSTF_CLOUD;
4469  }
4470  }/* end case for modis */
4471 
4472  else if (l2rec->l1rec->l1file->sensorID == VIIRSN || l2rec->l1rec->l1file->sensorID == VIIRSJ1) {
4473  /* if we run this more than once I'll add a switch, but I think it's just here to make a 'previous' picture */
4474  v5viirs = 0;
4475  /* v5 viirs */
4476  if (v5viirs == 1) {
4477  if (l2rec->l1rec->solz[ip] < solznight) {
4478  /* day VIIRS SST tree test */
4479  /* v5.33 day tree */
4480 
4481  /* sd6.4 version used class files that contained values in K */
4482  /* sd7.4 version uses gsfc files with values in C */
4483  if (dBt_11_12 < 0.23435 ||
4484  Bt12 < (270.78 - CtoK)) {
4485  flags_sst[ip] |= SSTF_CLOUD;
4486  } else if (Bt37_maxmin[ip] >= 0.40919 &&
4487  (Bt37_maxmin[ip] >= 1.0765 ||
4488  Bt12_maxmin[ip] >= 0.8412 ||
4489  Bt12_maxmin[ip] < 0.3193)) {
4490  flags_sst[ip] |= SSTF_CLOUD;
4491  }
4492  } else {
4493  /* night VIIRS SST tree test */
4494  /* v5.33 night tree */
4495  if (dBt_37_11 < 0.18654) {
4496  flags_sst[ip] |= SSTF_CLOUD;
4497  } else if (dSST_SST3 < -1.021674) {
4498  if (dBt_11_12 < 2.187215 || dBt_37_11 >= 5.46243) {
4499  flags_sst[ip] |= SSTF_CLOUD;
4500  }
4501  } else if (Bt37_maxmin[ip] >= 0.48349) {
4502  if (dBt_11_12 < 0.32251 || dSST_SST3 >= 0.4414579) {
4503  flags_sst[ip] |= SSTF_CLOUD;
4504  } else if (Bt37_maxmin[ip] >= 1.012265) {
4505  if (dBt_11_12 < 0.951295 || Bt12_maxmin[ip] < 0.69827) {
4506  flags_sst[ip] |= SSTF_CLOUD;
4507  }
4508  } else {
4509  if (l2rec->l1rec->senz[ip] > 56.3) {
4510  if (dBt_37_11 < 1.101845 || dBt_11_12 >= 1.308475) {
4511  flags_sst[ip] |= SSTF_CLOUD;
4512  }
4513  }
4514  }
4515  }
4516  }
4517  } else {
4518 
4519  /* v6 viirs tree */
4520  if (l2rec->l1rec->solz[ip] >= solznight) {
4521  /* night VIIRS v6 SST tree test */
4522  treesum[ip] = 0.413;
4523  if (dBt_37_11 < 0.115) {
4524  treesum[ip] += -1.931;
4525  } else {
4526  treesum[ip] += 0.496;
4527  if (dSST_SST3 < -0.465) {
4528  treesum[ip] += -0.942;
4529  if (dSST_SST3 < -1.102) {
4530  treesum[ip] += -0.735;
4531  if (dBt_37_11 < 1.686) {
4532  treesum[ip] += 0.651;
4533  } else {
4534  treesum[ip] += -0.82;
4535  }
4536  } else {
4537  treesum[ip] += 0.207;
4538  }
4539  } else {
4540  treesum[ip] += 0.309;
4541  if (dSST_SST3 < 0.78) {
4542  treesum[ip] += 0.043;
4543  if (dSST_SST3 < -0.158) {
4544  treesum[ip] += -0.392;
4545  } else {
4546  treesum[ip] += 0.194;
4547  }
4548  } else {
4549  treesum[ip] += -1.069;
4550  }
4551  if (dBt_37_11 < 0.306) {
4552  treesum[ip] += -0.692;
4553  } else {
4554  treesum[ip] += 0.077;
4555  }
4556  }
4557  if (Bt12 < 287.407) {
4558  treesum[ip] += -0.339;
4559  } else {
4560  treesum[ip] += 0.426;
4561  }
4562  }
4563  if (sstq[csstbox][ip] < (270.414 - CtoK)) {
4564  treesum[ip] += -3.879;
4565  } else {
4566  treesum[ip] += 0.047;
4567  if (Bt37_stdev[ip] < 0.247 && Bt37_stdev[ip] > (BAD_FLT + 0.1)) {
4568  treesum[ip] += 0.176;
4569  } else {
4570  treesum[ip] += -0.187;
4571  }
4572  }
4573  if (treesum[ip] <= 0.0) {
4574  /* set bit to show failed tree test */
4575  flags_sst[ip] |= SSTF_CLOUD;
4576  }
4577  /* end case for v6 viirs night */
4578  } else if (l2rec->l1rec->glint_coef[ip] <= glintmax) {
4579  /* day not glint VIIRS v6 SST tree test */
4580  treesum[ip] = 0.0;
4581 
4582  /* add ice test here for now. May want to add it to set_qual instead of setting cloud flag */
4583  /* only for viirs, day, not glint */
4584  // check for bad rhotred? if (rhotRED[ip] > 0.5 && LtRED > (BAD_FLT+0.1) && LtRED < (1000.0-0.1) &&
4585 
4586  /* instead of lat > 40, calc sub solar point and +/- 30 degrees to get lat bounds that vary by season */
4587  /* units of x=radians=(deg/day)*day/(deg/rad)
4588  * subsolar is negative in the southern hemisphere
4589  */
4590  xdoy = day + 284.0;
4591  xrad = (360.0 / 365.0) * xdoy / RADEG;
4592  subsolar = 23.45 * sin(xrad);
4593 
4594  if (rhotRED[ip] > 0.3 &&
4595  rhot16[ip] >= 0.006 && rhot16[ip] < 0.1 &&
4596  ((l2rec->l1rec->lat[ip] > (subsolar + 30.0)) || (l2rec->l1rec->lat[ip] < (subsolar - 30.0)))) {
4597  /* set cloud bit to show ice (or wispy clouds?) */
4598  flags_sst[ip] |= SSTF_CLOUD;
4599  /* don't really need to do the tree test now, but it doesn't hurt */
4600  }
4601 
4602  /* rho's are bad when Lt's are -32767 or 0.0 (make rho <=0) or 1000 (make rho large) */
4603  // 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)
4604  /* seadas (1907,3148) should be qual 0 not 3; in l2gen it's (3199-1907,3231-3148)=(1292,83) */
4605  /* seadas (1935,3147) should be qual 0 not 0; in l2gen it's (3199-1935,3231-3147)=(1264,84) */
4606  if (rhot16[ip] < 0.16 && Lt16 > (BAD_FLT + 0.1) && Lt16 < (1000.0 - 0.1)) {
4607  treesum[ip] += 0.805;
4608  if (rhotNIR7[ip] < 0.062 && LtNIR7 > (BAD_FLT + 0.1) && LtNIR7 < (1000.0 - 0.1)) {
4609  treesum[ip] += 0.393;
4610  if (rhoCirrus < 0.004 && rhoCirrus > (BAD_FLT + 0.1)) {
4611  treesum[ip] += 0.287;
4612  if (Tdeflong < 0.002) {
4613  treesum[ip] += -0.681;
4614  } else {
4615  treesum[ip] += 0.026;
4616  if (rhotNIR7[ip] < 0.039 && LtNIR7 > (BAD_FLT + 0.1) && LtNIR7 < (1000.0 - 0.1)) {
4617  treesum[ip] += 0.364;
4618  } else {
4619  treesum[ip] += -0.21;
4620  }
4621  }
4622  } else {
4623  treesum[ip] += -1.244;
4624  }
4625  } else {
4626  treesum[ip] += -0.0572;
4627  /* default min is 32767, but gets set to -32767 if no valid values in box */
4628  if (rhot16_min[ip] < 0.032 && rhot16_min[ip] > (BAD_FLT + 0.1)) {
4629  treesum[ip] += 0.455;
4630  } else {
4631  treesum[ip] += -0.395;
4632  }
4633  }
4634  if (l2rec->l1rec->senz[ip] < 64.994) {
4635  treesum[ip] += 0.216;
4636  if (rhoCirrus < 0.007 && rhoCirrus > (BAD_FLT + 0.1)) {
4637  treesum[ip] += 0.065;
4638  } else {
4639  treesum[ip] += -1.077;
4640  }
4641  } else {
4642  treesum[ip] += -0.708;
4643  }
4644  } else {
4645  treesum[ip] += -1.755;
4646  if (rhot16[ip] < 0.266 && Lt16 > (BAD_FLT + 0.1) && Lt16 < (1000.0 - 0.1)) {
4647  treesum[ip] += 0.642;
4648  } else {
4649  treesum[ip] += -0.19;
4650  if (rhotRED_maxmin[ip] < 0.103 && rhotRED_maxmin[ip] > (BAD_FLT + 0.1)) {
4651  treesum[ip] += 0.425;
4652  } else {
4653  treesum[ip] += -0.195;
4654  }
4655  }
4656  if (dBt_11_12 < 0.235) {
4657  treesum[ip] += -0.189;
4658  } else {
4659  treesum[ip] += 0.411;
4660  }
4661  if (l2rec->l1rec->wv[ip] < 2.946) {
4662  treesum[ip] += 0.038;
4663  } else {
4664  treesum[ip] += -1.137;
4665  }
4666  }
4667  if (Bt11_maxmin[ip] < 0.762 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4668  treesum[ip] += 0.156;
4669  } else {
4670  treesum[ip] += -0.188;
4671  }
4672  if (l2rec->l1rec->wv[ip] < 1.315) {
4673  treesum[ip] += 0.327;
4674  } else {
4675  treesum[ip] += -0.054;
4676  if (sstq[csstbox][ip] < (278.171 - CtoK)) {
4677  treesum[ip] += -0.679;
4678  } else {
4679  treesum[ip] += 0.05;
4680  }
4681  }
4682  } else {
4683  /* day glint VIIRS v6 SST tree test */
4684  treesum[ip] = 0.0;
4685  /* seadas (1907,3148) should be qual 0 not 3; in l2gen it's (3199-1907,3231-3148)=(1292,83) */
4686  /* seadas (1935,3147) should be qual 0 not 0; in l2gen it's (3199-1935,3231-3147)=(1264,84) */
4687  if (rhotRED[ip] > 0.065) {
4688  /* high glint tree */
4689  treesum[ip] = 0.0;
4690  /* default min is 32767, but gets set to -32767 if no valid values in box */
4691  if (Bt85_min[ip] < (287.451 - CtoK)) {
4692  treesum[ip] += -0.812;
4693  if (l2rec->l1rec->lat[ip] < 32.315) {
4694  treesum[ip] += -0.296;
4695  if (Tdeflong < 0.001) {
4696  treesum[ip] += -1.109;
4697  if (l2rec->l1rec->lat[ip] < -30.0) {
4698  treesum[ip] += 0.525;
4699  } else {
4700  treesum[ip] += -1.827;
4701  }
4702  } else {
4703  treesum[ip] += 0.44;
4704  }
4705  if (l2rec->l1rec->wv[ip] < 2.065) {
4706  treesum[ip] += 0.49;
4707  } else {
4708  treesum[ip] += -1.104;
4709  if (Bt85 < (287.059 - CtoK)) {
4710  treesum[ip] += -1.071;
4711  } else {
4712  treesum[ip] += 0.727;
4713  }
4714  }
4715  } else {
4716  treesum[ip] += 0.669;
4717  }
4718  if (dBt_37_11 < 17.594) {
4719  treesum[ip] += 0.452;
4720  } else {
4721  treesum[ip] += -0.76;
4722  }
4723  } else {
4724  treesum[ip] += 0.858;
4725  if (l2rec->l1rec->lon[ip] < -69.475) {
4726  treesum[ip] += 0.512;
4727  } else {
4728  treesum[ip] += -0.176;
4729  if (l2rec->l1rec->lat[ip] < 1.01) {
4730  treesum[ip] += 0.64;
4731  } else {
4732  treesum[ip] += -0.176;
4733  if (tmonth < 6) {
4734  treesum[ip] += 0.69;
4735  } else {
4736  treesum[ip] += -0.305;
4737  }
4738  }
4739  }
4740  if (dBt_37_11 < 9.655) {
4741  treesum[ip] += 0.562;
4742  } else {
4743  treesum[ip] += -0.173;
4744  }
4745  if (Bt85 < (290.343 - CtoK)) {
4746  treesum[ip] += -0.39;
4747  } else {
4748  treesum[ip] += 0.243;
4749  if (Tdeflong < 0.003) {
4750  treesum[ip] += -0.817;
4751  } else {
4752  treesum[ip] += 0.267;
4753  }
4754  }
4755  if (l2rec->l1rec->lat[ip] < 22.465) {
4756  treesum[ip] += -0.206;
4757  } else {
4758  treesum[ip] += 0.523;
4759  }
4760  }
4761  if (Bt11_maxmin[ip] < 1.189 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
4762  treesum[ip] += 0.059;
4763  } else {
4764  treesum[ip] += -1.22;
4765  }
4766  } else {
4767  /* low glint tree */
4768  treesum[ip] = 0.0;
4769  /* default min is 32767, but gets set to -32767 if no valid values in box */
4770  if (rhotNIR7_min[ip] < 0.104 && rhotNIR7_min[ip] > (BAD_FLT + 0.1)) {
4771  treesum[ip] += 0.91;
4772  if (rhotRED[ip] < 0.086 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4773  treesum[ip] += 0.518;
4774  /* default min is 32767, but gets set to -32767 if no valid values in box */
4775  if (rhotRED_min[ip] < 0.067 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4776  treesum[ip] += 0.558;
4777  } else {
4778  treesum[ip] += -0.263;
4779  /* don't check for > BAD_FLT here because we want the negative treesum if Lt < 0.0 (or =1000.0) */
4780  if (rhot16[ip] < 0.06 || Lt16 > (1000.0 - 0.1)) {
4781  treesum[ip] += -0.231;
4782  } else {
4783  treesum[ip] += 1.712;
4784  }
4785  /* don't check for > BAD_FLT here because we want the negative treesum if Lt < 0.0 (or =1000.0) */
4786  if (rhot16[ip] < 0.046 || Lt16 > (1000.0 - 0.1)) {
4787  treesum[ip] += -1.353;
4788  } else {
4789  treesum[ip] += 0.352;
4790  }
4791  }
4792  } else {
4793  treesum[ip] += -0.585;
4794  if (dBt_37_12 < 12.951) {
4795  treesum[ip] += -0.905;
4796  } else {
4797  treesum[ip] += 0.187;
4798  if (rhotRED[ip] < 0.098 && LtRED > (BAD_FLT + 0.1) && LtRED < (1000.0 - 0.1)) {
4799  treesum[ip] += 0.549;
4800  } else {
4801  treesum[ip] += -0.484;
4802  }
4803  }
4804  }
4805  } else {
4806  treesum[ip] += -1.819;
4807  /* default min is 32767, but gets set to -32767 if no valid values in box */
4808  if (rhotRED_min[ip] < 0.206 && rhotRED_min[ip] > (BAD_FLT + 0.1)) {
4809  treesum[ip] += 0.467;
4810  } else {
4811  treesum[ip] += -1.18;
4812  if (rhotRED_maxmin[ip] < 0.037 && rhotRED_maxmin[ip] > (BAD_FLT + 0.1)) {
4813  treesum[ip] += 1.747;
4814  } else {
4815  treesum[ip] += -1.79;
4816  }
4817  }
4818  if (l2rec->l1rec->wv[ip] < 1.705) {
4819  treesum[ip] += 0.434;
4820  } else {
4821  treesum[ip] += -0.645;
4822  }
4823  }
4824  if (rhoCirrus < 0.002 && rhoCirrus > (BAD_FLT + 0.1)) {
4825  treesum[ip] += 0.23;
4826  if (l2rec->l1rec->lat[ip] < 32.13) {
4827  treesum[ip] += -0.067;
4828  if (sstq[csstbox][ip] < (300.034 - CtoK)) {
4829  treesum[ip] += -0.146;
4830  } else {
4831  treesum[ip] += 0.824;
4832  }
4833  } else {
4834  treesum[ip] += 0.758;
4835  }
4836  } else {
4837  treesum[ip] += -1.153;
4838  /* default min is 32767, but gets set to -32767 if no valid values in box */
4839  if (rhoCirrus_min[ip] < 0.005 && rhoCirrus_min[ip] > (BAD_FLT + 0.1)) {
4840  treesum[ip] += 0.28;
4841  } else {
4842  treesum[ip] += -0.939;
4843  }
4844  }
4845  if (l2rec->l1rec->senz[ip] < 33.204) {
4846  treesum[ip] += 0.2;
4847  } else {
4848  treesum[ip] += -0.331;
4849  }
4850  }
4851  }
4852  if (treesum[ip] <= 0.0) {
4853  flags_sst[ip] |= SSTF_CLOUD;
4854  }
4855 
4856  } /* end case for viirs v5 vs v6 */
4857 
4858  }/* end case for viirs */
4859  break;
4860  } /* end switch statement for tree tests */
4861 
4862  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
4863  /* Check stray sunlight test, M2B2 (shared bit with SST4DIFF) */
4864  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
4865 
4866  if (l2rec->l1rec->lat[ip] < 0.0
4867  && ((ASCEND == 1 && l2rec->l1rec->pixnum[ip] <= cpix)
4868  || (ASCEND == 0 && l2rec->l1rec->pixnum[ip] > cpix))
4869  && l2rec->l1rec->senz[ip] > 45.00) {
4870  flags_sst[ip] |= SSTF_SUNLIGHT;
4871  }
4872  }
4873 
4874  } /* End of pixel loop */
4875  return;
4876 }
4877 
4878 /* ----------------------------------------------------------------------------------- */
4879 /* set_flags_sst4() - set quality flags for short wave sea surface temperature */
4880 /* */
4881 /* B. Franz, SAIC, August 2005. */
4882 
4883 /* ----------------------------------------------------------------------------------- */
4884 void set_flags_sst4(l2str *l2rec) {
4885 
4886  extern l1qstr l1que;
4887  int32_t npix = l2rec->l1rec->npix;
4888  int32_t ip, ipbir;
4889  float Bt37, Bt39, Bt40, Bt85;
4890  float Bt11, Bt12;
4891  float dBt_37_12, dBt_40_11, dBt_85_11;
4892  float dBt_34;
4893  float dSST4_ref; /* sst4 - ref */
4894  float dSST_SST4; /* sst - sst4 */
4895  statstr statrec;
4896 
4897  /* check each pixel in scan */
4898  for (ip = 0; ip < npix; ip++) {
4899 
4900  d3940ref[ip] = sstbad;
4901 
4902  /* SST not processed */
4903  if (sstmasked(l2rec->l1rec->flags, ip)) {
4904  flags_sst4[ip] |= SSTF_ISMASKED;
4905  continue;
4906  }
4907 
4908  ipbir = ip * NBANDSIR;
4909 
4910  Bt37 = l2rec->l1rec->Bt[ipbir + ib37]; /* modis chan 20, avhrr cen3, viirs m12 */
4911  Bt39 = l2rec->l1rec->Bt[ipbir + ib39]; /* modis chan 22 */
4912  Bt40 = l2rec->l1rec->Bt[ipbir + ib40]; /* modis chan 23, viirs m13 */
4913  Bt85 = l2rec->l1rec->Bt[ipbir + ib85]; /* modis chan 29, viirs m14 */
4914  Bt11 = l2rec->l1rec->Bt[ipbir + ib11]; /* modis chan 31, avhrr cen4, viirs m15 */
4915  Bt12 = l2rec->l1rec->Bt[ipbir + ib12]; /* modis chan 32, avhrr cen5, viirs m16 */
4916 
4917  /* aqua detector 0 channel 23 has a problem so average detectors 9 and 1 instead */
4918  if ((l2rec->l1rec->l1file->sensorID == MODISA)
4919  && l2rec->l1rec->detnum == 0) {
4920  /* current row is center of box */
4921  int32_t is = l1que.nq / 2;
4922  if (btavg(ip, is, ib40, NBANDSIR, &statrec) > 0) {
4923  Bt40 = statrec.avg;
4924  }
4925  }
4926  dBt_37_12 = Bt37 - Bt12;
4927  dBt_40_11 = Bt40 - Bt11;
4928  dBt_85_11 = Bt85 - Bt11;
4929 
4930  /* if BT could not be computed (radiance out-of-range) */
4931  if (Bt39 < BT_LO + 0.1 || Bt39 > BT_HI - 0.1 || Bt40 < BT_LO + 0.1
4932  || Bt40 > BT_HI - 0.1) {
4933  flags_sst4[ip] |= SSTF_BTBAD;
4934  continue;
4935  }
4936 
4937  /* check BT range (don't care about Bt37 for sst4) */
4938  if (Bt39 < Btmin || Bt39 > Btmax || Bt40 < Btmin || Bt40 > Btmax40)
4939  flags_sst4[ip] |= SSTF_BTRANGE;
4940 
4941  /* check BT diff */
4942  dBt_34 = Bt39 - Bt40;
4943  if (dBt_34 < dBt4min || dBt_34 > dBt4max)
4944  flags_sst4[ip] |= SSTF_BTDIFF;
4945 
4946  /* check BT diff against BT reference */
4947 
4948  /* pathfinder v6 */
4949  d3940ref[ip] = btrefdiffv6(ip, Bt39, Bt40, l2rec);
4950 
4951  if (d3940ref[ip] < dBtrefmin || d3940ref[ip] > dBtrefmax)
4952  flags_sst4[ip] |= SSTF_BT4REFDIFF;
4953 
4954  /* check SST range */
4955  if (sst4q[csstbox][ip] < SSTmin || sst4q[csstbox][ip] > SSTmaxn)
4956  flags_sst4[ip] |= SSTF_SSTRANGE;
4957 
4958  /* check SST difference with reference */
4959  if (sst4q[csstbox][ip] < sstbad + 1.0 || l2rec->l1rec->sstref[ip] < sstbad + 1.0) {
4960  dSST4_ref = sstbad;
4961  } else {
4962  dSST4_ref = sst4q[csstbox][ip] - l2rec->l1rec->sstref[ip];
4963  }
4964  // the "coldonly" and equatorial aerosol tests are to be run by default, but if SSTMODS is set, don't
4965  if ((evalmask & SSTMODS) == 0) {
4966  /* evaluate change to cold-test only */
4967  /* set the flag bit if sst4 is too much colder than reference */
4968  if (dSST4_ref < -SSTdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
4969  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4970  if (dSST4_ref < -input->sstrefdif &&
4971  l2rec->l1rec->lat[ip] >= equatorialSouth && l2rec->l1rec->lat[ip] <= equatorialNorth &&
4972  l2rec->l1rec->lon[ip] >= equatorialWest && l2rec->l1rec->lon[ip] <= equatorialEast) {
4973  /* tighter test for avhrr between 10S and 30N and -105 to 105 longitude */
4974  /* equatorial aerosol test */
4975  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4976  }
4977  /* set the flag bit if sst4 is too much warmer than reference at night */
4978  if (l2rec->l1rec->solz[ip] >= solznight) {
4979  if (fabs(dSST4_ref) > SSTdiff)
4980  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4981  if (fabs(dSST4_ref) > SSTvdiff)
4982  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4983  }
4984  /* is sst4 way too much colder than reference? */
4985  if (dSST4_ref < -SSTvdiff || l2rec->l1rec->sstref[ip] < sstbad + 1.0)
4986  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4987  } else {
4988  if (l2rec->l1rec->solz[ip] >= SOLZNIGHT) {
4989  if (fabs(dSST4_ref) > SSTdiff)
4990  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4991  } else {
4992  if (dSST4_ref < -SSTdiff || dSST4_ref > (SSTdiff + 1))
4993  flags_sst4[ip] |= SSTF_SSTREFDIFF;
4994  }
4995  if (fabs(dSST4_ref) > SSTvdiff)
4996  flags_sst4[ip] |= SSTF_SSTREFVDIFF;
4997  }
4998  /* check SST4 difference with 11um SST */
4999  dSST_SST4 = sstq[csstbox][ip] - sst4q[csstbox][ip];
5000  if (sstq[csstbox][ip] > sstbad + 1.0) {
5001  if (fabs(dSST_SST4) < SST4diff1)
5002  flags_sst4[ip] |= SSTF_SST4DIFF;
5003  if (fabs(dSST_SST4) < SST4diff2)
5004  flags_sst4[ip] |= SSTF_SST4VDIFF;
5005  }
5006 
5007  /* check sensor zenith limits */
5008  if (l2rec->l1rec->senz[ip] > hisenz)
5009  flags_sst4[ip] |= SSTF_HISENZ;
5010  if (l2rec->l1rec->senz[ip] > vhisenz)
5011  flags_sst4[ip] |= SSTF_VHISENZ;
5012  // flag 2 edge pixels as SSTF_VHISENZ so quality gets set to 3
5013  if (l2rec->l1rec->pixnum[ip] < 2 || l2rec->l1rec->pixnum[ip] > (fullscanpix - 3))
5014  flags_sst4[ip] |= SSTF_VHISENZ;
5015  // set the last 4 pixels of the scan for Terra to VHISENZ
5016  // as there is an apparent calibration issue with the BT12 for those pixels
5017  if ((l2rec->l1rec->l1file->sensorID == MODIST) && (l2rec->l1rec->pixnum[ip] > 1349))
5018  flags_sst4[ip] |= SSTF_VHISENZ;
5019 
5020  // /* if 11um SST is cold, check for clouds (collect 5) */
5021  // /* sstcloud checks to make sure it's day and not glint */
5022  // if (sstq[csstbox][ip] > sstbad+1.0 && l2rec->l1rec->sstref[ip] > sstbad+1.0 &&
5023  // sstq[csstbox][ip]-l2rec->l1rec->sstref[ip] <= -1.0)
5024  // if (sstcloud(ip,cldbox,cldbox,cldthresh) == 1)
5025  // flags_sst4[ip] |= SSTF_REDNONUNIF;
5026 
5027  /* check homogeneity of BT */
5028 
5029  if (Bt39_maxmin[ip] > Bt39unif1)
5030  flags_sst4[ip] |= SSTF_BTNONUNIF;
5031  if (Bt39_maxmin[ip] > Bt39unif2)
5032  flags_sst4[ip] |= SSTF_BTVNONUNIF;
5033 
5034  if (Bt40_maxmin[ip] > Bt40unif1)
5035  flags_sst4[ip] |= SSTF_BTNONUNIF;
5036  if (Bt40_maxmin[ip] > Bt40unif2)
5037  flags_sst4[ip] |= SSTF_BTVNONUNIF;
5038  /* end of homogeneity checks */
5039 
5040  if (l2rec->l1rec->solz[ip] >= solznight) {
5041  if (l2rec->l1rec->l1file->sensorID == MODIST) {
5042  /* night TERRA SST tree test */
5043  treesum[ip] = 0.0;
5044  if (dBt_37_12 < -0.053) {
5045  treesum[ip] += -1.257;
5046  if (l2rec->l1rec->lat[ip] < 33.195) {
5047  treesum[ip] += -0.278;
5048  if (l2rec->l1rec->lat[ip] < -40.185) {
5049  treesum[ip] += 0.619;
5050  } else {
5051  treesum[ip] += -0.711;
5052  if (Bt37 < 6.477) {
5053  treesum[ip] += -3.733;
5054  } else {
5055  treesum[ip] += -0.111;
5056  }
5057  }
5058  } else {
5059  treesum[ip] += 0.333;
5060  }
5061  if (Bt37 < 9.372) {
5062  treesum[ip] += -0.292;
5063  } else {
5064  treesum[ip] += 0.764;
5065  }
5066  } else {
5067  treesum[ip] += 0.430;
5068  if (Bt11_maxmin[ip] < 0.486 && Bt11_maxmin[ip] > (BAD_FLT + 0.1)) {
5069  treesum[ip] += 0.628;
5070  if (Bt40_stdev[ip] < 0.146) {
5071  treesum[ip] += 0.177;
5072  } else {
5073  treesum[ip] += -0.723;
5074  }
5075  } else {
5076  treesum[ip] += -0.450;
5077  }
5078  if (dSST_SST4 < -0.878) {
5079  treesum[ip] += -1.353;
5080  if (dSST_SST4 < -1.533) {
5081  treesum[ip] += -1.439;
5082  } else {
5083  treesum[ip] += 0.346;
5084  }
5085  } else {
5086  treesum[ip] += 0.219;
5087  if (Bt37_stdev[ip] < 0.448) {
5088  treesum[ip] += 0.290;
5089  if (dSST_SST4 < -0.422) {
5090  treesum[ip] += -0.504;
5091  } else {
5092  treesum[ip] += 0.268;
5093  }
5094  } else {
5095  treesum[ip] += -0.484;
5096  }
5097  if (Bt12 < 16.736) {
5098  treesum[ip] += -0.285;
5099  if (dBt_40_11 < -2.199) {
5100  treesum[ip] += 0.518;
5101  } else {
5102  treesum[ip] += -0.316;
5103  if (Bt12 < 11.896) {
5104  treesum[ip] += -0.527;
5105  } else {
5106  treesum[ip] += 0.400;
5107  }
5108  }
5109  } else {
5110  treesum[ip] += 0.500;
5111  }
5112  if (dSST_SST4 < 1.183) {
5113  treesum[ip] += 0.051;
5114  } else {
5115  treesum[ip] += -0.898;
5116  }
5117  }
5118  }
5119 
5120  } else if (l2rec->l1rec->l1file->sensorID == MODISA) {
5121  /* MODIS AQUA SST4 Night decision Tree test */
5122  /* night AQUA SST tree test */
5123  /* use tree that starts with 0, not -0.029 */
5124  treesum[ip] += 0.0;
5125  if (dBt_37_12 < 0.117) {
5126  treesum[ip] += -1.279;
5127  if (l2rec->l1rec->lat[ip] < 33.035) {
5128  treesum[ip] += -0.289;
5129  if (l2rec->l1rec->lat[ip] < -42.235) {
5130  treesum[ip] += 0.747;
5131  } else {
5132  treesum[ip] += -0.564;
5133  if (Bt37 < 13.117) {
5134  treesum[ip] += -0.580;
5135  } else {
5136  treesum[ip] += 0.638;
5137  }
5138  }
5139  } else {
5140  treesum[ip] += 0.355;
5141  }
5142  if (Bt37 < 9.447) {
5143  treesum[ip] += -0.307;
5144  } else {
5145  treesum[ip] += 0.747;
5146  }
5147  } else {
5148  treesum[ip] += 0.470;
5149  if (Bt11_stdev[ip] < 0.155) {
5150  treesum[ip] += 0.690;
5151  if (Bt37_maxmin[ip] < 0.524 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
5152  treesum[ip] += 0.150;
5153  } else {
5154  treesum[ip] += -0.794;
5155  }
5156  } else {
5157  treesum[ip] += -0.430;
5158  }
5159  if (dSST_SST4 < -0.787) {
5160  treesum[ip] += -1.404;
5161  if (dSST_SST4 < -1.253) {
5162  treesum[ip] += -1.086;
5163  } else {
5164  treesum[ip] += 0.388;
5165  }
5166  } else {
5167  treesum[ip] += 0.197;
5168  if (Bt37_maxmin[ip] < 1.229 && Bt37_maxmin[ip] > (BAD_FLT + 0.1)) {
5169  treesum[ip] += 0.287;
5170  if (dSST_SST4 < -0.383) {
5171  treesum[ip] += -0.531;
5172  } else {
5173  treesum[ip] += 0.279;
5174  }
5175  } else {
5176  treesum[ip] += -0.497;
5177  }
5178  if (Bt12 < 16.353) {
5179  treesum[ip] += -0.300;
5180  if (dBt_40_11 < -2.171) {
5181  treesum[ip] += 0.516;
5182  } else {
5183  treesum[ip] += -0.395;
5184  if (sst4q[csstbox][ip] < 16.212) {
5185  treesum[ip] += -0.694;
5186  } else {
5187  treesum[ip] += 0.392;
5188  }
5189  }
5190  } else {
5191  treesum[ip] += 0.451;
5192  }
5193  }
5194  // maybe add check to make sure Bt85 is valid? if ((Bt85 +/- 1000)
5195  if (dBt_85_11 < -1.577) {
5196  treesum[ip] += 0.088;
5197  } else {
5198  treesum[ip] += -0.408;
5199  }
5200  }
5201  }
5202  if (treesum[ip] <= 0.0) {
5203  flags_sst4[ip] |= SSTF_CLOUD;
5204  }
5205  } /* end case for sst4 */
5206  } /* end pixel loop */
5207  return;
5208 }
5209 
5210 /* ----------------------------------------------------------------------------------- */
5211 /* set_qual_sst() - set quality levels for long wave sea surface temperature */
5212 /* */
5213 /* B. Franz, SAIC, August 2005. */
5214 
5215 /* ----------------------------------------------------------------------------------- */
5216 void set_qual_sst(l2str *l2rec) {
5217  int32_t ip, ipb;
5218  float rhot;
5219 
5220  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
5221 
5222  if (l2rec->l1rec->l1file->sensorID == AVHRR) {
5223  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5224  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5225 
5226  qual_sst[ip] = 8;
5227 
5228  } else if ((flags_sst[ip] & SSTF_BTRANGE) == 0
5229  && (flags_sst[ip] & SSTF_BTVNONUNIF) == 0
5230  && (flags_sst[ip] & SSTF_VHISENZ) == 0
5231  && (flags_sst[ip] & SSTF_SSTRANGE) == 0
5232  && (flags_sst[ip] & SSTF_SUNLIGHT) == 0) {
5233 
5234  if ((flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5235  && (flags_sst[ip] & SSTF_BTNONUNIF) == 0
5236  && (flags_sst[ip] & SSTF_HISENZ) == 0
5237  && (flags_sst[ip] & SSTF_GLINT) == 0
5238  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5239  qual_sst[ip] = 0;
5240  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0
5241  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5242  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5243  qual_sst[ip] = 1;
5244  } else if ((flags_sst[ip] & SSTF_HISENZ) == 0
5245  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0
5246  && (flags_sst[ip] & SSTF_CLOUD) == 0) {
5247  qual_sst[ip] = 2;
5248  } else if ((flags_sst[ip] & SSTF_CLOUD) == 0
5249  && (flags_sst[ip] & SSTF_SSTREFDIFF) == 0) {
5250  qual_sst[ip] = 3;
5251  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0
5252  && (flags_sst[ip] & SSTF_HISENZ) == 0) {
5253  qual_sst[ip] = 4;
5254  } else if ((flags_sst[ip] & SSTF_BTNONUNIF) == 0) {
5255  qual_sst[ip] = 5;
5256  } else if ((flags_sst[ip] & SSTF_HISENZ) == 0) {
5257  qual_sst[ip] = 6;
5258  } else {
5259  qual_sst[ip] = 7;
5260  }
5261  if (l2rec->l1rec->ssttype[ip] == 0 && qual_sst[ip] < 4) {
5262  qual_sst[ip] = 4;
5263  }
5264  } else {
5265  qual_sst[ip] = 7;
5266  }
5267  } else {
5268 
5269  if (l2rec->l1rec->solz[ip] < solznight) {
5270 
5271  /* daytime 11um SST quality level */
5272 
5273  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5274  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5275 
5276  qual_sst[ip] = 4;
5277 
5278  } else if ((flags_sst[ip] & SSTF_VHISENZ) > 0
5279  || (flags_sst[ip] & SSTF_BTRANGE) > 0
5280  || (flags_sst[ip] & SSTF_SSTRANGE) > 0
5281  || (flags_sst[ip] & SSTF_BTVNONUNIF) > 0
5282  || (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0
5283  || (flags_sst[ip] & SSTF_CLOUD) > 0
5284  || l2rec->l1rec->ssttype[ip] == 0) {
5285 
5286  qual_sst[ip] = 3;
5287 
5288  } else if (((input->viirsnv7 >= 1)
5289  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5290  || (flags_sst[ip] & SSTF_SSTREFDIFF) > 0
5291  || (flags_sst[ip] & SSTF_REDNONUNIF) > 0) {
5292 
5293  qual_sst[ip] = 2;
5294 
5295  } else if (((input->viirsnv7 <= 0)
5296  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5297  || (l2rec->l1rec->glint_coef[ip] > glintmax)
5298  || (flags_sst[ip] & SSTF_HISENZ) > 0) {
5299  /* if this changes then change comp_sses_sstv6mv */
5300  /* because we use qual 0 sses stats for qual 1 due to glint only */
5301 
5302  qual_sst[ip] = 1;
5303 
5304  } else {
5305 
5306  qual_sst[ip] = 0;
5307 
5308  }
5309 
5310  // Kay thinks this was supposed to be taken out for modis also, but definitely don't want it for viirs?
5311  //if (l2rec->l1rec->l1file->sensorID == MODIST || l2rec->l1rec->l1file->sensorID == MODISA) {
5312  /* Reduce quality if red-band reflectance is high and sst is cold (collect 5) */
5313  /* only in non glint areas */
5314  if (haveRed && l2rec-> l1rec->glint_coef[ip] <= glintmax) {
5315  ipb = ip * nbvis + ibred;
5316  /* Sue: should it divide by csolz here as it does everywhere else?:
5317  * rhot = PI * l2rec->l1rec->Lt[ipb] / l2rec->Fo[ibred] / l2rec->csolz[ip];
5318  * but then the 0.05 limit needs to change?
5319  */
5320  rhot = PI * l2rec->l1rec->Lt[ipb] / l2rec->l1rec->Fo[ibred];
5321  if (qual_sst[ip] < 3 && sstq[csstbox][ip] - l2rec->l1rec->sstref[ip] < -1.0
5322  && rhot > 0.05) {
5323  qual_sst[ip]++;
5324  if ((input->viirsnv7 >= 1)
5325  && (qual_sst[ip] == 1)) {
5326  /* viirsnv7 test want only hisenz to be qual 1 so degrade q0 to q2 */
5327  qual_sst[ip]++;
5328  }
5329 
5330  }
5331  }
5332  //}
5333  } else {
5334  /* nighttime 11um SST quality level */
5335 
5336  if ((flags_sst[ip] & SSTF_ISMASKED) > 0
5337  || (flags_sst[ip] & SSTF_BTBAD) > 0) {
5338 
5339  qual_sst[ip] = 4;
5340 
5341  } else if ((flags_sst[ip] & SSTF_BTRANGE) > 0
5342  || (flags_sst[ip] & SSTF_SSTRANGE) > 0
5343  || (flags_sst[ip] & SSTF_BT4REFDIFF) > 0
5344  || (flags_sst[ip] & SSTF_BTVNONUNIF) > 0
5345  || (flags_sst[ip] & SSTF_CLOUD) > 0
5346  || (flags_sst[ip] & SSTF_VHISENZ) > 0
5347  || (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0) {
5348  /* BTVNONUNIF should be qual 3 as it is for day */
5349  /* don't need to check for haveSST4 because BT4REFDIF is only set if haveSST4 */
5350 
5351  qual_sst[ip] = 3;
5352 
5353  } else if ((flags_sst[ip] & SSTF_SSTREFDIFF) > 0
5354  || (flags_sst[ip] & SSTF_SST4VDIFF) > 0
5355  || ((input->viirsnv7 >= 1)
5356  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5357  || ((input->viirsnv7 >= 1)
5358  && ((flags_sst[ip] & SSTF_SST4DIFF) > 0))
5359  || ((input->viirsnv7 >= 1)
5360  && ((flags_sst[ip] & SSTF_SST3DIFF) > 0))
5361  || (flags_sst[ip] & SSTF_SST3VDIFF) > 0) {
5362  /* for now sst4vdiff and sst3vdiff are the same bit */
5363  /* BTNONUNIF should be qual 2 as it is for day */
5364 
5365  qual_sst[ip] = 2;
5366 
5367  } else if (((input->viirsnv7 <= 0)
5368  && ((flags_sst[ip] & SSTF_BTNONUNIF) > 0))
5369  || ((input->viirsnv7 <= 0)
5370  && ((flags_sst[ip] & SSTF_SST4DIFF) > 0))
5371  || ((input->viirsnv7 <= 0)
5372  && ((flags_sst[ip] & SSTF_SST3DIFF) > 0))
5373  || (flags_sst[ip] & SSTF_HISENZ) > 0) {
5374  /* for now sst4diff and sst3diff are the same bit */
5375 
5376  qual_sst[ip] = 1;
5377 
5378  } else {
5379 
5380  qual_sst[ip] = 0;
5381  }
5382 
5383  /* reduce quality if 4um BTs show non-uniformity (collect 5) */
5384  if (haveSST4 && sst4q[csstbox][ip] > sstbad + 1.0) {
5385  if (qual_sst[ip] < 3
5386  && (flags_sst4[ip] & SSTF_BTNONUNIF) > 0)
5387  qual_sst[ip]++;
5388  if ((input->viirsnv7 >= 1) && (qual_sst[ip] == 1)) {
5389  /* viirsnv7 test want only hisenz to be qual 1 so degrade q0 to q2 */
5390  qual_sst[ip]++;
5391  }
5392  }
5393 
5394  } /* end night sst quality */
5395 
5396  } /* end modis and viirs sst quality */
5397  } /* end pixel loop */
5398  return;
5399 }
5400 
5401 /* ----------------------------------------------------------------------------------- */
5402 /* set_qual_sst4() - set quality levels for short wave sea surface temperature */
5403 /* */
5404 /* B. Franz, SAIC, August 2005. */
5405 
5406 /* ----------------------------------------------------------------------------------- */
5407 void set_qual_sst4(l2str *l2rec) {
5408  int32_t ip;
5409 
5410  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
5411 
5412  if (l2rec->l1rec->solz[ip] < solznight) {
5413 
5414  /* daytime 4um SST quality level */
5415 
5416  if ((flags_sst4[ip] & SSTF_ISMASKED) > 0
5417  || (flags_sst4[ip] & SSTF_BTBAD) > 0) {
5418 
5419  qual_sst4[ip] = 4;
5420 
5421  } else {
5422 
5423  qual_sst4[ip] = 3;
5424 
5425  }
5426 
5427  } else {
5428 
5429  /* nighttime 4um SST quality level */
5430 
5431  if ((flags_sst4[ip] & SSTF_ISMASKED) > 0
5432  || (flags_sst4[ip] & SSTF_BTBAD) > 0) {
5433 
5434  qual_sst4[ip] = 4;
5435 
5436  } else if ((flags_sst4[ip] & SSTF_BTRANGE) > 0
5437  || (flags_sst4[ip] & SSTF_SSTRANGE) > 0
5438  || (flags_sst4[ip] & SSTF_BT4REFDIFF) > 0
5439  || (flags_sst4[ip] & SSTF_BTVNONUNIF) > 0
5440  || (flags_sst4[ip] & SSTF_CLOUD) > 0
5441  || (flags_sst4[ip] & SSTF_VHISENZ) > 0
5442  || (flags_sst4[ip] & SSTF_SSTREFVDIFF) > 0) {
5443  /* BTVNONUNIF and VHISENZ are very bad */
5444 
5445  qual_sst4[ip] = 3;
5446 
5447  } else if ((flags_sst4[ip] & SSTF_SSTREFDIFF) > 0
5448  || (flags_sst4[ip] & SSTF_SST4VDIFF) > 0) {
5449  /* I'm pretty sure this SSTF_SSTREFDIFF should be in flags_sst4, not flags_sst */
5450 
5451  qual_sst4[ip] = 2;
5452 
5453  } else if ((flags_sst4[ip] & SSTF_BTNONUNIF) > 0
5454  || (flags_sst4[ip] & SSTF_SST4DIFF) > 0
5455  || (flags_sst4[ip] & SSTF_HISENZ) > 0) {
5456 
5457  qual_sst4[ip] = 1;
5458 
5459  } else {
5460 
5461  qual_sst4[ip] = 0;
5462  }
5463  }
5464 
5465  }
5466  return;
5467 }
5468 
5469 /* ----------------------------------------------------------------------------------- */
5470 /* nlsst() - nonlinear sst, long wave sea surface temperature */
5471 /* */
5472 /* B. Franz, SAIC, August 2005. */
5473 
5474 /* ----------------------------------------------------------------------------------- */
5475 float nlsst(float Bt11, float Bt12, int32_t is, float sstref,
5476  float **bounds, float **coef, float sstrefoffday,
5477  float sstrefoffnight, int32_t ip, int32_t xsatid,
5478  size_t nlatbands, size_t ncoefficients) {
5479  extern l1qstr l1que;
5480  static float dBtlo = 0.5;
5481  static float dBthi = 0.9;
5482 
5483  float lat = l1que.r[is].lat[ip];
5484  float mu = l1que.r[is].csenz[ip];
5485 
5486  int32_t sensorID = l1que.r[is].l1file->sensorID;
5487  int mside = l1que.r[is].mside;
5488 
5489  float dBt = Bt11 - Bt12; /* aka "thing" */
5490  float sstlo, ssthi, lsst;
5491  float satzdir;
5492  int i, j;
5493 
5494  /* it doesn't hurt, but V5 and viirs osi-saf do not have latband coeffs so they don't use 'ic' */
5495  /* bounds go from -90 to 90 so we want the first set */
5496  /* that has max lat greater than pixel lat */
5497  int ic = nlatbands - 1;
5498  while ((ic > 0) && (lat <= bounds[ic - 1][1] + latwin)) {
5499  ic--;
5500  }
5501 
5502  /* convert sstref to K or whatever offset is required for the coeffs - usually 0.0 except for VIIRS */
5503  if (l1que.r[is].solz[ip] >= solznight)
5504  sstref = sstref + sstrefoffnight;
5505  else
5506  sstref = sstref + sstrefoffday;
5507 
5508  /* sstref is supposed to warm the retrieved temperature, not cool it */
5509  /* this is the scaling factor for temp deficit term */
5510  /* not a problem for deg F or K, but deg C can go negative and we don't want that */
5511  if (sstref < 0.0) {
5512  sstref = 0.0;
5513  }
5514  // VIIRS coefficients made for brightness temperatures are in Kelvin...so add the offset
5515  if (sensorID == VIIRSN || sensorID == VIIRSJ1) {
5516  Bt11 += CtoK;
5517  }
5518  /* this WILL now work if working on a subset of the data (spixl or epixl were specified) */
5519  satzdir = (l1que.r[is].pixnum[ip] < fullscanpix / 2) ? -1.0 : 1.0;
5520 
5521  if (isV5) {
5522  /* Noaa-16 uses different bounds than all the other AVHRR satellites */
5523  if (xsatid == NO16) {
5524  dBtlo = 1.0;
5525  dBthi = 1.4;
5526  }
5527  if(nlatbands < 2) {
5528  printf("ERROR - need nlatbands >= 2, found %ld\n", nlatbands);
5529  exit(EXIT_FAILURE);
5530  }
5531  if(ncoefficients < 4) {
5532  printf("ERROR - need ncoefficients >= 4, found %ld\n", ncoefficients);
5533  exit(EXIT_FAILURE);
5534  }
5535  if (dBt <= dBtlo)
5536  lsst = coef[0][0] + coef[0][1] * Bt11 + coef[0][2] * dBt * sstref + coef[0][3] * dBt * ((1.0 / mu) - 1.0);
5537  else if (dBt >= dBthi)
5538  lsst = coef[1][0] + coef[1][1] * Bt11 + coef[1][2] * dBt * sstref + coef[1][3] * dBt * ((1.0 / mu) - 1.0);
5539  else {
5540  sstlo = coef[0][0] + coef[0][1] * Bt11 + coef[0][2] * dBt * sstref + coef[0][3] * dBt * ((1.0 / mu) - 1.0);
5541  ssthi = coef[1][0] + coef[1][1] * Bt11 + coef[1][2] * dBt * sstref + coef[1][3] * dBt * ((1.0 / mu) - 1.0);
5542  lsst = sstlo + ((dBt - dBtlo) / (dBthi - dBtlo))*(ssthi - sstlo);
5543  }
5544 
5545  } else if (input->viirsnosisaf == 1) {
5546  if (l1que.r[is].solz[ip] >= solznight)
5547  ic = 2;
5548  else
5549  ic = 0;
5550  /* OSI-SAF only has one set of coeffs for day and one for night, not two */
5551  /* OSI-SAF equation: Ts = b0 + (b1 + bt2 Stheta) T11 + [b3 + b4 (Ts0 - 273.15) + b5 Stheta] dT11-12 + b6 Stheta */
5552  /* Ts0 is in degrees Kelvin. Our sstref is in Deg C so don't need to subract 273.15 */
5553  if(ncoefficients < 7) {
5554  printf("ERROR - need ncoefficients >= 7, found %ld\n", ncoefficients);
5555  exit(EXIT_FAILURE);
5556  }
5557  lsst = coef[ic][0]
5558  + (coef[ic][1] + coef[ic][2] * ((1.0 / mu) - 1.0)) * Bt11
5559  + (coef[ic][3] + coef[ic][4] * sstref + coef[ic][5] * ((1.0 / mu) - 1.0)) * dBt
5560  + coef[ic][6] * ((1.0 / mu) - 1.0);
5561 
5562  } else if (input->viirsnv7 >= 1) {
5563  /* v7 high senz latband coeffs */
5564  /* choose coeffs by latband */
5565 
5566  if (ncoefficients < 7) {
5567  printf("ERROR - need ncoefficients >= 7, found %ld\n", ncoefficients);
5568  exit(EXIT_FAILURE);
5569  }
5570  if (lat < bounds[ic][1] - latwin || ic == (nlatbands - 1)) {
5571  lsst = coef[ic][0] + coef[ic][1] * Bt11 + coef[ic][2] * dBt * sstref
5572  + coef[ic][3] * dBt * ((1.0 / mu) - 1.0)
5573  + coef[ic][4] * ((1.0 / mu) - 1.0)
5574  + coef[ic][5] * pow(((1.0 / mu) - 1.0), coef[ic][6]);
5575 
5576  } else {
5577  dBtlo = bounds[ic][1] - latwin;
5578  dBthi = bounds[ic][1] + latwin;
5579  sstlo = coef[ic][0] + coef[ic][1] * Bt11 + coef[ic][2] * dBt * sstref
5580  + coef[ic][3] * dBt * ((1.0 / mu) - 1.0)
5581  + coef[ic][4] * ((1.0 / mu) - 1.0)
5582  + coef[ic][5] * pow(((1.0 / mu) - 1.0), coef[ic][6]);
5583  ssthi = coef[ic + 1][0] + coef[ic + 1][1] * Bt11
5584  + coef[ic + 1][2] * dBt * sstref
5585  + coef[ic + 1][3] * dBt * ((1.0 / mu) - 1.0)
5586  + coef[ic + 1][4] * ((1.0 / mu) - 1.0)
5587  + coef[ic + 1][5] * pow(((1.0 / mu) - 1.0), coef[ic + 1][6]);
5588 
5589  lsst = sstlo + ((lat - dBtlo) / (dBthi - dBtlo)) * (ssthi - sstlo);
5590 
5591  }
5592 
5593  } else {
5594  /* v6 latband coeffs */
5595 
5596  /* Only coeff 4 for terra and aqua, they're zeroed for others */
5597  /* Only coeffs 5-6 for terra, aqua and viirs; they're zeroed for avhrr */
5598  static float* C0;
5599  static float* C1;
5600  static int last_ic = -1;
5601 
5602  if (C0 == NULL) {
5603  C0 = (float *) calloc(7, sizeof(float));
5604  C1 = (float *) calloc(7, sizeof(float));
5605  }
5606  if (last_ic != ic) {
5607  for (i = 0; i < ncoefficients; i++){
5608  j = i;
5609  if (((sensorID == VIIRSN || sensorID == VIIRSJ1) || (sensorID == AVHRR)) && (i > 3)){
5610  j += 1;
5611  }
5612  C0[j] = coef[ic][i];
5613  if (ic < nlatbands -1)
5614  C1[j] = coef[ic + 1][i];
5615  }
5616  }
5617  if (lat < bounds[ic][1] - latwin || ic == (nlatbands - 1)) {
5618  lsst = C0[0]
5619  + C0[1] * Bt11
5620  + C0[2] * dBt * sstref
5621  + C0[3] * dBt * ((1.0 / mu) - 1.0)
5622  + C0[4] * mside
5623  + C0[5] * (l1que.r[is].senz[ip] * satzdir)
5624  + C0[6] * pow((l1que.r[is].senz[ip] * satzdir), 2);
5625  } else {
5626  dBtlo = bounds[ic][1] - latwin;
5627  dBthi = bounds[ic][1] + latwin;
5628  sstlo = C0[0]
5629  + C0[1] * Bt11
5630  + C0[2] * dBt * sstref
5631  + C0[3] * dBt * ((1.0 / mu) - 1.0)
5632  + C0[4] * mside
5633  + C0[5] * (l1que.r[is].senz[ip] * satzdir)
5634  + C0[6] * pow((l1que.r[is].senz[ip] * satzdir), 2);
5635  ssthi = C1[0]
5636  + C1[1] * Bt11
5637  + C1[2] * dBt * sstref
5638  + C1[3] * dBt * ((1.0 / mu) - 1.0)
5639  + C1[4] * mside
5640  + C1[5] * (l1que.r[is].senz[ip] * satzdir)
5641