OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
DtAlgOcean.cpp
Go to the documentation of this file.
1 /*******************************************************************************
2  *
3  * NAME: DtAlgOcean.cpp
4  *
5  * DESCRIPTION: Object class that provides data structures and processes that
6  * compute processings for a given DDProcess object class.
7  *
8  * Created on: November 3, 2016
9  * Author: Sam Anderson, DT
10  *
11  * Modified:
12  *
13  *******************************************************************************/
14 
15 #include <math.h>
16 #include <cstring>
17 #include <fstream>
18 
19 #include <DDProcess.h>
20 #include "darktarget/DtLutNetcdf.h"
21 #include "darktarget/DtMask.h"
22 #include "darktarget/DtAlgorithm.h"
23 
24 using namespace std;
25 
26 /*******************************************************************************
27 *
28 * DESCRIPTION:
29 * The MODIS ocean aerosol product consists of aerosol optical thickness
30 * and size parameters estimates derived on a 10x10 (1-km) pixel spatial
31 * array. The measured radiances in a wide spectral range (0.47-2.13
32 * microns) are inverted assuming a bi-lognormal size distribution.
33 * The volume and the mean particle size for each log-normal mode are
34 * determined. When fully developed, the aerosol ocean algorithm with
35 * use the seven MODIS bands: 1, 2, 3, 4, 5, 6, and 7.
36 *
37 * INPUT PARAMETERS: NONE
38 *
39 * OUTPUT PARAMETERS:
40 * SDS_ref_STD Standard deviation of reflectances at 7 bands
41 * SDSTAU_best Optical thickness for best solution
42 * SDSTAUS_best Optical thickness contribution small particles for best solution
43 * SDSTAUB_best Optical thickness contribution large particles for best solution
44 * SDSTAU_average Optical thickness for best solution
45 * SDSTAUS_average Optical thickness contribution small particles for best solution
46 * SDSTAUB_average Optical thickness contribution large particles for best solution
47 * SDS_Least_error Least square error estimated
48 * SDS_small_weighting Small mode weighting factor
49 * SDS_sol_INDX_small Solution Number index small particles
50 * SDS_sol_INDX_large Solution Number index large particles
51 * SDSASSY_best Asymmetry_Factor at 7 bands for best solution
52 * SDSASSY_average Asymmetry_Factor at 7 bands for average solution
53 * SDSBAcK_best Backscattering ratio at 7 bands of best solution
54 * SDSBAcK_average Backscattering ratio at 7 bands of average solution
55 * SDS_effrad Effective_Radius at 0.55 micron of both solutions
56 * SDS_AOT_model Ratio of optical depth of small mode vs effective optical depth at 550
57 * SDS_RefF_best Normalized reflected_flux at 7 bands of best solution
58 * SDS_RefF_average Normalized reflected_flux at 7 bands of average solution
59 * SDS_TranF_best Normalized Transmitted_flux at 7 bands of best solution
60 * SDS_TranF_average Normalized Transmitted_flux at 7 bands of average solution
61 * SDS_SCAT_ANGLE_OCEAN Scattering angle ocean
62 * SDS_QCONTROL Quality control SDS array
63 * SDS_NUMPIXELS Number of Pixels used for 0.55 micron
64 * SDS_ccn cloud_Fraction in percentage
65 * SDS_mass_conc Mass concentration
66 * SDS_angs_coeff1 Angstrom Exponent for 0.550 and 0.865 miron
67 * SDS_angs_coeff2 Angstrom exponent for 0.865 and 2.130 micron
68 *
69 * REVISION HISTORY:
70 * $Log: Process_ocean_V2.f,v $
71 * Revision 2.1 1996/05/28 15:34:05 vlin
72 * Updated to write out one scancube at a time
73 *
74 * TEAM-UNIQUE HEADER:
75 *
76 * *REFERENcES AND CREDITS:
77 *
78 * WRITTEN BY:
79 * SHANA MATTOO E-MAIL:mattoo@climate.gsfc.nasa.gov
80 * APPLIED RESEARCH CORPORATION PHONE: (301) 614-6214
81 * NASA/GODDARD SPAcE FLIGHT cENTER FAX: (301) 614-6307
82 * CODE 913
83 * GREENBELT, MD 20771
84 *
85 * DESIGN NOTES:
86 *
87 * Internals Variables:
88 *
89 * Reflectances of Wavelengths used modis,ch1-ch7 identified by REFW*
90 *
91 * REFW670
92 * REFW865
93 * REFW488
94 * REFW550
95 * REFW1240
96 * REFW1610
97 * REFW2250
98 * CLDMSK cloud mask 0=cloudy,1=clear
99 * MTHET0 Measured solar zenith angle in deg.
100 * MTHET Measured viewangle from ground in deg.
101 * MPHI Measured azimuth in deg.
102 * NWAV Number of wavelengths.
103 * NAOT Number of optical thicknesses.
104 * NTH0 Number of solar zenith angles.
105 * NTHET Number of view angles.
106 * NPHI Number of azimuth angles.
107 * PHc,JPHI Azimuth angles.
108 * THET View angles.
109 * THET0 Solar zenith angles.
110 * TVALUE Transmission factor.
111 * AINT Reflectance from look-up table.
112 * TAUA optical thicknesses.
113 * Numdata Number of input data sets.
114 *
115 *************************************************************************/
116 
117 /**************************************************************************
118  * NAME: DtAlgOcean()
119  *
120  * DESCRIPTION: Class Constructor
121  *
122  *************************************************************************/
123 
124 DtAlgOcean::DtAlgOcean()
125 {
126  cm_ = nullptr;
127  memset(refl_rayl_, 0, NWAV*sizeof(float));
128  memset(refl_big_, 0, NWAV*NUM_CASES_BIG*NAOT*sizeof(float));
129  memset(refl_small_, 0, NWAV*NUM_CASES_SMALL*NAOT*sizeof(float));
130  memset(albedo_R_small_tau_, 0, NWAV*NUM_CASES_SMALL*NAOT*sizeof(float));
131  memset(albedo_R_big_tau_, 0, NWAV*NUM_CASES_BIG*NAOT*sizeof(float));
132  memset(albedo_T_small_tau_, 0, NWAV*NUM_CASES_SMALL*NAOT*sizeof(float));
133  memset(albedo_T_big_tau_, 0, NWAV*NUM_CASES_BIG*NAOT*sizeof(float));
134  memset(refl_, 0, NWAV*sizeof(float));
135  memset(sdev_, 0, NWAV*sizeof(float));
136  memset(numData_, 0, NWAV*sizeof(int));
137  memset(good_pixels_, 0, NWAV*sizeof(int));
138  memset(tau_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
139  memset(tau_small_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
140  memset(tau_big_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
141  memset(backscatter_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
142  memset(assym_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
143  memset(refl_flux_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
144  memset(trans_flux_, 0, NWAV*NUMCASES*NUMCASEB*sizeof(float));
145  memset(angstrom_exp_, 0, NUM_STATS*NUMCASES*NUMCASEB*sizeof(float));
146  memset(solution_index_, 0, NUM_SIZES*NUMCASES*NUMCASEB*sizeof(float));
147  memset(ccn_, 0, NUMCASES*NUMCASEB*sizeof(float));
148  memset(eff_radius_, 0, NUMCASES*NUMCASEB*sizeof(float));
149  memset(eff_variance_, 0, NUMCASES*NUMCASEB*sizeof(float));
150  memset(mass_con_ocean_, 0, NUMCASES*NUMCASEB*sizeof(float));
151  memset(xmin_, 0, 2*NUMCASES*NUMCASEB*sizeof(float));
152  memset(funmin_, 0, 2*NUMCASES*NUMCASEB*sizeof(float));
153  memset(tau_X55_, 0, 2*NUMCASES*NUMCASEB*sizeof(float));
154  memset(error_, 0, 2*NWAV*NUMCASES*NUMCASEB*sizeof(float));
155  memset(refl_avg_, 0, NWAV*sizeof(float));
156  memset(sdev_avg_, 0, NWAV*sizeof(float));
157  memset(tau_avg_, 0, NUM_STATS*NWAV*sizeof(float));
158  memset(tau_small_avg_, 0, NUM_STATS*NWAV*sizeof(float));
159  memset(tau_big_avg_, 0, NUM_STATS*NWAV*sizeof(float));
160  memset(backscatter_avg_, 0, NUM_STATS*NWAV*sizeof(float));
161  memset(assym_avg_, 0, NUM_STATS*NWAV*sizeof(float));
162  memset(refl_flux_avg_, 0, NUM_STATS*NWAV*sizeof(float));
163  memset(trans_flux_avg_, 0, NUM_STATS*NWAV*sizeof(float));
164  memset(angstrom_exp_avg_, 0, NUM_STATS*NUM_SIZES*sizeof(float));
165  memset(solution_index_avg_, 0, NUM_STATS*NUM_SIZES*sizeof(float));
166  memset(ccn_avg_, 0, NUM_STATS*sizeof(float));
167  memset(eff_radius_avg_, 0, NUM_STATS*sizeof(float));
168  memset(eff_variance_avg_, 0, NUM_STATS*sizeof(float));
169  memset(mass_con_ocean_avg_, 0, NUM_STATS*sizeof(float));
170  memset(xmin_avg_, 0, NUM_STATS*sizeof(float));
171  memset(funmin_avg_, 0, NUM_STATS*sizeof(float));
172  memset(tau_X55_avg_, 0, NUM_STATS*sizeof(float));
173  qcontrol_= 0;
174  qcontrol_exclude_= 0;
175  qcontrol_special_= 0;
176  qcontrol_cirrus_= 0;
177  quality_dust_flag_glint_= 0;
178  quality_dust_flag_off_glint_= 0;
179  memset(quality_to_pass_, 0, 2*sizeof(short));
180  memset(quality_flag_, 0, 12*sizeof(short));
181  memset(sds_refl_, 0, NWAV*sizeof(float));
182  memset(sds_refl_sdev_, 0, NWAV*sizeof(float));
183  memset(sds_numPixels_, 0, NWAV*sizeof(short));
184  memset(sds_tau_best_, 0, NWAV*sizeof(float));
185  memset(sds_tau_avg_, 0, NWAV*sizeof(float));
186  memset(sds_tau_big_best_, 0, NWAV*sizeof(float));
187  memset(sds_tau_big_avg_, 0, NWAV*sizeof(float));
188  memset(sds_tau_small_best_, 0, NWAV*sizeof(float));
189  memset(sds_tau_small_avg_, 0, NWAV*sizeof(float));
190  memset(sds_assy_best_, 0, NWAV*sizeof(float));
191  memset(sds_assy_avg_, 0, NWAV*sizeof(float));
192  memset(sds_back_best_, 0, NWAV*sizeof(float));
193  memset(sds_back_avg_, 0, NWAV*sizeof(float));
194  memset(sds_reff_best_, 0, NWAV*sizeof(float));
195  memset(sds_reff_avg_, 0, NWAV*sizeof(float));
196  memset(sds_tranf_best_, 0, NWAV*sizeof(float));
197  memset(sds_tranf_avg_, 0, NWAV*sizeof(float));
198  sds_Tau_Land_Ocean_= 0;
199  sds_Tau_Land_Ocean_img_= 0;
200  memset(sds_Small_Weighting_, 0, NUM_SIZES*sizeof(float));
201  memset(sds_Least_Error_, 0, NUM_SIZES*sizeof(float));
202  memset(sds_tau_X55_, 0, NUM_SIZES*sizeof(float));
203  memset(sds_EffRad_, 0, NUM_SIZES*sizeof(float));
204  memset(sds_EffVar_, 0, NUM_SIZES*sizeof(float));
205  memset(sds_Sol_Index_Small_, 0, NUM_SIZES*sizeof(short));
206  memset(sds_Sol_Index_Large_, 0, NUM_SIZES*sizeof(short));
207  memset(sds_Angs_Coeff1_, 0, NUM_SIZES*sizeof(float));
208  memset(sds_Angs_Coeff2_, 0, NUM_SIZES*sizeof(float));
209  memset(sds_Mass_Conc_, 0, NUM_SIZES*sizeof(float));
210  memset(sds_CCN_, 0, NUM_SIZES*sizeof(float));
211  memset(sds_AOT_model_, 0, (NUMCASES+NUMCASEB)*sizeof(float));
212  sds_land_ocean_quality_= 0;
213 }
214 
215 /**************************************************************************
216  * NAME: ~DtAlgOcean()
217  *
218  * DESCRIPTION: Class Destructor
219  *
220  *************************************************************************/
221 
222 DtAlgOcean::~DtAlgOcean()
223 {
224  if (cm_ != 0) {
225  delete cm_;
226  }
227 }
228 
229 /**************************************************************************
230  * NAME: initialize()
231  *
232  * DESCRIPTION: Virtual function initializes data and object classes for
233  * processing operations.
234  *
235  *************************************************************************/
236 
237 int DtAlgOcean::initialize( map<string, ddata*> imap )
238 {
239  int status = DTDB_SUCCESS;
240 
241  qcontrol_special_ = 0;
242 
243  DtAlgorithm::initialize( imap );
244 
245  DtLutNetcdf* lutgen = new DtLutNetcdf();
246  status = lutgen->read_ocean_aerosol_lut( lut_ );
247  delete lutgen;
248 
249  cm_ = new DtCloudMaskOcean(this);
250 
251  return status;
252 }
253 
254 /**************************************************************************
255  * NAME: process()
256  *
257  * DESCRIPTION: Dark target ocean algorithm process.
258  *
259  *************************************************************************/
260 
261 map<string, ddata*> DtAlgOcean::process(vector<size_t> start, vector<size_t> count,
262  map<string, ddata*> imap)
263 {
264  get_inputs(start, count, imap);
265 
266  set_fill_out();
267  if (rfl_[(size_t)rhot_band::W550]<0.0) {
268  l2_flags_ += (unsigned int) flags::BOWTIEDEL;
269  return set_fills();
270  }
271 
272  if((solz_ >= MINMTHET0) && (solz_ <= MAXMTHET0) &&
273  (senz_ >= MINMTHET) && (senz_ <= MAXMTHET) &&
274  (raa_ >= MINMPHI) && (raa_ <= MAXMPHI)) {
275 
276  rfld_[D412] = rfl_[(size_t)rhot_band::W410];
277  rfld_[D488] = rfl_[(size_t)rhot_band::W490];
278  rfld_[D550] = rfl_[(size_t)rhot_band::W550];
279  rfld_[D670] = rfl_[(size_t)rhot_band::W670];
280  rfld_[D865] = rfl_[(size_t)rhot_band::W865];
281  rfld_[D1240] = rfl_[(size_t)rhot_band::W1240];
282  rfld_[D1380] = rfl_[(size_t)rhot_band::W1380];
283  rfld_[D1610] = rfl_[(size_t)rhot_band::W1610];
284  rfld_[D2250] = rfl_[(size_t)rhot_band::W2250];
285  if (bgascorrect_) {
286  compute_gas_correction();
287  rfld_[D412] *= gasc_[D412];
288  rfld_[D488] *= gasc_[D488];
289  rfld_[D550] *= gasc_[D550];
290  rfld_[D670] *= gasc_[D670];
291  rfld_[D865] *= gasc_[D865];
292  rfld_[D1240] *= gasc_[D1240];
293  rfld_[D1610] *= gasc_[D1610];
294  rfld_[D2250] *= gasc_[D2250];
295  }
296 
297  if (cloud_mask_ == DFILL_UBYTE) {
298  cm_->compute(cmask_);
299  cloud_mask_ = (cmask_+1)%2;
300  } else {
301  cmask_ = (cloud_mask_+1)%2;
302  }
303  if (!bcloudmask_) cmask_ = 1;
304  if (bcloudmask_ && cloud_mask_) {
305  l2_flags_ += (unsigned int) flags::CLDICE;
306  return set_fills();
307  }
308 
309  compute_glint_angle();
310  index_geometry( solz_, senz_, raa_ );
311  interpolate_rayleigh();
312  compute_avg_refl();
313  store_reflectance();
314 
315  if((!bglintmask_) ||
316  ((refl_[D550] > 0) && (refl_[D670] > 0) &&
317  (refl_[D865] > 0) && (refl_[D2250] > 0) &&
318  (qcontrol_special_== 0 || qcontrol_special_== 2))) {
319 
320  interpolate_angles();
321 
322  int iSol = 0;
323  for (int iSmall=0; iSmall<NUMCASES; iSmall++) {
324  for (int iBig=0; iBig<NUMCASEB; iBig++) {
325 
326 // compute_minimum_baseline(iBig, iSmall, iSol);
327  compute_minimum(iBig, iSmall, iSol);
328  compute_tau_flux(iBig, iSmall, iSol);
329  solution_index_[ISMALL][iSol] = iSmall;
330  solution_index_[IBIG][iSol] = iBig;
331  iSol++;
332  }
333  }
334 
335  average_output();
336  store_output();
337 
338  qcontrol_ = (tau_avg_[1][0] > MAXTAU) ? -29 : qcontrol_;
339  qcontrol_ = (tau_avg_[1][0] < 0) ? 9 : qcontrol_;
340  qcontrol_ = (tau_avg_[1][0] <= MINTAU) ? -28 : qcontrol_;
341  }
342  else if (qcontrol_special_ == 3) qcontrol_ = -21;
343  }
344  else qcontrol_ = -26;
345 
346  assign_quality();
347 
348  // write to generic outputs
349 
350  qual_flag_ = quality_to_pass_[0];
351  aerosol_type_ = (NUMCASEB-1)*sds_Sol_Index_Large_[1] +
352  sds_Sol_Index_Small_[1];
353  error_flag_ = abs(quality_to_pass_[1]);
354 
355  scatter_ang_ = 180.0 - glint_angle_;
356  glint_ang_ = glint_angle_;
357  sse_ = sds_Least_Error_[1];
358  fmf_ = sds_Small_Weighting_[1];
359  aot_550_ = sds_tau_X55_[1];
360  ae1_ = sds_Angs_Coeff1_[1];
361  ae2_ = sds_Angs_Coeff2_[1];
362  ndv_ = DFILL_FLOAT;
363  chlor_ = DFILL_FLOAT;
364  aot_[0] = DFILL_FLOAT;
365  for ( int iWav=0; iWav<NOWL; iWav++ ) {
366  aot_[iWav+1] = sds_tau_avg_[iWav];
367  }
368  for ( int iWav=0; iWav<NLWL+1; iWav++ ) {
369  sr_[iWav] = DFILL_FLOAT;
370  ssa_[iWav] = DFILL_FLOAT;
371  }
372 
373  return set_outputs();
374 }
375 
376 /**************************************************************************
377  * NAME: index_geometry()
378  *
379  * DESCRIPTION: Reset LUT indices based on measured geometry.
380  *
381  *************************************************************************/
382 
383 int DtAlgOcean::index_geometry( float sza, float azim, float phi)
384 {
385  int status = DTDB_SUCCESS;
386 
387  float del = (sza <= lut_.THET0[4]) ? 12.0 : 6.0;
388  SZA_0_ = (sza <= lut_.THET0[4]) ? (int)(sza/del) : (int)(sza/del)-4;
389  SZA_1_ = SZA_0_+1;
390  THE_0_ = (int)(azim/6.0);
391  THE_1_ = THE_0_+1;
392  PHI_0_ = (int)(phi/12.0);
393  PHI_1_ = PHI_0_+1;
394 
395  return status;
396 }
397 
398 /**************************************************************************
399  * NAME: interpolate_rayleigh()
400  *
401  * DESCRIPTION: Subroutine interpolate_rayleigh interpolates the lookup
402  * Rayleigh reflectances to the measured geometry.
403  *
404  *************************************************************************/
405 
406 int DtAlgOcean::interpolate_rayleigh()
407 {
408  int status = DTDB_SUCCESS;
409 
410  int lutindex = 0;
411  while ((wind_[++lutindex] < ws_) &&
412  (lutindex < (WIND_LUT_ENTRIES-1))) {}
413  short iws1 = lutindex-1;
414  short iws2 = lutindex;
415  float ws_lut[WIND_LUT_ENTRIES];
416  memset(ws_lut, 0, WIND_LUT_ENTRIES*sizeof(float));
417  float phi_out[2] = {0.0,0.0};
418  float the_out[2] = {0.0,0.0};
419  float sza_out[2] = {0.0,0.0};
420  for (int iWav=0; iWav<NWAV; iWav++) {
421  int numWS = 0;
422  for (int iWS=iws1; iWS <= iws2; iWS++ ) {
423  ws_lut[numWS] = wind_[iWS];
424  int numSza = 0;
425  for (int iSza=SZA_0_; iSza<=SZA_1_; iSza++) {
426  int numThe = 0;
427  for (int iThe=THE_0_; iThe<=THE_1_; iThe++) {
428  int numPhi = 0;
429  for (int iPhi=PHI_0_; iPhi<=PHI_1_; iPhi++) {
430  numPhi++;
431  }
432  interp_extrap(numPhi, raa_, &lut_.PHC[PHI_0_],
433  &lut_.REF_RAYALL[iWav][iWS][iSza][iThe][PHI_0_],
434  phi_out[numThe]);
435  numThe++;
436  }
437  interp_extrap(numThe, senz_, &lut_.THET[THE_0_],
438  phi_out, the_out[numSza]);
439  numSza++;
440  }
441  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
442  the_out, sza_out[numWS]);
443  numWS++;
444  }
445  interp_extrap(numWS, ws_, ws_lut, sza_out, refl_rayl_[iWav]);
446 
447  refl_rayl_[iWav] *= M_PI/cos(DEGtoRAD*solz_);
448  }
449 
450  return status;
451 }
452 
453 /**************************************************************************
454  * NAME: interpolate_angles()
455  *
456  * DESCRIPTION: Subroutine interpolates the lookup reflectances to the
457  * measured geometry.
458  *
459  *************************************************************************/
460 
461 int DtAlgOcean::interpolate_angles()
462 {
463  int status = DTDB_SUCCESS;
464 
465  int lutindex = 0;
466  while ((wind_[++lutindex] < ws_) &&
467  (lutindex < (WIND_LUT_ENTRIES-1))) {}
468  short iws1 = lutindex-1;
469  short iws2 = lutindex;
470  float ws_lut[WIND_LUT_ENTRIES];
471  memset(ws_lut, 0, WIND_LUT_ENTRIES*sizeof(float));
472 
473  float phi_out[2][2] = {{0.0,0.0},{0.0,0.0}};
474  float the_out[2][2] = {{0.0,0.0},{0.0,0.0}};
475  float sza_out[2][2] = {{0.0,0.0},{0.0,0.0}};
476  float ws_out = 0.0;
477 
478  for (int iWav=0; iWav<NWAV; iWav++) {
479  // INTERPOLATE FOR SMALL CASES
480  for (int iCase=0; iCase<NUMCASES; iCase++) {
481  for (int iTau=0; iTau<NAOT; iTau++) {
482  int numWS = 0;
483  for (int iWS=iws1; iWS <= iws2; iWS++) {
484  ws_lut[numWS] = wind_[iWS];
485  int numSza = 0;
486  for (int iSza=SZA_0_; iSza<=SZA_1_; iSza++) {
487  numSza++;
488  }
489  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
490  &lut_.ALBEDO_R_SMALL[iWav][iCase][iTau][iWS][SZA_0_],
491  sza_out[0][numWS]);
492  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
493  &lut_.ALBEDO_T_SMALL[iWav][iCase][iTau][iWS][SZA_0_],
494  sza_out[1][numWS]);
495  numWS++;
496  }
497  interp_extrap(numWS, ws_, ws_lut,
498  sza_out[0], albedo_R_small_tau_[iWav][iCase][iTau]);
499  interp_extrap(numWS, ws_, ws_lut,
500  sza_out[1], albedo_T_small_tau_[iWav][iCase][iTau]);
501  }
502  }
503  // INTERPOLATE FOR BIG CASES
504  for (int iCase=0; iCase<NUMCASEB; iCase++) {
505  for (int iTau=0; iTau<NAOT; iTau++) {
506  int numWS = 0;
507  for (int iWS=iws1; iWS <= iws2; iWS++) {
508  ws_lut[numWS] = wind_[iWS];
509  int numSza = 0;
510  for (int iSza=SZA_0_; iSza<=SZA_1_; iSza++) {
511  numSza++;
512  }
513  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
514  &lut_.ALBEDO_R_BIG[iWav][iCase][iTau][iWS][SZA_0_],
515  sza_out[0][numWS]);
516  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
517  &lut_.ALBEDO_T_BIG[iWav][iCase][iTau][iWS][SZA_0_],
518  sza_out[1][numWS]);
519  numWS++;
520  }
521  interp_extrap(numWS, ws_, ws_lut,
522  sza_out[0], albedo_R_big_tau_[iWav][iCase][iTau]);
523  interp_extrap(numWS, ws_, ws_lut,
524  sza_out[1], albedo_T_big_tau_[iWav][iCase][iTau]);
525  }
526  }
527  }
528 
529  for (int iWav = 0; iWav < NWAV; iWav++) {
530  for (int iCase = 0; iCase < NUMCASES; iCase++) {
531  for (int iTau = 1; iTau < NAOT; iTau++) {
532  int numWS = 0;
533  for (int iWS = iws1; iWS <= iws2; iWS++) {
534  ws_lut[numWS] = wind_[iWS];
535  int numSza = 0;
536  for (int iSza = SZA_0_; iSza <= SZA_1_; iSza++) {
537  int numThe = 0;
538  for (int iThe = THE_0_; iThe <= THE_1_; iThe++) {
539  int numPhi = 0;
540  for (int iPhi = PHI_0_; iPhi <= PHI_1_; iPhi++) {
541  numPhi++;
542  }
543  interp_extrap(numPhi, raa_, &lut_.PHC[PHI_0_],
544  &lut_.AINTS[iWav][iCase][iTau][iWS][iSza][iThe][PHI_0_],
545  phi_out[ISMALL][numThe]);
546  numThe++;
547  }
548  interp_extrap(numThe, senz_, &lut_.THET[THE_0_],
549  phi_out[ISMALL], the_out[ISMALL][numSza]);
550  numSza++;
551  }
552  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
553  the_out[ISMALL], sza_out[ISMALL][numWS]);
554  numWS++;
555  }
556  interp_extrap(numWS, ws_, ws_lut, sza_out[ISMALL],
557  ws_out);
558  refl_small_[iWav][iCase][iTau] = ws_out * M_PI
559  / cos( solz_*DEGtoRAD);
560  }
561  refl_small_[iWav][iCase][0] = refl_rayl_[iWav];
562  }
563  for (int iCase = 0; iCase < NUMCASEB; iCase++) {
564  for (int iTau = 1; iTau < NAOT; iTau++) {
565  int numWS = 0;
566  for (int iWS = iws1; iWS <= iws2; iWS++) {
567  ws_lut[numWS] = wind_[iWS];
568  int numSza = 0;
569  for (int iSza = SZA_0_; iSza <= SZA_1_; iSza++) {
570  int numThe = 0;
571  for (int iThe = THE_0_; iThe <= THE_1_; iThe++) {
572  int numPhi = 0;
573  for (int iPhi = PHI_0_; iPhi <= PHI_1_; iPhi++) {
574  numPhi++;
575  }
576  interp_extrap(numPhi, raa_, &lut_.PHC[PHI_0_],
577  &lut_.AINTB[iWav][iCase][iTau][iWS][iSza][iThe][PHI_0_],
578  phi_out[IBIG][numThe]);
579  numThe++;
580  }
581  interp_extrap(numThe, senz_, &lut_.THET[THE_0_],
582  phi_out[IBIG], the_out[IBIG][numSza]);
583  numSza++;
584  }
585  interp_extrap(numSza, solz_, &lut_.THET0[SZA_0_],
586  the_out[IBIG], sza_out[IBIG][numWS]);
587  numWS++;
588  }
589  interp_extrap(numWS, ws_, ws_lut, sza_out[IBIG],
590  ws_out);
591  refl_big_[iWav][iCase][iTau] = ws_out * M_PI
592  / cos( solz_*DEGtoRAD);
593  }
594  refl_big_[iWav][iCase][0] = refl_rayl_[iWav];
595  }
596  }
597 
598  return status;
599 }
600 
601 /**************************************************************************
602  * NAME: compute_tau_flux()
603  *
604  * DESCRIPTION: This subroutine computes CCN, assymetry factor,
605  * backscattering ratio, effective radius, effective
606  * variance, and optical thicknesses for large and small
607  * mode.
608  *
609  *************************************************************************/
610 
611 int DtAlgOcean::compute_tau_flux( int iBig, int iSmall, int iSol )
612 {
613  int status = DTDB_SUCCESS;
614 
615  int lutindex = 0;
616  while ((wind_[++lutindex] < ws_) &&
617  (lutindex < (WIND_LUT_ENTRIES-1))) {}
618  short iws1 = lutindex-1;
619  short iws2 = lutindex;
620  int numWS = 0;
621 
622 // Since we see the glint effect around tau of 0.5 no retrievals will be done
623 // less than tau of 0.5 in glint area. It is set to negative value so that later on it is
624 // filled with fill values
625  if((tau_X55_[0][iSol] <= 0.7) && (quality_dust_flag_glint_ == 1)) {
626  tau_X55_[0][iSol] = -0.02;
627  }
628 // If tau at 0.55um is greater than maxtau(ie 5.00) taux55 is set to maxtau for
629 // glint and off glint heavy dust episodes
630  if(((tau_X55_[0][iSol] > MAXTAU) && (quality_dust_flag_glint_== 1)) ||
631  ((tau_X55_[0][iSol] > MAXTAU) &&
632  (quality_dust_flag_off_glint_== 1))) {
633  tau_X55_[0][iSol]= MAXTAU;
634  }
635 // Compute normalized extinction coefficient
636  float tau_X55_small = tau_X55_[0][iSol]*xmin_[0][iSol];
637  float tau_X55_big = tau_X55_[0][iSol]- tau_X55_small;
638 
639  float ws_lut[WIND_LUT_ENTRIES];
640  float rsmall[WIND_LUT_ENTRIES];
641  float rbig[WIND_LUT_ENTRIES];
642  memset(ws_lut, 0, WIND_LUT_ENTRIES*sizeof(float));
643  memset(rsmall, 0, WIND_LUT_ENTRIES*sizeof(float));
644  memset(rbig, 0, WIND_LUT_ENTRIES*sizeof(float));
645 
646  for (int iWav=0; iWav<NWAV; iWav++) {
647  numWS = 0;
648  for (int iWS=iws1; iWS <= iws2; iWS++) {
649  ws_lut[numWS] = wind_[iWS];
650  rsmall[numWS] =
651  lut_.EXTSMALL[iWav][iSmall][iWS]/lut_.EXTSMALL[D550][iSmall][iWS];
652  rbig[numWS] =
653  lut_.EXTBIG[iWav][iBig][iWS]/lut_.EXTBIG[D550][iBig][iWS];
654  numWS++;
655  }
656  float extsmall = 0.0;
657  interp_extrap(numWS,ws_,ws_lut,rsmall,extsmall);
658  float extbig = 0.0;
659  interp_extrap(numWS,ws_,ws_lut,rbig,extbig);
660 
661  tau_[iWav][iSol] = tau_X55_small*extsmall + tau_X55_big*extbig;
662  tau_small_[iWav][iSol] = tau_X55_small*extsmall;
663  tau_big_[iWav][iSol] = tau_X55_big*extbig;
664  }
665 
666  if ((tau_[D550][iSol] > 0) && (tau_[D865][iSol] > 0)) {
667  float tlog = log(tau_[D550][iSol]/tau_[D865][iSol]) /
668  log(lut_.WAVE[D865]/lut_.WAVE[D550]);
669  angstrom_exp_[0][iSol] = tlog;
670  }
671  if ((tau_[D865][iSol] > 0) && (tau_[D2250][iSol] > 0)) {
672  float tlog = log(tau_[D865][iSol]/tau_[D2250][iSol]) /
673  log(lut_.WAVE[D2250]/lut_.WAVE[D865]);
674  angstrom_exp_[1][iSol] = tlog;
675  }
676 
677 // Interpolate for large and small cases for total optical depth at 0.55
678  float albedo_R_small[NWAV][NUM_CASES_SMALL];
679  float albedo_R_big[NWAV][NUM_CASES_BIG];
680  float albedo_T_small[NWAV][NUM_CASES_SMALL];
681  float albedo_T_big[NWAV][NUM_CASES_BIG];
682  memset(albedo_R_small, 0, NWAV*NUM_CASES_SMALL*sizeof(float));
683  memset(albedo_R_big, 0, NWAV*NUM_CASES_BIG*sizeof(float));
684  memset(albedo_T_small, 0, NWAV*NUM_CASES_SMALL*sizeof(float));
685  memset(albedo_T_big, 0, NWAV*NUM_CASES_BIG*sizeof(float));
686 
687  for (int iWav=0; iWav<NWAV; iWav++) {
688  interp_extrap(NAOT,tau_[D550][iSol],lut_.TAUAS[D550][iSmall],
689  albedo_R_small_tau_[iWav][iSmall],albedo_R_small[iWav][iSmall]);
690  interp_extrap(NAOT,tau_[D550][iSol],lut_.TAUAS[D550][iSmall],
691  albedo_T_small_tau_[iWav][iSmall],albedo_T_small[iWav][iSmall]);
692  interp_extrap(NAOT,tau_[D550][iSol],lut_.TAUAB[D550][iBig],
693  albedo_R_big_tau_[iWav][iBig],albedo_R_big[iWav][iBig]);
694  interp_extrap(NAOT,tau_[D550][iSol],lut_.TAUAB[D550][iBig],
695  albedo_T_big_tau_[iWav][iBig],albedo_T_big[iWav][iBig]);
696  }
697 
698 // Compute CCN, asymmetry factor, backscattering ratio,
699 // effective radius, effective variance, and optical thicknesses for
700 // large and small mode.
701  float wout = 0.0;
702  interp_extrap(numWS, ws_, &ws_lut[iws1],
703  &lut_.EXTSMALL[D550][iSmall][iws1], wout);
704  float nSmall = tau_X55_small/wout;
705  interp_extrap(numWS, ws_, &ws_lut[iws1],
706  &lut_.EXTBIG[D550][iBig][iws1], wout);
707  float nBig = tau_X55_big/wout;
708  interp_extrap(numWS, ws_, &ws_lut[iws1],
709  &lut_.CCNSMALL[iSmall][iws1], wout);
710  ccn_[iSol] = nSmall* wout;
711 
712  eff_radius_[iSol] = 0.0;
713  eff_variance_[iSol] = 0.0;
714  mass_con_ocean_[iSol] = 0.0;
715 
716  if (nSmall > 0.0 ||nBig > 0.0 ) {
717  float m[NUM_SIZES][NUM_MOMENTS];
718  memset(m, 0, NUM_SIZES*NUM_MOMENTS*sizeof(float));
719  for ( int iM=1; iM < NUM_MOMENTS; iM++) {
720  interp_extrap(numWS, ws_, &ws_lut[iws1],
721  &lut_.MOMENTSSMALL[iSmall][iM][iws1], m[ISMALL][iM]);
722  interp_extrap(numWS, ws_, &ws_lut[iws1],
723  &lut_.MOMENTSBIG[iBig][iM][iws1], m[IBIG][iM]);
724  }
725  eff_radius_[iSol] = (nSmall*m[ISMALL][2]+ nBig*m[IBIG][2]) /
726  (nSmall*m[ISMALL][1]+ nBig*m[IBIG][1]);
727  mass_con_ocean_[iSol] = nSmall*m[ISMALL][2]+ nBig*m[IBIG][2];
728  eff_variance_[iSol] = (nSmall*m[ISMALL][3]+ nBig*m[IBIG][3])/
729  (pow(eff_radius_[iSol],2)*(nSmall*m[ISMALL][1]+nBig*4))-1.0;
730  }
731 
732  for (int iWav=0; iWav<NWAV; iWav++) {
733  float al[NUM_SIZES] = {0.0,0.0};
734  float bs[NUM_SIZES] = {0.0,0.0};
735  float as[NUM_SIZES] = {0.0,0.0};
736  interp_extrap(numWS, ws_, &ws_lut[iws1],
737  &lut_.ALBEDOSMALL[iWav][iSmall][iws1], al[ISMALL]);
738  interp_extrap(numWS, ws_, &ws_lut[iws1],
739  &lut_.ALBEDOBIG[iWav][iBig][iws1], al[IBIG]);
740  interp_extrap(numWS, ws_, &ws_lut[iws1],
741  &lut_.BACKSCTTSMALL[iWav][iSmall][iws1], bs[ISMALL]);
742  interp_extrap(numWS, ws_, &ws_lut[iws1],
743  &lut_.BACKSCTTBIG[iWav][iBig][iws1], bs[IBIG]);
744  interp_extrap(numWS, ws_, &ws_lut[iws1],
745  &lut_.ASSYMSMALL[iWav][iSmall][iws1], as[ISMALL]);
746  interp_extrap(numWS, ws_, &ws_lut[iws1],
747  &lut_.ASSYMBIG[iWav][iBig][iws1], as[IBIG]);
748 // compute backscattering ratio
749  if ((tau_small_[iWav][iSol] > 0.0) ||
750  (tau_big_[iWav][iSol] > 0.0)) {
751  backscatter_[iWav][iSol] =
752  (al[ISMALL]*tau_small_[iWav][iSol]*bs[ISMALL] +
753  al[IBIG]*tau_big_[iWav][iSol]*bs[IBIG]) /
754  (al[ISMALL]*tau_small_[iWav][iSol] +
755  al[IBIG]*tau_big_[iWav][iSol]);
756  assym_[iWav][iSol] =
757  (al[ISMALL]*tau_small_[iWav][iSol]*as[ISMALL] +
758  al[IBIG]*tau_big_[iWav][iSol]*as[IBIG]) /
759  (al[ISMALL]*tau_small_[iWav][iSol] +
760  al[IBIG]*tau_big_[iWav][iSol]);
761  refl_flux_[iWav][iSol] =
762  (tau_small_[iWav][iSol]* albedo_R_small[iWav][iSmall]
763  + tau_big_[iWav][iSol]*albedo_R_big[iWav][iBig])/
764  (tau_small_[iWav][iSol] + tau_big_[iWav][iSol]);
765 // subtract rayleigh and apply angle correction because it is divided in LUT
766  refl_flux_[iWav][iSol] =
767  (refl_flux_[iWav][iSol] - albedo_R_small_tau_[iWav][0][0]) *
768  cos(solz_*DEGtoRAD);
769 // compute transmittance
770  trans_flux_[iWav][iSol]=
771  (tau_small_[iWav][iSol]*albedo_T_small[iWav][iSmall]
772  + tau_big_[iWav][iSol]*albedo_T_big[iWav][iBig])/
773  (tau_small_[iWav][iSol]+tau_big_[iWav][iSol]);
774  }
775  }
776 
777  return status;
778 }
779 
780 /**************************************************************************
781  * NAME: compute_avg_refl()
782  *
783  * DESCRIPTION: This subroutine averages the reflectances for each
784  * 10*10 pixel square and finds standard deviation for
785  * averaged Reflectances
786  *
787  *************************************************************************/
788 
789 int DtAlgOcean::compute_avg_refl()
790 {
791  int status = DTDB_SUCCESS;
792 
793  float var[NWAV];
794  float dvar[NWAV];
795  float refl_interm[NWAV][GRIDX*GRIDY];
796  float array_interm[NWAV][GRIDX*GRIDY];
797  memset(var, 0, NWAV*sizeof(float));
798  memset(dvar, 0, NWAV*sizeof(float));
799  memset(refl_interm, 0, NWAV*GRIDX*GRIDY*sizeof(float));
800  memset(array_interm, 0, NWAV*GRIDX*GRIDY*sizeof(float));
801 
802  qcontrol_exclude_ = 0;
803  qcontrol_special_ = 0;
804  quality_dust_flag_glint_ = 0;
805  quality_dust_flag_off_glint_ = 0;
806  for (int iWav=0; iWav<NWAV; iWav++) {
807  good_pixels_[iWav] = 0;
808  }
809 
810  qcontrol_cirrus_ = 0;
811 
812  for (int iWav = 0; iWav<NWAV; iWav++) {
813  numData_[iWav]= 0;
814  for (int ix = 0; ix < GRIDX*GRIDY; ix++) {
815  refl_interm[iWav][ix] = DFILL_FLOAT;
816  }
817  }
818 
819  compute_glint_angle();
820  compute_scatter_angle(scatter_angle_);
821 
822  DtSedimentMask* sm = new DtSedimentMask(this);
823  short smask = 1;
824 
825  for ( size_t iy=0; iy<GRIDY; iy++ ) {
826  for ( size_t ix = 0; ix < GRIDX; ix++ ) {
827 
828  sm->compute(smask);
829 
830  if( rfld_[D1240] > 0) {
831  if ((cmask_ == 0) &&
832  (rfld_[D670] > 1.50*refl_rayl_[D670]) &&
833  ((rfld_[D1380] / rfld_[D1240] > 0.10) &&
834  (rfld_[D1380] / rfld_[D1240] < 0.30)) &&
835  ((rfld_[D1380] > 0.01) && (rfld_[D1380] <= 0.03))) {
836  qcontrol_cirrus_ = 1;
837  }
838  }
839  if (cmask_*smask > 0) {
840  refl_interm[D488][numData_[D488]] = rfld_[D488];
841  numData_[D488] += 1;
842  refl_interm[D550][numData_[D550]] = rfld_[D550];
843  numData_[D550] += 1;
844  refl_interm[D670][numData_[D670]] = rfld_[D670];
845  numData_[D670] += 1;
846  refl_interm[D865][numData_[D865]] = rfld_[D865];
847  numData_[D865] += 1;
848  refl_interm[D1240][numData_[D1240]] = rfld_[D1240];
849  numData_[D1240] += 1;
850  refl_interm[D1610][numData_[D1610]] = rfld_[D1610];
851  numData_[D1610] += 1;
852  refl_interm[D2250][numData_[D2250]] = rfld_[D2250];
853  numData_[D2250] += 1;
854  }
855  }
856  }
857 
858  delete sm;
859 
860 // Sort in ascending order and get the index for wavelength at 0.865
861  int index[GRIDX*GRIDY];
862 
863  sort_index(numData_[D865], &refl_interm[D865][0], index);
864 
865 // Reject 1/4 of brightest and darkest cloud-free pixels to eliminate
866 // the residual shadows and sub_cloud pixel based on wavelength of 865 nm
867  int n10 = numData_[D865]/4;
868  int n40 = (numData_[D865]-n10);
869 // if(n10 == 0) n10 = 1;
870 // if(n40 == 0) n40 = 1;
871  int q1_pixels = (int) LINE*LINE*0.025;
872  int q2_pixels = (int) LINE*LINE*0.050;
873  int q3_pixels = (int) LINE*LINE*0.075;
874 // Throw away 1/4 of brightest & darkest pixels in all wavelengths
875  if ((n40-n10) > q1_pixels) {
876  for (int iWav=0; iWav < NWAV; iWav++) {
877  for (int ix=n10; ix < n40; ix++) {
878  if ((refl_interm[iWav][index[ix]] > 0.0) &&
879  (refl_interm[iWav][index[ix]]<=1.3)) {
880  array_interm[iWav][good_pixels_[iWav]] =
881  refl_interm[iWav][index[ix]];
882  good_pixels_[iWav] += 1;
883  }
884  }
885  }
886  for (int iWav=0; iWav < NWAV; iWav++) {
887  if (good_pixels_[iWav] > 0) {
888  float refl_val = 0.0;
889  float std_val = 0.0;
890  mean_std(good_pixels_[iWav], array_interm[iWav],
891  refl_val, std_val);
892  refl_[iWav] = refl_val;
893  sdev_[iWav] = std_val;
894  var[iWav] = std_val/refl_val;
895  }
896  else {
897  refl_[iWav] = DFILL_FLOAT;
898  sdev_[iWav] = DFILL_FLOAT;
899  var[iWav] = DFILL_FLOAT;
900  }
901  }
902  dvar[D670] = -(1.0/(log(lut_.WAVE[D670]/lut_.WAVE[D865])))*
903  log((1.+ var[D670])/(1.+var[D865]));
904  dvar[D1240] = -(1./(log(lut_.WAVE[D1240]/lut_.WAVE[D865])))*
905  log((1.+ var[D1240])/(1.+var[D865]));
906  dvar[D1610] = -(1./(log(lut_.WAVE[D1610]/lut_.WAVE[D865])))*
907  log((1.+ var[D1610])/(1.+var[D865]));
908  dvar[D2250] = -(1./(log(lut_.WAVE[D2250]/lut_.WAVE[D865])))*
909  log((1.+ var[D2250])/(1.+var[D865]));
910 
911 // if GLINT_ANGLE <= glint threshold store the
912 // reflectances standard deviation and number of pixels only.
913  // Apply glint mask
914  float grthresh = GLINT_REFL_THRESHOLD/cos(solz_*DEGtoRAD);
915 
916  if ((glint_angle_ > 0.75*GLINT_ANGLE_THRESHOLD) &&
917  (glint_angle_ <= GLINT_ANGLE_THRESHOLD)) {
918  qcontrol_ = 10;
919  }
920  if( glint_refl_< grthresh &&
921  ((refl_[D488] > 0 && refl_[D670]> 0) &&
922  (refl_[D488]/refl_[D670] <= 0.75))) {
923  quality_dust_flag_off_glint_ = 1;
924  }
925  if( glint_refl_>= grthresh &&
926  ((refl_[D488] > 0 && refl_[D670] > 0) &&
927  (refl_[D488]/refl_[D670] <= 0.95))) {
928  quality_dust_flag_glint_ = 1;
929  }
930  if( glint_refl_ < grthresh ||
931  quality_dust_flag_off_glint_ == 1 ||
932  quality_dust_flag_glint_ == 1) {
933 // If there is valid data in 865 nm channel and at least one channel has valid data
934  int total_good = good_pixels_[D550]+good_pixels_[D670]+
935  good_pixels_[D1240]+good_pixels_[D1610]+good_pixels_[D2250];
936  if((good_pixels_[D865] > q1_pixels) &&
937  (total_good > q3_pixels)) {
938 // If Reflectance at wave 865 nm is greater than rhomin1
939 // the aerosol content is enough to process
940  if (refl_[D865] > (1.10*refl_rayl_[D865])) {
941 // Quality control flags are set for the number of cloud free pixels used
942 // qcontrol_=1 means that box is less than 10% cloud free.
943  numData_[D550] = good_pixels_[D865];
944  if ( good_pixels_[D865] < q2_pixels) {
945  qcontrol_=1;
946  }
947  if ( good_pixels_[D865] > q2_pixels) {
948  qcontrol_ = 0;
949  }
950 // If reflectance at wave 865 nm is less than rhomin2, then
951 // the aerosol content is not enough to derive size-distribution, but
952 // is enough to derive optical thickness ... process
953  if ( refl_[D865] < (1.50*refl_rayl_[D865])) {
954  qcontrol_= 3;
955  qcontrol_special_= 2;
956  }
957 // Test if channel 1. 24 is good to process
958  if (refl_[D1240]*cos(solz_*DEGtoRAD) < 3.600E-04) {
959  qcontrol_exclude_= 6;
960  }
961 // Test if channel 2.13 and 1.64 um are good to process
962  if (refl_[D1610]*cos(solz_*DEGtoRAD) < 3.600E-04) {
963  qcontrol_exclude_= 3;
964  }
965  if (refl_[D2250]*cos(solz_*DEGtoRAD) < 3.110E-04) {
966  qcontrol_exclude_= 4;
967  }
968  if (refl_[D2250]*cos(solz_*DEGtoRAD) < 3.110E-04 &&
969  refl_[D1610]*cos(solz_*DEGtoRAD) < 3.600E-04) {
970  qcontrol_exclude_= 5;
971  }
972 // Process only for valid data
973  qcontrol_ = qcontrol_exclude_;
974  if ( qcontrol_ > 0) {
975 // The aerosol type and content are variable
976  if (var[D865] > 0.05 &&
977  (fabs(dvar[D670]) > 0.15 ||
978  fabs(dvar[D1610]) > 0.15 ||
979  fabs(dvar[D2250]) > 0.15)) {
980  qcontrol_= 6;
981  }
982 // The aerosol content are variable not the aerosol content
983  if (var[D865] > 0.05 &&
984  (fabs(dvar[D670]) < 0.15 ||
985  fabs(dvar[D1610]) < 0.15 ||
986  fabs(dvar[D2250]) < 0.15)) {
987  qcontrol_= 7;
988  }
989  }
990  }
991  else {
992 // If Reflectance at wavelength 865 nm is less than rhomin1
993 // the aerosol content is not enough to process
994  qcontrol_= 17;
995  qcontrol_special_= 1;
996  }
997  }
998  else {
999 // If there is no valid data in all channels
1000  qcontrol_ = -20;
1001  qcontrol_special_ = -2;
1002  }
1003  }
1004  else {
1005 // Glint and only Reflectances, sdev and number of pixels will be stored
1006  qcontrol_special_= 3;
1007  numData_[D550] = good_pixels_[D865];
1008  qcontrol_= 11;
1009  }
1010  }
1011  else {
1012 // else for number of cloud-free pixels in box is too cloudy
1013 // If cloud free pixels are less than 10 no processing is performed .
1014  qcontrol_= -22;
1015  qcontrol_special_= -2;
1016  }
1017  if (qcontrol_ < 0) {
1018  refl_[D488] = DFILL_FLOAT;
1019  refl_[D550] = DFILL_FLOAT;
1020  refl_[D670] = DFILL_FLOAT;
1021  refl_[D865] = DFILL_FLOAT;
1022  refl_[D1240] = DFILL_FLOAT;
1023  refl_[D1610] = DFILL_FLOAT;
1024  refl_[D2250] = DFILL_FLOAT;
1025  }
1026 
1027  return status;
1028 }
1029 
1030 /**************************************************************************
1031  * NAME: store_results()
1032  *
1033  * This subroutine stores all the output variables to be written to output file.
1034  *
1035  ***************************************************************************/
1036 
1037 int DtAlgOcean::store_reflectance()
1038 {
1039  int status = DTDB_SUCCESS;
1040 
1041  float maxrefl = 1.2;
1042  float minrefl = 0.0;
1043  for (int iWav=0; iWav < NWAV; iWav++) {
1044  if (refl_[iWav] >= minrefl) {
1045  if (refl_[iWav] <= maxrefl) {
1046  sds_refl_[iWav] = refl_[iWav];
1047  sds_refl_sdev_[iWav] = sdev_[iWav];
1048  } else {
1049  sds_refl_[iWav] = maxrefl;
1050  sds_refl_sdev_[iWav] = sdev_[iWav];
1051  }
1052  } else {
1053  sds_refl_[iWav] = minrefl;
1054  sds_refl_sdev_[iWav] = 0.0;
1055  }
1056  if (good_pixels_[iWav] >= 0) {
1057  sds_numPixels_[iWav] = good_pixels_[iWav];
1058  } else {
1059  sds_numPixels_[iWav] = 0;
1060  }
1061  }
1062 
1063  return status;
1064 }
1065 
1066 /****************************************************************************
1067  * NAME: store_output()
1068  *
1069  * This subroutine stores all the output variables to be written to output file.
1070  *
1071  ***************************************************************************/
1072 
1073 int DtAlgOcean::store_output()
1074 {
1075  int status = DTDB_SUCCESS;
1076 
1077  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1078  if ((tau_avg_[IBEST][iWav] > -0.1) /* && (tau_avg_[IBEST][iWav] < 5.0) */) {
1079  sds_tau_best_[iWav] = tau_avg_[IBEST][iWav];
1080  }
1081  if ((tau_avg_[IAVG][iWav] > -0.1) /* && (tau_avg_[IBEST][iWav] < 5.0) */) {
1082  sds_tau_avg_[iWav] = tau_avg_[IAVG][iWav];
1083  }
1084  if (tau_small_avg_[IBEST][iWav] > 0.0) {
1085  sds_tau_small_best_[iWav] = tau_small_avg_[IBEST][iWav];
1086  }
1087  if (tau_small_avg_[IAVG][iWav] > 0.0) {
1088  sds_tau_small_avg_[iWav] = tau_small_avg_[IAVG][iWav];
1089  }
1090  if (tau_big_avg_[IBEST][iWav] > 0.0) {
1091  sds_tau_big_best_[iWav] = tau_big_avg_[IBEST][iWav];
1092  }
1093  if (tau_big_avg_[IAVG][iWav] > 0.0) {
1094  sds_tau_big_avg_[iWav] = tau_big_avg_[IAVG][iWav];
1095  }
1096  if (assym_avg_[IBEST][iWav] > 0.0) {
1097  sds_assy_best_[iWav] = assym_avg_[IBEST][iWav];
1098  }
1099  if (assym_avg_[IAVG][iWav] > 0.0) {
1100  sds_assy_avg_[iWav] = assym_avg_[IAVG][iWav];
1101  }
1102  if (backscatter_avg_[IBEST][iWav] > 0.0) {
1103  sds_back_best_[iWav] = backscatter_avg_[IBEST][iWav];
1104  }
1105  if (backscatter_avg_[IAVG][iWav] > 0.0) {
1106  sds_back_avg_[iWav] = backscatter_avg_[IAVG][iWav];
1107  }
1108  if (refl_flux_avg_[IBEST][iWav] > 0.0) {
1109  sds_reff_best_[iWav] = refl_flux_avg_[IBEST][iWav];
1110  }
1111  if (refl_flux_avg_[IAVG][iWav] > 0.0) {
1112  sds_reff_avg_[iWav] = refl_flux_avg_[IAVG][iWav];
1113  }
1114  if (trans_flux_avg_[IBEST][iWav] > 0.0) {
1115  sds_tranf_best_[iWav] = trans_flux_avg_[IBEST][iWav];
1116  }
1117  if (trans_flux_avg_[IAVG][iWav] > 0.0) {
1118  sds_tranf_avg_[iWav] = trans_flux_avg_[IAVG][iWav];
1119  }
1120  }
1121  for ( int iS=0; iS<NUM_STATS; iS++ ) {
1122  if (mass_con_ocean_avg_[iS] > 0.0) {
1123  sds_Mass_Conc_[iS] = mass_con_ocean_avg_[iS];
1124  }
1125  if (eff_radius_avg_[iS] > 0.0) {
1126  sds_EffRad_[iS] = eff_radius_avg_[iS];
1127  }
1128  if (ccn_avg_[iS] > 0.0) {
1129  sds_CCN_[iS] = ccn_avg_[iS];
1130  }
1131  if (angstrom_exp_avg_[iS][0] > -1.0) {
1132  sds_Angs_Coeff1_[iS] = angstrom_exp_avg_[iS][0];
1133  }
1134  if (angstrom_exp_avg_[iS][1] > -1.0) {
1135  sds_Angs_Coeff2_[iS] = angstrom_exp_avg_[iS][1];
1136  }
1137  if (funmin_[iS][0] >= 0.0) {
1138  sds_Least_Error_[iS] = funmin_avg_[iS];
1139  }
1140  if (tau_X55_[iS][0] >= 0.0) {
1141  sds_tau_X55_[iS] = tau_X55_avg_[iS];
1142  }
1143  if (xmin_[iS][0] >= 0.0) {
1144  sds_Small_Weighting_[iS] = xmin_avg_[iS];
1145  }
1146  if (solution_index_avg_[iS][0] >= 0.0) {
1147  sds_Sol_Index_Small_[iS] = (short) solution_index_avg_[iS][0];
1148  }
1149  if (solution_index_avg_[iS][1] >= 0.0) {
1150  sds_Sol_Index_Large_[iS] = (short) solution_index_avg_[iS][1];
1151  }
1152  }
1153  sds_land_ocean_quality_ = quality_to_pass_[0];
1154 
1155  return status;
1156 }
1157 
1158 /****************************************************************************
1159  * NAME: average_output()
1160  *
1161  * Subroutine sorts according to minimum error and averages all variables of output
1162  *
1163  ***************************************************************************/
1164 
1165 int DtAlgOcean::average_output()
1166 {
1167  int status = DTDB_SUCCESS;
1168 
1169  float sdev = 0.0;
1170  int index[NUMCASES*NUMCASEB];
1171  memset(index, 0, NUMCASES*NUMCASEB*sizeof(int));
1172 
1173  sort_index(NUMCASES*NUMCASEB, funmin_[0], index);
1174 
1175  for ( int iWav=0; iWav<NWAV; iWav++) {
1176  tau_avg_[IBEST][iWav] = tau_[iWav][index[0]];
1177  tau_small_avg_[IBEST][iWav] = tau_small_[iWav][index[0]];
1178  tau_big_avg_[IBEST][iWav] = tau_big_[iWav][index[0]];
1179  backscatter_avg_[IBEST][iWav] = backscatter_[iWav][index[0]];
1180  assym_avg_[IBEST][iWav] = assym_[iWav][index[0]];
1181  refl_flux_avg_[IBEST][iWav] = refl_flux_[iWav][index[0]];
1182  trans_flux_avg_[IBEST][iWav] = trans_flux_[iWav][index[0]];
1183  }
1184  angstrom_exp_avg_[IBEST][0] = angstrom_exp_[0][index[0]];
1185  angstrom_exp_avg_[IBEST][1] = angstrom_exp_[1][index[0]];
1186  solution_index_avg_[IBEST][ISMALL] = solution_index_[ISMALL][index[0]];
1187  solution_index_avg_[IBEST][IBIG] = solution_index_[IBIG][index[0]];
1188  ccn_avg_[IBEST] = ccn_[index[0]];
1189  eff_radius_avg_[IBEST] = eff_radius_[index[0]];
1190  eff_variance_avg_[IBEST] = eff_variance_[index[0]];
1191  mass_con_ocean_avg_[IBEST] = mass_con_ocean_[index[0]];
1192  xmin_avg_[IBEST] = xmin_[0][index[0]];
1193  funmin_avg_[IBEST] = funmin_[0][index[0]];
1194  tau_X55_avg_[IBEST] = tau_X55_[0][index[0]];
1195 
1196  float weight[NUMCASES*NUMCASEB];
1197  memset(weight, 0, NUMCASES*NUMCASEB*sizeof(float));
1198  int num_good = 0;
1199  for ( int iSol=0; iSol < NUMCASES*NUMCASEB; iSol++ ) {
1200  if(100.0*funmin_[0][index[iSol]] <= Threshold_LSQ_Error) {
1201  weight[index[iSol]] = 1.0;
1202  num_good++;
1203  }
1204  else weight[index[iSol]] = 0.0;
1205  }
1206  if(num_good < 2) {
1207  qcontrol_ = 8;
1208  for ( int iSol=0; iSol<NUMCASES*NUMCASEB; iSol++ ) {
1209  if (iSol < 3) {
1210  if(funmin_[0][index[iSol]] > 0.0) {
1211  weight[index[iSol]] =
1212  1.0/pow(100.0*funmin_[0][index[iSol]],2);
1213  }
1214  }
1215  else weight[index[iSol]] = 0.0;
1216  }
1217  }
1218  else qcontrol_ = 0;
1219 
1220  for ( int iWav=0; iWav<NWAV; iWav++) {
1221  mean_std_weighted( NUMCASES*NUMCASEB, tau_[iWav],
1222  tau_avg_[IAVG][iWav], sdev, weight );
1223  mean_std_weighted( NUMCASES*NUMCASEB, tau_small_[iWav],
1224  tau_small_avg_[IAVG][iWav], sdev, weight );
1225  mean_std_weighted( NUMCASES*NUMCASEB, tau_big_[iWav],
1226  tau_big_avg_[IAVG][iWav], sdev, weight );
1227  mean_std_weighted( NUMCASES*NUMCASEB, backscatter_[iWav],
1228  backscatter_avg_[IAVG][iWav], sdev, weight );
1229  mean_std_weighted( NUMCASES*NUMCASEB, assym_[iWav],
1230  assym_avg_[IAVG][iWav], sdev, weight );
1231  mean_std_weighted( NUMCASES*NUMCASEB, refl_flux_[iWav],
1232  refl_flux_avg_[IAVG][iWav], sdev, weight );
1233  mean_std_weighted( NUMCASES*NUMCASEB, trans_flux_[iWav],
1234  trans_flux_avg_[IAVG][iWav], sdev, weight );
1235  }
1236  mean_std_weighted( NUMCASES*NUMCASEB, angstrom_exp_[0],
1237  angstrom_exp_avg_[IAVG][0], sdev, weight );
1238  mean_std_weighted( NUMCASES*NUMCASEB, angstrom_exp_[1],
1239  angstrom_exp_avg_[IAVG][1], sdev, weight );
1240  mean_std_weighted( NUMCASES*NUMCASEB, ccn_,
1241  ccn_avg_[IAVG], sdev, weight );
1242  mean_std_weighted( NUMCASES*NUMCASEB, eff_radius_,
1243  eff_radius_avg_[IAVG], sdev, weight );
1244  mean_std_weighted( NUMCASES*NUMCASEB, eff_variance_,
1245  eff_variance_avg_[IAVG], sdev, weight );
1246  mean_std_weighted( NUMCASES*NUMCASEB, mass_con_ocean_,
1247  mass_con_ocean_avg_[IAVG], sdev, weight );
1248  mean_std_weighted( NUMCASES*NUMCASEB, xmin_[0],
1249  xmin_avg_[IAVG], sdev, weight );
1250  mean_std_weighted( NUMCASES*NUMCASEB, funmin_[0],
1251  funmin_avg_[IAVG], sdev, weight );
1252  mean_std_weighted( NUMCASES*NUMCASEB, tau_X55_[0],
1253  tau_X55_avg_[IAVG], sdev, weight );
1254  solution_index_avg_[IAVG][ISMALL] = solution_index_[ISMALL][index[0]];
1255  solution_index_avg_[IAVG][IBIG] = solution_index_[IBIG][index[0]];
1256 
1257  return status;
1258 }
1259 
1260 /***************************************************************************
1261  * NAME: compute_minimum()
1262  *
1263  * This function computes a minimum value for aerosol retrieval algorithm.
1264  * It finds the minimum using a direct least mean square calculation.
1265  *
1266  ***************************************************************************/
1267 
1268 int DtAlgOcean::compute_minimum(int iBig, int iSmall, int iSol)
1269 {
1270  int status = DTDB_SUCCESS;
1271 
1272  bool bUseWav[NWAV];
1273  // For summation do not use 1st wavelength
1274  bUseWav[D488] = false;
1275  bUseWav[D550] = true;
1276  bUseWav[D670] = true;
1277  bUseWav[D865] = true;
1278  // if wavelength 1.24 or 1.64 or 2.13 is small do not use for inversion
1279  bUseWav[D1240] = (qcontrol_exclude_== 6) ? false : true;
1280  bUseWav[D1610] = (qcontrol_exclude_== 3 ||
1281  qcontrol_exclude_== 5) ? false : true;
1282  bUseWav[D2250] = (qcontrol_exclude_== 4 ||
1283  qcontrol_exclude_== 5) ? false : true;
1284 
1285  // Compute for tau using radiance value of 0.86 um and tau values of wav550
1286  float tau_X55[2] = {0.0,0.0};
1287  interp_extrap(NAOT,refl_[D865],refl_small_[D865][iSmall],
1288  lut_.TAUAS[D550][0],tau_X55[ISMALL]);
1289  interp_extrap(NAOT,refl_[D865],refl_big_[D865][iBig],
1290  lut_.TAUAB[D550][0],tau_X55[IBIG]);
1291 
1292  float rbdif[NWAV];
1293  float sbdif[NWAV];
1294  float trefl[NWAV][2];
1295  float scale[NWAV];
1296  memset(rbdif, 0, NWAV*sizeof(float));
1297  memset(sbdif, 0, NWAV*sizeof(float));
1298  memset(trefl, 0, NWAV*2*sizeof(float));
1299  memset(scale, 0, NWAV*sizeof(float));
1300  float rb2 = 0;
1301  float sb2 = 0;
1302  float rbsb = 0;
1303  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1304  float denom = refl_[iWav]-refl_rayl_[iWav] + 0.01;
1305  denom = (denom < 0.01) ? 0.01 : denom;
1306  scale[iWav] = 1.0/denom/denom;
1307  interp_extrap(NAOT,tau_X55[ISMALL],lut_.TAUAS[D550][0],
1308  refl_small_[iWav][iSmall],trefl[iWav][ISMALL]);
1309  interp_extrap(NAOT,tau_X55[IBIG],lut_.TAUAB[D550][0],
1310  refl_big_[iWav][iBig],trefl[iWav][IBIG]);
1311  if (bUseWav[iWav]) {
1312  rbdif[iWav] = refl_[iWav] - trefl[iWav][IBIG];
1313  sbdif[iWav] = trefl[iWav][ISMALL] - trefl[iWav][IBIG];
1314  rb2 += rbdif[iWav]*rbdif[iWav]*scale[iWav];
1315  sb2 += sbdif[iWav]*sbdif[iWav]*scale[iWav];
1316  rbsb += rbdif[iWav]*sbdif[iWav]*scale[iWav];
1317  }
1318  }
1319 
1320  float xm = rbsb / sb2;
1321  xm = (xm > 1) ? 1.0 : xm;
1322  xm = (xm < 0) ? 0.0 : xm;
1323  float asum = 0.0;
1324  float sum_good_pixels = 0.0;
1325  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1326  float mrfl = (xm*trefl[iWav][ISMALL]+(1.0-xm)*trefl[iWav][IBIG]);
1327  error_[0][iWav][iSol] = (refl_[iWav] - mrfl) / refl_[iWav];
1328  if (bUseWav[iWav]) {
1329  asum += error_[0][iWav][iSol]*error_[0][iWav][iSol];
1330  sum_good_pixels += good_pixels_[iWav];
1331  }
1332  }
1333  xmin_[0][iSol] = xm;
1334  funmin_[0][iSol] = asum/(NWAV-2);
1335  tau_X55_[0][iSol] = xm*tau_X55[ISMALL] + (1-xm)*tau_X55[IBIG];
1336 
1337  return status;
1338 }
1339 
1340 /***************************************************************************
1341  * NAME: compute_minimum_baseline()
1342  *
1343  * This function computes a minimum value for aerosol retrieval algorithm.
1344  * It finds the minimum of a function by using an interval halving search.
1345  *
1346  ***************************************************************************/
1347 
1348 int DtAlgOcean::compute_minimum_baseline(int iBig, int iSmall, int iSol)
1349 {
1350  int status = DTDB_SUCCESS;
1351 
1352  float iFlag = 0;
1353  float a[2][5] = {{0.0, 0.0, 0.0, 0.0, 0.0},{0.0, 0.0, 0.0, 0.0, 0.0}};
1354  float x[5] = {0.0, 0.5, 1.0, 0.0, 0.0};
1355  float f[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
1356  f[0] = fun_tau(x[0], iBig, iSmall, iSol);
1357  f[1] = fun_tau(x[1], iBig, iSmall, iSol);
1358  f[2] = fun_tau(x[2], iBig, iSmall, iSol);
1359 
1360  for ( int i=0; i<10; i++ ) {
1361  x[3] = (x[0]+x[1])/2.0;
1362  x[4] = (x[1]+x[2])/2.0;
1363  f[3] = fun_tau(x[3], iBig, iSmall, iSol);
1364  f[4] = fun_tau(x[4], iBig, iSmall, iSol);
1365  std::copy(x, x+5, a[0]);
1366  std::copy(f, f+5, a[1]);
1367 
1368  sort_inplace(5,a[1],a[0]);
1369 
1370  float rmin = a[0][0];
1371  if( rmin <= x[3]) {
1372  iFlag++;
1373  x[2] = x[1];
1374  x[1] = x[3];
1375  f[2] = f[1];
1376  f[1] = f[3];
1377  }
1378  if((fabs(rmin-x[1]) < fabs(x[1])*0.000001) && (rmin > x[3])) {
1379  iFlag++;
1380  x[0] = x[3];
1381  x[2] = x[4];
1382  f[0] = f[3];
1383  f[2] = f[4];
1384  }
1385  if((rmin >= x[4]) &&
1386  (fabs(rmin-x[1]) > fabs(x[1])*0.000001) && (rmin > x[3])) {
1387  iFlag++;
1388  x[0] = x[1];
1389  x[1] = x[4];
1390  f[0] = f[1];
1391  f[1] = f[4];
1392  }
1393  }
1394  if( iFlag > 1) {
1395 
1396  sort_inplace(3, f, x);
1397 
1398  xmin_[1][iSol] = x[0];
1399  funmin_[1][iSol] = f[0];
1400  }
1401  else {
1402  xmin_[1][iSol] = DFILL_FLOAT;
1403  funmin_[1][iSol] = DFILL_FLOAT;
1404  }
1405 
1406  return status;
1407 }
1408 
1409 /***************************************************************************
1410  * NAME: fun_tau()
1411  *
1412  * DESCFRIPTION: Function to be minimized
1413  *
1414  ***************************************************************************/
1415 
1416 float DtAlgOcean::fun_tau( float xmin, int iBig, int iSmall, int iSol )
1417 {
1418  float xs[NWAV][NAOT];
1419  float ys[NWAV][NAOT];
1420  memset(xs, 0, NWAV*NAOT*sizeof(float));
1421  memset(ys, 0, NWAV*NAOT*sizeof(float));
1422 // Compute function to be minimized.
1423  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1424  for ( int iTau=0; iTau<NAOT; iTau++ ) {
1425  xs[iWav][iTau] = lut_.TAUAS[D550][0][iTau];
1426  ys[iWav][iTau]= xmin*refl_small_[iWav][iSmall][iTau] +
1427  (1.-xmin)*refl_big_[iWav][iBig][iTau];
1428  }
1429  }
1430 
1431 // Compute for tau using radiance value of 0.86 um and tau values of wav550
1432  interp_extrap(NAOT,refl_[D865],ys[D865],xs[D550],tau_X55_[1][iSol]);
1433 
1434 // For TAUX55 compute reflectance for all wavelengths.
1435  float alxwav[NWAV];
1436  memset(alxwav, 0, NWAV*sizeof(float));
1437  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1438  interp_extrap(NAOT,tau_X55_[1][iSol],xs[D550],ys[iWav],alxwav[iWav]);
1439  error_[1][iWav][iSol] = refl_[iWav] - alxwav[iWav];
1440  }
1441 
1442  bool bUseWav[NWAV];
1443 // For summation do not use 1st wavelength
1444  bUseWav[D488] = false;
1445  bUseWav[D550] = true;
1446  bUseWav[D670] = true;
1447  bUseWav[D865] = true;
1448 // if wavelength 1.24 or 1.64 or 2.13 is small do not use for inversion
1449  bUseWav[D1240] = (qcontrol_exclude_== 6) ? false : true;
1450  bUseWav[D1610] = (qcontrol_exclude_== 3 ||
1451  qcontrol_exclude_== 5) ? false : true;
1452  bUseWav[D2250] = (qcontrol_exclude_== 4 ||
1453  qcontrol_exclude_== 5) ? false : true;
1454 
1455  float asum = 0.0;
1456  int sum_good_pixels = 0;
1457  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1458  if (bUseWav[iWav]) {
1459  float denom = (refl_[iWav]-refl_rayl_[iWav]) + 0.01;
1460  if( denom < 0.01) {
1461  denom = 0.01;
1462  }
1463  asum += pow(error_[1][iWav][iSol]/denom,2.0)*good_pixels_[iWav];
1464  sum_good_pixels += good_pixels_[iWav];
1465  }
1466  }
1467  float result = sqrt(asum/((float)sum_good_pixels));
1468 
1469  return result;
1470 }
1471 
1472 /***************************************************************************
1473  * NAME: assign_quality()
1474  *
1475  * DESCFRIPTION: This subroutine assigns quality flags and attributes.
1476  * Product Quality & Retrieval Processing QA flags over ocean
1477  * Quality of 3 and 4 (average solution) are repeated tion have
1478  * as best and average solution have same quality.
1479  *
1480  * quality_flag_[1]=0 NO retrieval ( NOT useful )
1481  * quality_flag_[1]=1 retrieval ( useful)
1482  *
1483  * For all non retrieval boxes quality_flag_ is set to zero to indicate
1484  * not useful and quality_flag_[4] is set to values of qcontrol_ to
1485  * indicate the cause of Retrieval. quality_flag_[5] is set to 11 for
1486  * no retrieval.
1487  *
1488  * Estimated quality of aerosol parameters
1489  * quality_flag_[1]=3 very Good
1490  * quality_flag_[1]=2 Good
1491  * quality_flag_[1]=1 Marginal
1492  * quality_flag_[1]=0 Bad
1493  * quality_flag_[3]=3 very Good
1494  * quality_flag_[3]=2 Good
1495  * quality_flag_[3]=1 Marginal
1496  * quality_flag_[3]=0 Bad
1497  *
1498  ***************************************************************************/
1499 
1500 int DtAlgOcean::assign_quality()
1501 {
1502  int status = DTDB_SUCCESS;
1503 
1504  quality_to_pass_[0] = 0;
1505  quality_to_pass_[1] = 0;
1506  if ( qcontrol_ < 0) {
1507  quality_flag_[0] = 0;
1508  quality_flag_[1] = 0;
1509  quality_flag_[2] = 0;
1510  quality_flag_[3] = 0;
1511  quality_flag_[4] = abs(qcontrol_);
1512  quality_flag_[5] = 15;
1513  }
1514  else {
1515  quality_flag_[0] = 1;
1516  quality_flag_[2] = 1;
1517  if( quality_dust_flag_glint_ == 1) {
1518  qcontrol_= 12;
1519  }
1520  if( qcontrol_cirrus_== 1) {
1521  qcontrol_= 13;
1522  }
1523  if( quality_dust_flag_off_glint_ == 1) {
1524  qcontrol_= 14;
1525  }
1526 // if qcontrol_ is 0 ( see doc.) then quality of retrieval is very good
1527  if( qcontrol_ == 0 ) {
1528  quality_flag_[1] = 3;
1529  quality_flag_[3] = 3;
1530  }
1531 // if qcontrol_ is 7 or 14 ( see doc.) then quality of retrieval is good
1532  if( qcontrol_ == 7 || qcontrol_ == 14) {
1533  quality_flag_[1] = 2;
1534  quality_flag_[3] = 2;
1535  }
1536 // if qcontrol_ is 1,3,4,6,8 or 10 then quality of retrieval is Average
1537  if( qcontrol_== 1 || qcontrol_ == 3 || qcontrol_ == 4
1538  || qcontrol_== 6 || qcontrol_ == 8 || qcontrol_== 10) {
1539  quality_flag_[1] = 1;
1540  quality_flag_[3] = 1;
1541  }
1542 // if qcontrol_ is 2,5,9,12,13 quality of retrieval is poor
1543  if( qcontrol_== 2 || qcontrol_== 5 || qcontrol_ == 9
1544  || qcontrol_ == 12 || qcontrol_== 13 ) {
1545  quality_flag_[1] = 0;
1546  quality_flag_[3] = 0;
1547  }
1548  quality_flag_[4] = 0;
1549  quality_flag_[5] = qcontrol_;
1550  if(qcontrol_ == 17) {
1551  quality_flag_[1] = 3;
1552  quality_flag_[3] = 3;
1553  }
1554  }
1555  quality_to_pass_[0] = quality_flag_[1];
1556  quality_to_pass_[1] = qcontrol_;
1557 
1558  return status;
1559 }
1560 
1561 /***************************************************************************
1562  * NAME: fill_values()
1563  *
1564  * DESCFRIPTION: This subroutine assigns fill values to output data
1565  * where appropriate.
1566  *
1567  ***************************************************************************/
1568 
1569 int DtAlgOcean::set_fill_out()
1570 {
1571  int status = DTDB_SUCCESS;
1572 
1573  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1574  sds_refl_[iWav] = DFILL_FLOAT;
1575  sds_refl_sdev_[iWav] = DFILL_FLOAT;
1576  sds_numPixels_[iWav] = DFILL_SHORT;
1577  }
1578  for ( int iWav=0; iWav<NWAV; iWav++ ) {
1579  sds_tau_best_[iWav] = DFILL_FLOAT;
1580  sds_tau_avg_[iWav] = DFILL_FLOAT;
1581  sds_tau_small_best_[iWav] = DFILL_FLOAT;
1582  sds_tau_small_avg_[iWav] = DFILL_FLOAT;
1583  sds_tau_big_best_[iWav] = DFILL_FLOAT;
1584  sds_tau_big_avg_[iWav] = DFILL_FLOAT;
1585  sds_assy_best_[iWav] = DFILL_FLOAT;
1586  sds_assy_avg_[iWav] = DFILL_FLOAT;
1587  sds_back_best_[iWav] = DFILL_FLOAT;
1588  sds_back_avg_[iWav] = DFILL_FLOAT;
1589  sds_reff_best_[iWav] = DFILL_FLOAT;
1590  sds_reff_avg_[iWav] = DFILL_FLOAT;
1591  sds_tranf_best_[iWav] = DFILL_FLOAT;
1592  sds_tranf_avg_[iWav] = DFILL_FLOAT;
1593  }
1594 
1595  for ( int iCase=0; iCase<NUM_SIZES; iCase++ ) {
1596  sds_Mass_Conc_[iCase] = DFILL_FLOAT;
1597  sds_EffRad_[iCase] = DFILL_FLOAT;
1598  sds_CCN_[iCase] = DFILL_FLOAT;;
1599  sds_Angs_Coeff1_[iCase] = DFILL_FLOAT;
1600  sds_Angs_Coeff2_[iCase] = DFILL_FLOAT;
1601  sds_Least_Error_[iCase] = DFILL_FLOAT;
1602  sds_Small_Weighting_[iCase] = DFILL_FLOAT;
1603  sds_Sol_Index_Small_[iCase] = DFILL_FLOAT;
1604  sds_Sol_Index_Large_[iCase] = DFILL_FLOAT;
1605  sds_tau_X55_[iCase] = DFILL_FLOAT;
1606  }
1607 
1608  return status;
1609 }
D412
@ D412
Definition: DtAlgorithm.h:24
index
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
DtAlgOcean::set_fill_out
int set_fill_out()
Definition: DtAlgOcean.cpp:1569
D865
@ D865
Definition: DtAlgorithm.h:20
status
int status
Definition: l1_czcs_hdf.c:32
DFILL_UBYTE
constexpr unsigned char DFILL_UBYTE
Definition: DDataset.hpp:28
DtAlgOcean::average_output
int average_output()
Definition: DtAlgOcean.cpp:1165
DtAlgOcean::compute_avg_refl
int compute_avg_refl()
Definition: DtAlgOcean.cpp:789
denom
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:99
D2250
@ D2250
Definition: DtAlgorithm.h:23
dtland.IBIG
int IBIG
Definition: dtland.py:33
DtAlgOcean::compute_minimum_baseline
int compute_minimum_baseline(int iBig, int iSmall, int iSol)
Definition: DtAlgOcean.cpp:1348
flags::BOWTIEDEL
@ BOWTIEDEL
DtAlgOcean::process
map< string, ddata * > process(vector< size_t > start, vector< size_t > count, map< string, ddata * > imap)
Definition: DtAlgOcean.cpp:261
DtAlgOcean::interpolate_angles
int interpolate_angles()
Definition: DtAlgOcean.cpp:461
as
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat as
Definition: HISTORY.txt:297
dtland.ISMALL
int ISMALL
Definition: dtland.py:32
DtAlgOcean::compute_tau_flux
int compute_tau_flux(int iBig, int iSmall, int iSol)
Definition: DtAlgOcean.cpp:611
D550
@ D550
Definition: DtAlgorithm.h:18
D670
@ D670
Definition: DtAlgorithm.h:19
DtSedimentMask::compute
int compute(short &mask)
Definition: DtMask.cpp:449
DtAlgOcean::store_output
int store_output()
Definition: DtAlgOcean.cpp:1073
f
double precision function f(R1)
Definition: tmd.lp.f:1454
DtAlgOcean::fun_tau
float fun_tau(float xmin, int iBig, int iSmall, int iSol)
Definition: DtAlgOcean.cpp:1416
hico.value_locate.x
list x
Definition: value_locate.py:64
DtAlgOcean::assign_quality
int assign_quality()
Definition: DtAlgOcean.cpp:1500
color_dtdb.scale
int scale
Definition: color_dtdb.py:459
M_PI
#define M_PI
Definition: pml_iop.h:15
bs
float * bs[MAX_BANDS]
Definition: pml_iop_tables.c:27
DtAlgOcean::compute_minimum
int compute_minimum(int iBig, int iSmall, int iSol)
Definition: DtAlgOcean.cpp:1268
D488
@ D488
Definition: DtAlgorithm.h:17
color_dtdb.result
result
Definition: color_dtdb.py:197
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
DFILL_FLOAT
constexpr float DFILL_FLOAT
Definition: DDataset.hpp:25
DtAlgOcean::store_reflectance
int store_reflectance()
Definition: DtAlgOcean.cpp:1037
DtAlgOcean::index_geometry
int index_geometry(float sza, float azim, float phi)
Definition: DtAlgOcean.cpp:383
fabs
#define fabs(a)
Definition: misc.h:93
copy
void copy(double **aout, double **ain, int n)
D1610
@ D1610
Definition: DtAlgorithm.h:22
DtAlgOcean::interpolate_rayleigh
int interpolate_rayleigh()
Definition: DtAlgOcean.cpp:406
MDN.benchmarks.utils.set_outputs
def set_outputs(output_keys)
Definition: utils.py:164
color_dtdb.ys
ys
Definition: color_dtdb.py:454
DtLutNetcdf::read_ocean_aerosol_lut
int read_ocean_aerosol_lut(dtOceanAerosolLUT &lo_lut)
Definition: DtLutNetcdf.cpp:1852
abs
#define abs(a)
Definition: misc.h:90
D1240
@ D1240
Definition: DtAlgorithm.h:21
i
int i
Definition: decode_rs.h:71
DtAlgorithm::initialize
virtual int initialize(map< string, ddata * > imap)
Definition: DtAlgorithm.cpp:76
a
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
DFILL_SHORT
constexpr short DFILL_SHORT
Definition: DDataset.hpp:27
viirs_ler_luts::sza
real sza
Definition: viirs_ler_luts.f95:46
DtAlgOcean::initialize
int initialize(map< string, ddata * > imap)
Definition: DtAlgOcean.cpp:237
D1380
@ D1380
Definition: DtAlgorithm.h:25
count
int count
Definition: decode_rs.h:79


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster