A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
calc_par.c
Go to the documentation of this file.
1 /******************************************************************************
2  * This routine computes daily Photosynthetically Active Radiation (PAR) just
3  * above the surface.
4  *
5  * The atmosphere and surface are modeled as a 2-layer system. The first layer
6  * contains molecules and aerosols and is located above the second layer, the
7  * cloud/surface layer. The TOA radiances are transformed into reflectance
8  * and corrected for gaseous absorption. The resulting reflectances are
9  * further corrected for molecular and aerosol scattering and combined to
10  * yield the instantaneous albedo of the cloud/surface layer in the PAR
11  * interval, A. This albedo is then used to compute instantaneous PAR. Daily
12  * PAR is finally obtained by integrating instantaneous PAR from sunrise to
13  * sunset, taking into account statistically diurnal effects on A.
14  *
15  * *** Important: This subroutine uses an external file containing
16  * aerosol optical properties.
17  *
18  * Output:
19  * PAR: Photosynthetically Active Radiation (mW/cm^2/um)
20  *
21  * Key changes in v 2.0:
22  * * wavelength integration weighting function added
23  *
24  * Key changes in v 1.8:
25  * * AsatMax renamed to AcldMax to indicate "cloud" emphasis
26  * * Use (Asat-AsBar) rather than Asat in the test that determines if
27  * Asat is too large
28  * * Added 670 section to aerosols table since we need that to choose
29  * the right model in GetPhase subroutine, + modified code in that
30  * routine
31  * * Changed TauCons from 20 to 15
32  *
33  * Key changes in v 1.7:
34  * * Corrected Tau term in spherical albedo expressions
35  * * Modified initial Asat definition with AsBar
36  * * Put 'cloudy' test after computation of new Asat
37  * * Replaced archaic intrinsic routines (COSD, SIND, ACOSD, JMOD, LONG)
38  *
39  * Key changes in v 1.6:
40  * * Angstrom exponents are expected to be primarily positive,
41  * that is, in the range (-0.1, 2.0) ... initial test was removed.
42  * * Surface albedo tables were added, with interpolation subroutines
43  * * TauAer(500nm) calculation was added (for Surface albedo tables)
44  * * Replaced AsBar calculation with call to As interpolation subroutine
45  * * Added test to estimate if pixel is clear or cloudy
46  * * Replaced AsAvg calculation with call to As interpolation subroutine
47  * * Removed AsBar and AsAvg from calculation of Aavg_spectral, in integration
48  * * Added Abar_spectral vs. AsAvg test for "function" integrand term
49  *
50  * Key changes in v 1.5:
51  * * addition of water vapor term in gaseous transmittance in integration
52  * * comparison test for Abar_spectral vs. AsAvg
53  *
54  * Authors: Robert Frouin <RFrouin@ucsd.edu>, scientific algorithms
55  * John McPherson <JMcPherson@ucsd.edu>, program structure
56  *
57  * Adapted from the original FORTRAN to C and better integrated with l2gen
58  * by Robert Lossing and Sean Bailey
59  *
60  * References:
61  *
62  * Frouin, R., and B. Chertock, 1992: A technique for global monitoring of net
63  * solar irradiance at the ocean surface. Part I: Model. J. Appl. Meteo., 31,
64  * 1056-1066.
65  *
66  * Frouin, R., D. W. Ligner, and C. Gautier, 1989: A Simple analytical formula
67  * to compute clear sky total and photosynthetically available solar irradiance
68  * at the ocean surface. J. Geophys. Res., 94, 9731-9742.
69  *
70  * Tanre, D., M. Herman, P.-Y. Deschamps, and A. De Leffe, 1979: Atmospheric
71  * modeling for Space measurements of ground reflectances, including
72  * bi-directional properties. Appl. Optics, 18, 21,3587-21,3597.
73  *
74  * Young, D. F., P. Minnis, D. R. Doelling, G. G. Gibson, and T. Wong:
75  * Temporal interpolation methods for the Clouds and the Earth's Radiant
76  * Energy System (CERES) Experiment. J. Appl. Meteo., 37, 572-590.
77  *
78  *****************************************************************************/
79 #include <math.h>
80 #include <stdio.h>
81 #include <stdlib.h>
82 #include <timeutils.h>
83 #include "genutils.h"
84 #include "l12_proto.h"
85 #include "l2_struc.h"
86 #include "par_utils.h"
87 
88 static int ini = 0;
89 luts_par luts_data;
90 float *grid_ozone[2];
91 float *grid_watvap[2];
92 float *grid_tg[4], *grid_td[4], *grid_rho[6];
93 float *grid_scalar_par[3];
94 float *grid_scalar_inst_par[3];
95 static int16_t year, month, mday, doy;
96 static double sec;
97 static size_t t_step;
98 static float observed_time;
99 static float *t_range;
100 const float A = 8.435e-03;
101 const float B = -1.225e-04;
102 const float C = 1.40e-04;
103 float Ps0 = EARTH_SURF_PRESSURE;
104 const float delta = 0.0095;
105 const float Gamma0 = 1.f / 3.f;
106 const size_t N_5nm = 60;
107 const float F0_5nm[] = {163.33, 167.94, 174.72, 175.26, 170.75, 159.98, 163.31, 177.76, 186.23, 199.60,
108  204.19, 203.33, 202.92, 200.22, 200.56, 204.09, 197.04, 190.83, 195.62, 193.83,
109  191.92, 193.99, 188.39, 181.65, 183.47, 188.85, 190.71, 187.82, 186.81, 187.72,
110  187.50, 185.28, 183.95, 184.37, 185.20, 185.36, 184.33, 180.52, 178.12, 177.98,
111  176.71, 175.58, 171.55, 170.52, 169.82, 167.42, 166.57, 165.00, 163.21, 160.81,
112  155.44, 153.10, 155.10, 154.64, 152.72, 150.99, 148.61, 146.63, 145.06, 142.50};
113 
114 const float wl_5nm[] = {402.5, 407.5, 412.5, 417.5, 422.5, 427.5, 432.5, 437.5, 442.5, 447.5, 452.5, 457.5,
115  462.5, 467.5, 472.5, 477.5, 482.5, 487.5, 492.5, 497.5, 502.5, 507.5, 512.5, 517.5,
116  522.5, 527.5, 532.5, 537.5, 542.5, 547.5, 552.5, 557.5, 562.5, 567.5, 572.5, 577.5,
117  582.5, 587.5, 592.5, 597.5, 602.5, 607.5, 612.5, 617.5, 622.5, 627.5, 632.5, 637.5,
118  642.5, 647.5, 652.5, 657.5, 662.5, 667.5, 672.5, 677.5, 682.5, 687.5, 692.5, 697.5};
119 static float weight_5nm[60];
120 static float Denom_5nm = 0;
121 
122 
123 void calc_scalar_inst_par(l2str *l2rec, int ip, float par_above_ins, float *par_scalar_ins) {
124  l1str *l1rec = l2rec->l1rec;
125  size_t dim_cot = luts_data.scalar_inst_luts.dim_cot;
126  const float *cot_grid = luts_data.scalar_inst_luts.cot;
127  float cf_planar_grid[dim_cot];
128  float cf_scalar_grid[dim_cot];
129  float second_derivative_cf_scalar_par[dim_cot];
130 
131  size_t dims[] = {
132  luts_data.scalar_inst_luts.dim_wind_speed,
133  luts_data.scalar_inst_luts.dim_solz,
134  dim_cot,
135  };
136  for (size_t i = 0; i < dim_cot; i++) {
137  float point[3] = {l1rec->ws[ip], l1rec->csolz[ip], cot_grid[i]};
138  cf_planar_grid[i] = interp3d(dims, point, grid_scalar_inst_par, luts_data.scalar_inst_luts.cf_pard_p);
139  cf_scalar_grid[i] = interp3d(dims, point, grid_scalar_inst_par, luts_data.scalar_inst_luts.cf_pard_m);
140  }
141  float point[2] = {l1rec->ws[ip], l1rec->csolz[ip]};
142 
143  spline(cf_planar_grid, cf_scalar_grid, dim_cot, 0, 0, second_derivative_cf_scalar_par);
144  // need to define interp2d
145  float inst_PARd_c = interpnd(2, dims, point, grid_scalar_inst_par, luts_data.scalar_inst_luts.pard_p_cs);
146  float inst_PARo_c = interpnd(2, dims, point, grid_scalar_inst_par, luts_data.scalar_inst_luts.pard_m_cs);
147  float inst_PARo_o = interpnd(2, dims, point, grid_scalar_inst_par, luts_data.scalar_inst_luts.pard_m_oc);
148  float inst_cf_obs = par_above_ins / inst_PARd_c / EINSTEIN * 3600 * 24;
149  float inst_cf_obs_spline;
150  splint(cf_planar_grid, cf_scalar_grid, second_derivative_cf_scalar_par, dim_cot, inst_cf_obs, &inst_cf_obs_spline);
151  // correct for 2d splie behavior when cloudy
152  if(inst_cf_obs < cf_planar_grid[1] && inst_cf_obs_spline < inst_cf_obs){
153  inst_cf_obs_spline = inst_cf_obs;
154  }
155  float inst_cf = inst_PARo_o/inst_PARo_c;
156  float slope = (inst_PARo_c - inst_PARo_o)/(1-inst_cf);
157  *par_scalar_ins = (slope * (inst_cf_obs_spline - inst_cf) + inst_PARo_o)/3600/24;
158 }
159 
160 void calc_scalar_par_mean_cosine(l2str *l2rec, int ip, float par_above, float par_c, float *scalar_par,
161  float *mean_cosine) {
162  if (par_c == BAD_FLT) {
163  *scalar_par = BAD_FLT;
164  *mean_cosine = BAD_FLT;
165  return;
166  }
167 
168  l1str *l1rec = l2rec->l1rec;
169  float point[] = {l1rec->ws[ip], (float) doy, l1rec->lat[ip]};
170  size_t dims[] = {
171  luts_data.scalar_luts.dim_wind_speed,
172  luts_data.scalar_luts.dim_doy,
173  luts_data.scalar_luts.dim_latitude,
174  };
175  float PARo_1 = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.PARo);
176  float PARo_2 = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.PARo_c);
177  float PARc_above = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.PARc_above);
178  float cf = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.CF);
179  float S1 = (PARo_2 - PARo_1) / (1. - cf);
180  const float cf_obs = par_above / PARc_above;
181  float PARo_est = S1 * (cf_obs - cf) + PARo_1;
182  *scalar_par = PARo_est;
183 
184  float mu1 = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.mu_cosine);
185  float mu2 = interp3d(dims, point, grid_scalar_par, luts_data.scalar_luts.mu_cosine_c);
186 
187  S1 = (mu2 - mu1) / (1. - cf);
188  float mu_est = S1 * (cf_obs - cf) + mu1;
189  *mean_cosine = mu_est;
190 }
191 
192 float calc_par(l2str *l2rec, int ip, int nbands, float *Lt, float taua, float angstrom, float *wl, float *fo,
193  float *ko3, float *taumolbar) {
194  int i, band;
195  float r[nbands], csolz, transg;
196  float weight[nbands];
197  float rp[nbands], transa, sa, ps0, del;
198  float asymfac, beta, tau[nbands], taumol[nbands], tauaer[nbands];
199  float gamma, cos2x, pmol;
200  float csenz, transav;
201 
202  float par0, integtrans, delt, saavg;
203  float tmstep, sz, cossz, tg, td;
204  float sola;
205  float asavg, func[11], ra[nbands], rc[nbands];
206 
207  float sumSa, sumrc, sumtdir, sumtdif, denom, tatemp, ravgsat;
208  float tdir, tdif;
209  float q, q4, taucons, q4tau, f, asat, asbar, abar, acldmax;
210  float fcsolz, fcsenz, fcossz;
211  float tg_oz, tg_wv;
212  float trise, tset;
213  float par;
214 
215  int widx490, widx510;
216  float ang500, ta500, slope;
217 
218  int cloudy = FALSE;
219 
220  static float kavg_oz = 5.2e-05;
221  static float kavg_wv = 0.002;
222  // need some memory allocation
223  float *aerphasefunc;
224  float *omega;
225  float *modelAngstrom;
226 
227  l1str *l1rec = l2rec->l1rec;
228 
229  int16_t year, month, mday, doy;
230  double sec;
231  unix2yds(l1rec->scantime, &year, &doy, &sec);
232  unix2ymds(l1rec->scantime, &year, &month, &mday, &sec);
233 
234  aerphasefunc = (float *) allocateMemory(nbands * sizeof(float), "aerphasefunc");
235  omega = (float *) allocateMemory(nbands * sizeof(float), "omega");
236  modelAngstrom = (float *) allocateMemory(nbands * sizeof(float), "modelAngstrom");
237 
238  int asat_is_max = 0;
239 
240  /****************************************************************************
241  * Perform data checks; assign defaults if necessary:
242  ** If data invalid and uncorrectable, PAR set to BAD_FLT missing value.
243  **
244  ***************************************************************************/
245 
246  if (l1rec->solz[ip] < 0.0f || l1rec->solz[ip] > 90.0f) {
247  printf("Error: subroutine CalcPAR. computation failed.\n");
248  printf("SolZen (outside valid range 0-90) = %f\n", l1rec->solz[ip]);
249  par = BAD_FLT;
250  return par;
251  }
252 
253  if (l1rec->senz[ip] < 0.0f || l1rec->senz[ip] > 90.0f) {
254  printf("Error: subroutine CalcPAR. computation failed.\n");
255  printf("ViewZen (outside valid range 0-90) = %f\n", l1rec->senz[ip]);
256  par = BAD_FLT;
257  return par;
258  }
259 
260  if (l1rec->delphi[ip] < -180.0f || l1rec->delphi[ip] > 180.0f) {
261  printf("Error: subroutine CalcPAR. computation failed.\n");
262  printf("delphi (outside valid range -180:180) = %f\n", l1rec->delphi[ip]);
263  par = BAD_FLT;
264  return par;
265  }
266 
267  float dobson = l1rec->oz[ip];
268 
269  csolz = l1rec->csolz[ip];
270  csenz = l1rec->csenz[ip];
271  cos2x = powf(cosf((l1rec->scattang[ip] * M_PI / 180.0)), 2.0);
272 
273  if (l1rec->scattang[ip] < 0.0 || l1rec->scattang[ip] > 180.0) {
274  printf("Error: subroutine CalcPAR. computation failed.\n");
275  printf("scatang (outside valid range 0-180) = %f\n", l1rec->scattang[ip]);
276  printf("scatang = -(cossz*zosvz)+(sinsz*sinvz*cosra)\n");
277  par = BAD_FLT;
278  return par;
279  }
280 
281  if (taua <= 0.0f) {
282  taua = 0.1f;
283  }
284 
285  float surfpress = l1rec->pr[ip];
286  if (surfpress <= 0.0f) {
287  surfpress = EARTH_SURF_PRESSURE;
288  }
289 
290  /* Need to use conventional Angstrom exponent.
291  if it's missing GetAerPhase should work and return a default value */
292 
293  GetAerPhase(l2rec, ip, nbands, angstrom, aerphasefunc, omega, modelAngstrom);
294  if (dobson < 0.0f) {
295  dobson = EstimateDobson(year, month, mday, l1rec->lat[ip]);
296  }
297 
298  float watvap = l1rec->wv[ip];
299  if (watvap < 0.0f) {
300  watvap = EstimateWatVap(year, month, mday, l1rec->lat[ip]);
301  }
302 
303  /*****************************************************************************
304  * Transform TOA radiances into reflectances: *
305  ****************************************************************************/
306  // TODO: consider passing in rhot instead of Lt
307  float Airmass = (1.f / csolz) + (1.f / csenz);
308 
309  for (band = 0; band < nbands; band++) {
310  r[band] = (M_PI * (Lt[band])) / (fo[band] * l1rec->fsol * csolz);
311  // Compute gaseous transmittance:
312  transg = expf((-ko3[band] * dobson * Airmass));
313  rp[band] = r[band] / transg;
314  }
315 
316  /****************************************************************************
317  * Compute reflectances of the cloud-surface subsystem:
318  **
319  ***************************************************************************/
320 
321  /* Setup for transmittance equation:*/
322 
323  ps0 = EARTH_SURF_PRESSURE;
324  del = 0.0095f;
325 
326  pmol = (2.0f * (1.0f - del) * 0.75f * (1.0f + cos2x) + 3.0f * del) / (2.0f + del);
327 
328  asymfac = 0.6666667f;
329  beta = 0.5f * (1.0f + asymfac);
330  gamma = 1.0f - asymfac;
331 
332  /******************************************************************************
333  * Compute diffuse atmospheric transmittance:
334  **
335  ******************************************************************************/
336  sumSa = 0.0f;
337  sumrc = 0.0f;
338  sumtdir = 0.0f;
339  sumtdif = 0.0f;
340  denom = 0.0f;
341 
342  weight[0] = (wl[0] - 400.0) + (wl[1] - wl[0]) / 2.0;
343  for (band = 1; band < nbands - 1; band++) {
344  weight[band] = (wl[band] - wl[band - 1]) / 2.0 + (wl[band + 1] - wl[band]) / 2.0;
345  }
346  weight[nbands - 1] = (700.0 - wl[nbands - 1]) + (wl[nbands - 1] - wl[nbands - 2]) / 2.0;
347 
348  for (band = 0; band < nbands; band++) {
349  taumol[band] = taumolbar[band] * (surfpress / ps0);
350 
351  tauaer[band] = taua * powf(wl[band] / (float) input->aer_wave_long, modelAngstrom[band]);
352 
353  tau[band] = taumol[band] + gamma * tauaer[band];
354 
355  ra[band] = ((taumol[band] * pmol) + (omega[band] * tauaer[band] * aerphasefunc[band])) /
356  (4.0f * csolz * csenz);
357 
358  /* Compute diffuse atmospheric transmittance:*/
359 
360  transa = expf(-(taumol[band] + tauaer[band]) / csolz) *
361  expf((0.52f * taumol[band] + beta * tauaer[band]) / csolz);
362  transav = expf(-(taumol[band] + tauaer[band]) / csenz) *
363  expf((0.52f * taumol[band] + beta * tauaer[band]) / csenz);
364 
365  // Get direct and diffuse transmittances
366  tdir = expf(-tau[band] / csolz);
367 
368  tdif =
369  (expf(-tau[band] / csolz)) * (expf((0.52f * taumol[band] + beta * tauaer[band]) / csolz) - 1.f);
370  sumtdir += (tdir * fo[band]);
371  sumtdif += (tdif * fo[band]);
372 
373  /* Compute spherical albedo of the atmosphere:*/
374  sa = expf(-(taumol[band] + gamma * tauaer[band])) * (0.92f * taumol[band] + gamma * tauaer[band]);
375 
376  sumSa += sa * fo[band] * weight[band];
377  denom += fo[band] * weight[band];
378 
379  rc[band] = (rp[band] - ra[band]) / ((transa * transav) + (sa * (rp[band] - ra[band])));
380  sumrc += (rc[band] * fo[band] * weight[band]);
381  }
382 
383  // Derive tauaer(500) from (469,555)
384  // indexes into modelAngstrom need to use bindex
385  switch (l1rec->l1file->sensorID) {
386  case MODIST:
387  case MODISA:
388  widx490 = windex(469., wl, nbands);
389  widx510 = windex(555., wl, nbands);
390  break;
391  case VIIRSN:
392  case VIIRSJ1:
393  case VIIRSJ2:
394  widx490 = windex(490., wl, nbands);
395  widx510 = windex(555., wl, nbands);
396  break;
397  case OCTS:
398  widx490 = windex(490., wl, nbands);
399  widx510 = windex(520., wl, nbands);
400  break;
401  default:
402  widx490 = windex(490., wl, nbands);
403  widx510 = windex(510., wl, nbands);
404  break;
405  }
406 
407  slope = (modelAngstrom[widx510] - modelAngstrom[widx490]) / (wl[widx510] - wl[widx490]);
408  ang500 = modelAngstrom[widx490] + slope * (wl[widx510] - wl[widx490]);
409  ta500 = taua * powf((500.0 / (float) input->aer_wave_long), ang500);
410 
411  /******************************************************************************
412  * Compute daily average PAR:
413  **
414  ******************************************************************************/
415 
416  /* Set up terms constant in the integration:*/
417  // Sa_bar:
418  saavg = sumSa / denom;
419  // R_bar:
420  ravgsat = sumrc / denom;
421 
422  /******************************************************************************
423  *Get simple bi-directional reflectance factor:
424  **
425  ******************************************************************************/
426 
427  q = 1.0f / (3.0f * (1.0f - 0.853f));
428  q4 = 4.0f * q;
429 
430  taucons = TAUCONS_LOW;
431  q4tau = q4 / (taucons + q4);
432 
433  fcsenz = (3.f / 7.f) * (1.0f + 2.0f * csenz);
434  fcsolz = (3.f / 7.f) * (1.0f + 2.0f * csolz);
435 
436  f = (1.0f - q4tau * fcsolz) /
437  (0.49f * ((1.f + 4.f * csolz * csenz) / (csolz + csenz)) - q4tau * fcsolz * fcsenz);
438 
439  /******************************************************************************
440  * Get As_bar:
441  **
442  ******************************************************************************/
443  asbar = interp_as_taulow(csolz, ta500);
444 
445  /******************************************************************************
446  * Get Albedo from R_bar:
447  **
448  ******************************************************************************/
449  asat = (ravgsat - asbar) * f + asbar;
450 
451  /******************************************************************************
452  * Test Asat to see if it's too large; adjust as necessary:
453  **
454  ******************************************************************************/
455 
456  acldmax = 1.0f - (q4 / (TAUCONS_HIGH + q4)) * fcsolz;
457 
458  if ((asat - asbar) >= acldmax) {
459  asat_is_max = 1;
460  taucons = TAUCONS_HIGH;
461  q4tau = q4 / (taucons + q4);
462  } else {
463  asat_is_max = 0;
464  taucons = q4 * ((fcsolz / (1.0f - (asat - asbar))) - 1.0f);
465 
466  /***************************************************************************
467  * Compute new Asat if needed:
468  **
469  ***************************************************************************/
470 
471  if (taucons > TAUCONS_LOW) {
472  q4tau = q4 / (taucons + q4);
473 
474  f = (1.0f - q4tau * fcsolz) /
475  (0.49f * ((1.f + 4.f * csolz * csenz) / (csolz + csenz)) - q4tau * fcsolz * fcsenz);
476  asat = (ravgsat - asbar) * f + asbar;
477 
478  /* New Test:*/
479 
480  if ((asat - asbar) >= acldmax) {
481  asat_is_max = 1;
482  taucons = TAUCONS_HIGH;
483  q4tau = q4 / (taucons + q4);
484  }
485  } else {
486  taucons = TAUCONS_LOW;
487  q4tau = q4 / (taucons + q4);
488  }
489  }
490 
491  // Test to estimate if pixel is clear or cloudy:
492  cloudy = FALSE;
493  if (asat > asbar)
494  cloudy = TRUE;
495 
496  /*****************************************************************************
497  * Get parameters used in integration time-steps:
498  **
499  *****************************************************************************/
500 
501  /* Get the rise and set time for this day (for integral):*/
502 
503  triseset(doy, l1rec->lon[ip], l1rec->lat[ip], &trise, &tset);
504 
505  /* Set up for 10 intervals*/
506  delt = (tset - trise) / 10.0f;
507 
508  /***************************************************************************
509  * Get parameters used in integration time-steps:
510  **
511  ***************************************************************************/
512 
513  for (i = 0; i < 11; i++) {
514  tmstep = trise + i * delt;
515  int intyear, intday;
516  intyear = year;
517  intday = doy;
518  sunangs_(&intyear, &intday, &tmstep, l1rec->lon + ip, l1rec->lat + ip, &sz, &sola);
519  cossz = cosf(sz * (float) M_PI / 180.0f);
520  fcossz = (3.f / 7.f) * (1.0f + 2.0f * cossz);
521 
522  /***************************************************************************
523  * Get gaseous transmittance term:
524  **
525  ***************************************************************************/
526  if (sz >= MAX_SOLZEN) {
527  tg_oz = 0.0f;
528  tg_wv = 0.0f;
529  } else {
530  // kavg_oz (5.2E-5) is valid for dobson in true dobson units, not
531  // the ozone units in l2gen
532  tg_oz = expf(-kavg_oz * powf((dobson * 1000. / cossz), 0.99f));
533  tg_wv = expf(-kavg_wv * powf((watvap / cossz), 0.87f));
534  }
535 
536  tg = tg_oz * tg_wv;
537  /***************************************************************************
538  * Get average diffuse atmospheric transmittance and As_bar terms: *
539  ***************************************************************************/
540  td = 0.0f;
541 
542  if (sz < MAX_SOLZEN) {
543  for (band = 0; band < nbands; band++) {
544  tatemp = expf(-(taumol[band] + tauaer[band]) / cossz) *
545  expf((0.52f * taumol[band] + beta * tauaer[band]) / cossz);
546  td += tatemp * fo[band] * weight[band];
547  }
548  td /= denom;
549  }
550 
551  if (cossz < 0.05) {
552  cossz = 0.05;
553  }
554  if (cloudy) {
555  asavg = interp_as_tauhigh(cossz);
556  } else {
557  asavg = interp_as_taulow(cossz, ta500);
558  }
559 
560  /*****************************************************************************
561  *Get A_bar term:
562  **
563  *****************************************************************************/
564 
565  if (asat_is_max) {
566  abar = 1.0f - q4tau * fcossz;
567  } else {
568  abar = asat * (1.0f - q4tau * fcossz) / (1.0f - q4tau * fcsolz);
569  }
570 
571  /*new test:*/
572  if (abar >= 1.0f) {
573  taucons = TAUCONS_HIGH;
574  q4tau = q4 / (taucons + q4);
575  abar = 1.0f - q4tau * fcossz;
576  }
577 
578  /****************************************************************************
579  * Get "function" integrand term:
580  **
581  ****************************************************************************/
582 
583  if (abar <= asavg) {
584  func[i] = cossz * tg * td / (1.0 - abar * saavg);
585  } else {
586  func[i] = cossz * tg * td * (1.0f - abar) / (1.0f - asavg) / (1.0f - abar * saavg);
587  }
588  }
589 
590  /* Use trapezoidal rule for integration:*/
591 
592  integtrans = 0.0f;
593 
594  for (i = 0; i < 10; i++) {
595  integtrans += delt * (func[i] + func[i + 1]) / 2.0f;
596  }
597 
598  /* Finally, determine the PAR:*/
599  // par0 determined as the mean value [400-700nm] for Thuillier 2003
600  // from run/data/common/Thuillier_F0.dat
601  par0 = PAR_0;
602  par = par0 * l1rec->fsol * integtrans / 24.0f;
603 
604  free(aerphasefunc);
605  free(omega);
606  free(modelAngstrom);
607 
608  return par;
609 }
610 
611 extern float *par_planar_a_inst, *par_planar_b_inst;
612 
613 float calc_par_impl_of_2023(l2str *l2rec, int ip, int nbands, float *Lt, float taua, float angstrom,
614  float *wl, float *fo, float *ko3, float *taumolbar, float *parb, float *parc) {
615  if (ini == 0) {
616  get_luts_data(l2rec, &luts_data);
617  ini = 1;
618  unix2yds(l2rec->l1rec->scantime, &year, &doy, &sec);
619  unix2ymds(l2rec->l1rec->scantime, &year, &month, &mday, &sec);
620  grid_ozone[0] = luts_data.ozonedims.days;
621  grid_ozone[1] = luts_data.ozonedims.latitude;
622  grid_watvap[0] = luts_data.watvapdims.days;
623  grid_watvap[1] = luts_data.watvapdims.latitude;
624 
625  grid_tg[0] = luts_data.tgdims.wavelength;
626  grid_tg[1] = luts_data.tgdims.air_mass;
627  grid_tg[2] = luts_data.tgdims.water_vapor_pressure;
628  grid_tg[3] = luts_data.tgdims.ozone_concentration;
629 
630  grid_td[0] = luts_data.tddims.wavelength;
631  grid_td[1] = luts_data.tddims.air_mass;
632  grid_td[2] = luts_data.tddims.angstrom_coefficients;
633  grid_td[3] = luts_data.tddims.optical_depth_at_550_nm;
634 
635  grid_rho[0] = luts_data.rhodims.solar_zenith_angle;
636  grid_rho[1] = luts_data.rhodims.view_zenith_angle;
637  grid_rho[2] = luts_data.rhodims.relative_azimuth_angle;
638  grid_rho[3] = luts_data.rhodims.angstrom_coefficients;
639  grid_rho[4] = luts_data.rhodims.optical_depth_at_550_nm;
640  grid_rho[5] = luts_data.rhodims.wavelength;
641 
642  grid_scalar_par[0] = luts_data.scalar_luts.wind_speed;
643  grid_scalar_par[1] = luts_data.scalar_luts.doy;
644  grid_scalar_par[2] = luts_data.scalar_luts.latitude;
645 
646  grid_scalar_inst_par[0] = luts_data.scalar_inst_luts.wind_speed;
647  grid_scalar_inst_par[1] = luts_data.scalar_inst_luts.cos_solz;
648  grid_scalar_inst_par[2] = luts_data.scalar_inst_luts.cot;
649  observed_time = sec / 3600;
650  {
651  size_t ntimes = l2rec->l1rec->cld_rad->ntimes;
652  t_step = ntimes;
653  t_range = l2rec->l1rec->cld_rad->timecldrange;
654  }
655 
656  weight_5nm[0] = (wl_5nm[0] - 400.0) + (wl_5nm[1] - wl_5nm[0]) / 2.0;
657  for (size_t j = 1; j < N_5nm - 1; j++)
658  weight_5nm[j] = (wl_5nm[j] - wl_5nm[j - 1]) / 2.0 + (wl_5nm[j + 1] - wl_5nm[j]) / 2.0;
659 
660  weight_5nm[N_5nm - 1] = (700.0 - wl_5nm[N_5nm - 1]) + (wl_5nm[N_5nm - 1] - wl_5nm[N_5nm - 2]) / 2.0;
661  for (size_t band = 0; band < N_5nm; band++)
662  Denom_5nm += F0_5nm[band] * weight_5nm[band];
663  // printf("Denom_5nm %f\n", Denom_5nm);
664  }
665  // calcs go here
666  l1str *l1rec = l2rec->l1rec;
667  float par = BAD_FLT;
668  // int iscan = l2rec->l1rec->iscan;
669  if (l1rec->solz[ip] < 0.0f || l1rec->solz[ip] > 90.0f) {
670  printf("Error: subroutine CalcPAR. computation failed.\n");
671  printf("SolZen (outside valid range 0-90) = %f\n", l1rec->solz[ip]);
672  return par;
673  }
674 
675  if (l1rec->senz[ip] < 0.0f || l1rec->senz[ip] > 90.0f) {
676  printf("Error: subroutine CalcPAR. computation failed.\n");
677  printf("ViewZen (outside valid range 0-90) = %f\n", l1rec->senz[ip]);
678  return par;
679  }
680 
681  if (l1rec->delphi[ip] < -180.0f || l1rec->delphi[ip] > 180.0f) {
682  printf("Error: subroutine CalcPAR. computation failed.\n");
683  printf("delphi (outside valid range -180:180) = %f\n", l1rec->delphi[ip]);
684  return par;
685  }
686 
687  float csolz = l1rec->csolz[ip];
688 
689  if (l1rec->scattang[ip] < 0.0 || l1rec->scattang[ip] > 180.0) {
690  printf("Error: subroutine CalcPAR. computation failed.\n");
691  printf("scatang (outside valid range 0-180) = %f\n", l1rec->scattang[ip]);
692  printf("scatang = -(cossz*zosvz)+(sinsz*sinvz*cosra)\n");
693  return par;
694  }
695  float weight[nbands];
696  weight[0] = (wl[0] - 400.0) + (wl[1] - wl[0]) / 2.0;
697  for (size_t band = 1; band < nbands - 1; band++) {
698  weight[band] = (wl[band] - wl[band - 1]) / 2.0 + (wl[band + 1] - wl[band]) / 2.0;
699  }
700  weight[nbands - 1] = (700.0 - wl[nbands - 1]) + (wl[nbands - 1] - wl[nbands - 2]) / 2.0;
701  // replace optical depth, surface pressure and angstrom coefficients with
702  // their default values
703 
704  if (taua <= 0.0f) {
705  taua = TAU_550_LOW_BOUND;
706  angstrom = ANGSTROM_DEFAULT_VALUE;
707  }
708  if (angstrom < -1.0f)
709  angstrom = ANGSTROM_DEFAULT_VALUE;
710  float surfpress = l1rec->pr[ip];
711  if (surfpress <= 0.0f) {
712  surfpress = EARTH_SURF_PRESSURE;
713  }
714 
715  float dobson = l1rec->oz[ip] * 1000;
716  float watvap = l1rec->wv[ip];
717  if (dobson < 0) {
718  float point[] = {doy, l1rec->lat[ip]};
719  size_t dims[] = {luts_data.ozonedims.dimdays, luts_data.ozonedims.dimlatitude};
720  dobson = interpnd(2, dims, point, grid_ozone, luts_data.lut_dobson);
721  }
722  if (watvap < 0) {
723  {
724  float point[] = {doy, l1rec->lat[ip]};
725  size_t dims[] = {luts_data.watvapdims.dimdays, luts_data.watvapdims.dimlatitude};
726  watvap = interpnd(2, dims, point, grid_watvap, luts_data.lut_watvap);
727  }
728  }
729 
730  float rho[nbands];
731  float airmass = kasten_equation(l1rec->solz[ip]);
732  float airmass_0 = kasten_equation(l1rec->senz[ip]);
733  {
734  size_t dims[] = {luts_data.tgdims.dimwavelength, luts_data.tgdims.dimair_mass,
735  luts_data.tgdims.dimwater_vapor_pressure, luts_data.tgdims.dimozone_concentration};
736  for (size_t band = 0; band < nbands; band++) {
737  rho[band] = (M_PI * (Lt[band])) / (fo[band] * l1rec->fsol * csolz);
738  if (rho[band] < 0.0 || rho[band] > 1.0) {
739  parb[ip] = BAD_FLT;
740  parc[ip] = BAD_FLT;
741  return BAD_FLT;
742  }
743  float point[] = {wl[band], airmass, watvap, dobson};
744  float tg_solzen = interp4d(dims, point, grid_tg,
745  luts_data.lut_tg); // the same so it seems
746  point[1] = airmass_0;
747  float tg_viewzen = interp4d(dims, point, grid_tg, luts_data.lut_tg);
748  rho[band] = rho[band] / tg_viewzen / tg_solzen;
749  }
750  }
751  float Ra[nbands], TauMol[nbands], TauAer[nbands], TransA[nbands], TransAv[nbands], Sa[nbands];
752  size_t dims[] = {luts_data.tddims.dimwavelength, luts_data.tddims.dimair_mass,
753  luts_data.tddims.dimangstrom_coefficients, luts_data.tddims.dimoptical_depth_at_550_nm};
754  size_t dims_rho[] = {
755  luts_data.rhodims.dimsolar_zenith_angle, luts_data.rhodims.dimview_zenith_angle,
756  luts_data.rhodims.dimrelative_azimuth_angle, luts_data.rhodims.dimangstrom_coefficients,
757  luts_data.rhodims.dimoptical_depth_at_550_nm, luts_data.rhodims.dimwavelength};
758 
759  float relaz = l1rec->sola[ip] - l1rec->sena[ip];
760  if (relaz < 0)
761  relaz += 360.0;
762  if (relaz > 180.0)
763  relaz = 360.0 - relaz;
764  for (size_t band = 0; band < nbands; band++) {
765  TauMol[band] = ((A / pow(wl[band] / 1000.0, 4)) + (B / pow(wl[band] / 1000.0, 5)) +
766  (C / pow(wl[band] / 1000.0, 6))) *
767  (surfpress / Ps0);
768  TauAer[band] = taua * pow((550.0 / wl[band]), angstrom);
769 
770  float point[] = {wl[band], airmass, angstrom, taua};
771  TransA[band] = interp4d(dims, point, grid_td, luts_data.lut_td);
772 
773 
774  point[1] = airmass_0;
775  TransAv[band] = interp4d(dims, point, grid_td, luts_data.lut_td);
776  // Compute spherical albedo of the atmosphere:
777  Sa[band] =
778  exp(-(TauMol[band] + Gamma0 * TauAer[band])) * (0.92 * TauMol[band] + Gamma0 * TauAer[band]);
779 
780  float point_rho[] = {l1rec->solz[ip], l1rec->senz[ip], relaz, angstrom, taua, wl[band]};
781  Ra[band] = interp6d(dims_rho, point_rho, grid_rho, luts_data.lut_rho);
782  }
783  float Sa_bar = 0.0f;
784  float Denom = 0.0f;
785  for (size_t band = 0; band < nbands; band++) {
786  Sa_bar += Sa[band] * fo[band] * weight[band];
787  Denom += fo[band] * weight[band];
788  }
789  Sa_bar /= Denom;
790 
791  float CF_obs, TauCld_obs;
792  float albe_obs[nbands];
793 
794  getcldalbe(l1rec->cld_rad->taucld[ip], l1rec->cld_rad->cfcld[ip], csolz, observed_time, t_range, albe_obs,
795  &TauCld_obs, &CF_obs, t_step, wl, nbands);
796  // averaged observed albedo
797  float AlbeCld_obs_bar = 0.0;
798  for (size_t band = 0; band < nbands; band++)
799  AlbeCld_obs_bar += albe_obs[band] / (float) nbands;
800  // Get simple bi-directional reflectance factor, needs a refer
801 
802  const float q = 1.0 / (3.0 * (1.0 - 0.853));
803  const float q4 = 4.0 * q;
804  float TauCons[nbands];
805  float q4tau[nbands];
806  float F[nbands];
807  const float cosvz = l1rec->csenz[ip];
808  const float fCosVZ = (3. / 7.) * (1.0 + 2.0 * cosvz);
809  const float fCosSolZen = (3. / 7.) * (1.0 + 2.0 * csolz);
810 
811  for (size_t band = 0; band < nbands; band++) {
812  if (TauCld_obs >= 5.0) {
813  TauCons[band] = TauCld_obs;
814  } else {
815  TauCons[band] = 5.0;
816  }
817 
818  q4tau[band] = q4 / (TauCons[band] + q4);
819  F[band] = (1.0 - q4tau[band] * fCosSolZen) /
820  (0.49 * ((1. + 4. * csolz * cosvz) / (csolz + cosvz)) - q4tau[band] * fCosSolZen * fCosVZ);
821  }
822  // # new parameterization of As
823 
824  float fr;
825  float As[nbands], Asat[nbands];
826  float As_bar = 0.0f, Asat_bar = 0.0f;
827  const float chl = 0.1;
828 
829  const float airmass2 = kasten_equation(l1rec->solz[ip]) * kasten_equation(l1rec->senz[ip]);
830  const float wind_speed = l1rec->ws[ip];
831  const float AsatMax = 1.0 - (q4 / (300.0 + q4)) * fCosSolZen;
832  // solving quadratic equation
833  float a, b, c;
834 
835  for (size_t band = 0; band < nbands; band++) {
836  fr = (1 - exp(-(TauMol[band] + TauAer[band]) / csolz) / TransA[band]) * (1. - CF_obs) + CF_obs;
837  As[band] = getosa(wl[band], l1rec->solz[ip], wind_speed, chl, fr, &luts_data);
838  a = -(1 / F[band] - 1) * exp(-(TauMol[band] + TauAer[band]) * airmass2) * Sa[band];
839  b = (1 / F[band] - 1) * exp(-(TauMol[band] + TauAer[band]) * airmass2) * Sa[band] * As[band] +
840  (1 / F[band] - 1) * exp(-(TauMol[band] + TauAer[band]) * airmass2) +
841  TransA[band] * TransAv[band] + Sa[band] * (rho[band] - Ra[band]);
842  c = -(1 / F[band] - 1) * exp(-(TauMol[band] + TauAer[band]) * airmass2) * As[band] -
843  (rho[band] - Ra[band]);
844  const float discr = b * b / 4. / a / a - c / a;
845  Asat[band] = As[band];
846  if (discr > 0.0 && isfinite(discr) && a != 0) {
847  const float d = -sqrt(discr) - b / 2. / a;
848  if ((d > 0) && (d < 1)) {
849  Asat[band] = d;
850  } else {
851  Asat[band] = sqrt(discr) - b / 2. / a;
852  }
853  }
854  if (Asat[band] < As[band])
855  Asat[band] = As[band];
856  As_bar += As[band] * fo[band] * weight[band] / Denom;
857  Asat_bar += Asat[band] * fo[band] * weight[band] / Denom;
858  }
859  // instantaneous PAR
860  float Func_inst = 0.0f, Funcb_inst = 0.0f;
861  {
862  float Tg_bar_inst = 0.0;
863  float Td_bar_inst = 0.0;
864  float Td_inst[nbands];
865  float Tg_inst[nbands];
866  float CosSZ = l1rec->csolz[ip];
867  float fCosSZ = (3. / 7.) * (1.0 + 2.0 * CosSZ);
868  for (size_t band = 0; band < nbands; band++) {
869  airmass = kasten_equation(l1rec->solz[ip]);
870  float point_tg[] = {wl[band], airmass, watvap, dobson};
871  float point_td[] = {wl[band], airmass, angstrom, taua};
872  size_t dims_tg[] = {luts_data.tgdims.dimwavelength, luts_data.tgdims.dimair_mass,
873  luts_data.tgdims.dimwater_vapor_pressure,
874  luts_data.tgdims.dimozone_concentration};
875  size_t dims_td[] = {luts_data.tddims.dimwavelength, luts_data.tddims.dimair_mass,
876  luts_data.tddims.dimangstrom_coefficients,
877  luts_data.tddims.dimoptical_depth_at_550_nm};
878  Tg_inst[band] = interp4d(dims_tg, point_tg, grid_tg, luts_data.lut_tg);
879  Td_inst[band] = interp4d(dims_td, point_td, grid_td, luts_data.lut_td);
880  Tg_bar_inst += Tg_inst[band] * fo[band] * weight[band] / Denom;
881  Td_bar_inst += Td_inst[band] * fo[band] * weight[band] / Denom;
882  }
883  float Abar = Asat_bar >= AsatMax ? AsatMax : Asat_bar;
884  if (Abar >= 1.0)
885  Abar = 1.0 - q4 / (300. + q4) * fCosSZ;
886  if (Abar < 0)
887  Abar = As_bar;
888  if (Abar <= As_bar)
889  Func_inst = CosSZ * Tg_bar_inst * Td_bar_inst / (1.0 - As_bar * Sa_bar);
890  else
891  Func_inst = CosSZ * Tg_bar_inst * Td_bar_inst * (1. - Abar) / (1.0 - As_bar) / (1.0 - Abar * Sa_bar);
892  Funcb_inst= CosSZ * Tg_bar_inst * Td_bar_inst * (1. - Abar) / (1.0 - Abar * Sa_bar); // # below surface
893 
894  const float fo_bar = PAR_0_2023; // 176.41f;
895  par_planar_a_inst[ip] = Func_inst * l1rec->fsol * fo_bar / 3600 / 24.0;
896  par_planar_b_inst[ip] = Funcb_inst * l1rec->fsol * fo_bar / 3600 / 24.0;
897 
898 
899  }
900 
901 
902  // time integration starts here
903  float Func_spectral[t_step][nbands], Funcb_spectral[t_step][nbands], Funcc_spectral[t_step][nbands];
904  float Func[t_step], Funcb[t_step], Funcc[t_step];
905  float trise, tset;
906  triseset(doy, l1rec->lon[ip], l1rec->lat[ip], &trise, &tset);
907  if (observed_time > trise && observed_time > tset) {
908  trise += 24;
909  tset += 24;
910  }
911  if (observed_time < trise && observed_time < tset) {
912  tset -= 24;
913  trise -= 24;
914  }
915  size_t st_int = t_step;
916  size_t end_int = 0;
917  for (size_t it = 0; it < t_step; it++) { // need to check
918  float sz = get_solz(doy, t_range[it], l1rec->lon[ip], l1rec->lat[ip]);
919  if (sz >= 90.0 || t_range[it] < trise || t_range[it] > tset) {
920  Func[it] = 0.0f;
921  Funcb[it] = 0.0f;
922  Funcc[it] = 0.0f;
923  for (size_t band = 0; band < nbands; band++) {
924  Func_spectral[it][band] = 0;
925  Funcb_spectral[it][band] = 0;
926  Funcc_spectral[it][band] = 0;
927  }
928 
929  continue;
930  }
931  if (st_int > it)
932  st_int = it;
933  if (end_int < it)
934  end_int = it;
935  float As2_bar = 0.0f, Tg_bar = 0.0f, Td_bar = 0.0f, AlbeCld_bar = 0.0f;
936  float Abar_spectral[nbands], Abar = 0.0f;
937  float Aavg_spectral[nbands], Aavg = 0.0f;
938  float Tg[nbands];
939  float Td[nbands];
940  // float temp1[nbands], temp2[nbands], temp3[nbands];
941  float As2[nbands];
942  float CosSZ = cos(sz * M_PI / 180.);
943  float albe_cld[nbands], cf, tau;
944  getcldalbe(l1rec->cld_rad->taucld[ip], l1rec->cld_rad->cfcld[ip], CosSZ, t_range[it], t_range,
945  albe_cld, &tau, &cf, t_step, wl,
946  nbands); // something is wrong with getcldalbedo
947  if (CosSZ < 0.05)
948  CosSZ = 0.05;
949  float fCosSZ = (3. / 7.) * (1.0 + 2.0 * CosSZ);
950  size_t dims_tg[] = {luts_data.tgdims.dimwavelength, luts_data.tgdims.dimair_mass,
951  luts_data.tgdims.dimwater_vapor_pressure,
952  luts_data.tgdims.dimozone_concentration};
953  size_t dims_td[] = {luts_data.tddims.dimwavelength, luts_data.tddims.dimair_mass,
954  luts_data.tddims.dimangstrom_coefficients,
955  luts_data.tddims.dimoptical_depth_at_550_nm};
956  airmass = kasten_equation(sz);
957  for (size_t band = 0; band < nbands; band++) {
958  float point_tg[] = {wl[band], airmass, watvap, dobson};
959  float point_td[] = {wl[band], airmass, angstrom, taua};
960  Tg[band] = interp4d(dims_tg, point_tg, grid_tg, luts_data.lut_tg);
961  Td[band] = interp4d(dims_td, point_td, grid_td, luts_data.lut_td);
962  fr = (1. - exp(-(TauMol[band] + TauAer[band]) / CosSZ) / Td[band]) * (1. - cf) + cf;
963  As2[band] = getosa(wl[band], sz, wind_speed, chl, fr, &luts_data);
964  As2_bar += As2[band] * fo[band] * weight[band] / Denom;
965  Tg_bar += Tg[band] * fo[band] * weight[band] / Denom;
966  Td_bar += Td[band] * fo[band] * weight[band] / Denom;
967  AlbeCld_bar += albe_cld[band] / (float) nbands;
968  }
969  // for (size_t band = 0; band < N_5nm; band++) {
970  // float point_tg[] = {wl_5nm[band], airmass, watvap, dobson};
971  // float point_td[] = {wl_5nm[band], airmass, angstrom, taua};
972  // Tg_bar +=
973  // interp4d(dims_tg, point_tg, grid_tg, luts_data.lut_tg) * F0_5nm[band] * weight_5nm[band];
974  // Td_bar +=
975  // interp4d(dims_td, point_td, grid_td, luts_data.lut_td) * F0_5nm[band] * weight_5nm[band];
976  // }
977  // Tg_bar /= Denom_5nm;
978  // Td_bar /= Denom_5nm;
979  for (size_t band = 0; band < nbands; band++) {
980  Aavg_spectral[band] = Asat[band] >= AsatMax ? AsatMax : Asat[band];
981  Abar_spectral[band] = Aavg_spectral[band] +
982  (albe_cld[band] + As2[band] - 2 * albe_cld[band] * As2[band]) -
983  (albe_obs[band] + As[band] - 2 * albe_obs[band] * As[band]);
984  if (Abar_spectral[band] >= 1.0)
985  Abar_spectral[band] = 1.0 - q4 / (300. + q4) * fCosSZ;
986  if (Abar_spectral[band] < 0)
987  Abar_spectral[band] = As2[band];
988 
989  if (Abar_spectral[band] <= As2[band])
990  Func_spectral[it][band] = CosSZ * Tg[band] * Td[band] / (1.0 - As2[band] * Sa[band]);
991  else
992  Func_spectral[it][band] = CosSZ * Tg[band] * Td[band] * (1. - Abar_spectral[band]) /
993  (1.0 - As2[band]) / (1.0 - Abar_spectral[band] * Sa[band]);
994  Funcb_spectral[it][band] = CosSZ * Tg[band] * Td[band] * (1.0 - Abar_spectral[band]) /
995  (1.0 - Abar_spectral[band] * Sa[band]); // # below surface
996 
997  Funcc_spectral[it][band] = CosSZ * Tg[band] * Td[band] / (1.0 - As2[band] * Sa[band]);
998  }
999  Aavg = Asat_bar >= AsatMax ? AsatMax : Asat_bar;
1000  Abar = Aavg + (AlbeCld_bar + As2_bar - 2 * AlbeCld_bar * As2_bar) -
1001  (AlbeCld_obs_bar + As_bar - 2 * AlbeCld_obs_bar * As_bar);
1002  if (Abar >= 1.0)
1003  Abar = 1.0 - q4 / (300. + q4) * fCosSZ;
1004  if (Abar < 0)
1005  Abar = As2_bar;
1006  if (Abar <= As2_bar)
1007  Func[it] = CosSZ * Tg_bar * Td_bar / (1.0 - As2_bar * Sa_bar);
1008  else
1009  Func[it] = CosSZ * Tg_bar * Td_bar * (1. - Abar) / (1.0 - As2_bar) / (1.0 - Abar * Sa_bar);
1010  Funcb[it] = CosSZ * Tg_bar * Td_bar * (1. - Abar) / (1.0 - Abar * Sa_bar); // # below surface
1011  Funcc[it] = CosSZ * Tg_bar * Td_bar / (1.0 - As2_bar * Sa_bar);
1012  }
1013 
1014  float PAR_spectral[nbands], PARb_spectral[nbands], PARc_spectral[nbands];
1015  float IntegTrans, IntegTransb, IntegTransc;
1016  float Par = 0, Parb = 0, Parc = 0;
1017  const float min_shift = t_range[st_int] - trise;
1018  const float max_shift = tset - t_range[end_int];
1019  size_t time_point_start = st_int;
1020  size_t time_point_end = end_int;
1021 
1022  // end debug print
1023  for (size_t band = 0; band < nbands; band++) {
1024  IntegTrans = 0;
1025  IntegTransb = 0;
1026  IntegTransc = 0;
1027  // trapezoid integration
1028  for (size_t it = st_int; it < end_int; it++) {
1029  IntegTrans += (Func_spectral[it][band] + Func_spectral[it + 1][band]) / 2.0;
1030  IntegTransb += (Funcb_spectral[it][band] + Funcb_spectral[it + 1][band]) / 2.0;
1031  IntegTransc += (Funcc_spectral[it][band] + Funcc_spectral[it + 1][band]) / 2.0;
1032  }
1033  IntegTrans += min_shift * Func_spectral[time_point_start][band] / 2.0 +
1034  max_shift * Func_spectral[time_point_end][band] / 2.0;
1035  IntegTransb += min_shift * Funcb_spectral[time_point_start][band] / 2.0 +
1036  max_shift * Funcb_spectral[time_point_end][band] / 2.0;
1037  IntegTransc += min_shift * Funcc_spectral[time_point_start][band] / 2.0 +
1038  max_shift * Funcc_spectral[time_point_end][band] / 2.0;
1039 
1040  PAR_spectral[band] = IntegTrans / 24.0 * fo[band] * l1rec->fsol;
1041  PARb_spectral[band] = IntegTransb / 24.0 * fo[band] * l1rec->fsol;
1042  PARc_spectral[band] = IntegTransc / 24.0 * fo[band] * l1rec->fsol;
1043  Par += PAR_spectral[band] * weight[band] / 300.0; // needs to modify to fit integration
1044  // step 300/5
1045  Parb += PARb_spectral[band] * weight[band] / 300.0;
1046  Parc += PARc_spectral[band] * weight[band] / 300.0;
1047  }
1048  IntegTrans = 0;
1049  IntegTransb = 0;
1050  IntegTransc = 0;
1051  for (size_t it = st_int; it < end_int; it++) {
1052  IntegTrans += (Func[it] + Func[it + 1]) / 2.0;
1053  IntegTransb += (Funcb[it] + Funcb[it + 1]) / 2.0;
1054  IntegTransc += (Funcc[it] + Funcc[it + 1]) / 2.0;
1055  }
1056 
1057  IntegTrans += min_shift * Func[time_point_start] / 2.0 + max_shift * Func[time_point_end] / 2.0;
1058  IntegTransb += min_shift * Funcb[time_point_start] / 2.0 + max_shift * Funcb[time_point_end] / 2.0;
1059  IntegTransc += min_shift * Funcc[time_point_start] / 2.0 + max_shift * Funcc[time_point_end] / 2.0;
1060  const float fo_bar = PAR_0_2023; // 176.41f;
1061  Par = IntegTrans / 24.0 * fo_bar * l1rec->fsol;
1062  Parb = IntegTransb / 24.0 * fo_bar * l1rec->fsol;
1063  Parc = IntegTransc / 24.0 * fo_bar * l1rec->fsol;
1064  if (l1rec->glint_coef[ip] > MAXGLINT) {
1065  Par = BAD_FLT;
1066  Parb = BAD_FLT;
1067  Parc = BAD_FLT;
1068  }
1069  parb[ip] = Parb;
1070  parc[ip] = Parc;
1071  return Par;
1072 }
MAX_SOLZEN
#define MAX_SOLZEN
Definition: par_utils.h:23
q
data_t q
Definition: decode_rs.h:74
r
int r
Definition: decode_rs.h:73
MAXGLINT
#define MAXGLINT
Definition: par_utils.h:20
j
int j
Definition: decode_rs.h:73
denom
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:93
delta
const float delta
Definition: calc_par.c:104
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
FALSE
#define FALSE
Definition: rice.h:164
calc_par
float calc_par(l2str *l2rec, int ip, int nbands, float *Lt, float taua, float angstrom, float *wl, float *fo, float *ko3, float *taumolbar)
Definition: calc_par.c:192
interp3d
float interp3d(const size_t *dims, const float *point, float **grid, const float *lut)
3-dimensional interpolation
Definition: par_utils.c:1795
calc_scalar_inst_par
void calc_scalar_inst_par(l2str *l2rec, int ip, float par_above_ins, float *par_scalar_ins)
Definition: calc_par.c:123
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
interp6d
float interp6d(const size_t *dims, const float *point, float **grid, const float *lut)
6-dimensional interpolation
Definition: par_utils.c:1675
VIIRSN
#define VIIRSN
Definition: sensorDefs.h:23
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
kasten_equation
float kasten_equation(float solz)
Kasten eq to compute airmass.
Definition: par_utils.c:58
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
TRUE
#define TRUE
Definition: rice.h:165
MODIST
#define MODIST
Definition: sensorDefs.h:18
grid_ozone
float * grid_ozone[2]
Definition: calc_par.c:90
Gamma0
const float Gamma0
Definition: calc_par.c:105
GetAerPhase
void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom, float *aerphasefunc, float *omega, float *modelAngstrom)
Definition: par_utils.c:680
input
instr * input
Definition: main_l2binmatch.cpp:27
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
F
const double F
l2_struc.h
l12_proto.h
grid_td
float * grid_td[4]
Definition: calc_par.c:92
A
const float A
Definition: calc_par.c:100
grid_rho
float * grid_rho[6]
Definition: calc_par.c:92
calc_par_impl_of_2023
float calc_par_impl_of_2023(l2str *l2rec, int ip, int nbands, float *Lt, float taua, float angstrom, float *wl, float *fo, float *ko3, float *taumolbar, float *parb, float *parc)
Definition: calc_par.c:613
grid_scalar_par
float * grid_scalar_par[3]
Definition: calc_par.c:93
par_planar_a_inst
float * par_planar_a_inst
Definition: get_par.c:38
TAUCONS_HIGH
#define TAUCONS_HIGH
Definition: par_utils.h:27
sunangs_
void sunangs_(int *, int *, float *, float *, float *, float *, float *)
Definition: sunangs.c:25
f
double precision function f(R1)
Definition: tmd.lp.f:1454
wl_5nm
const float wl_5nm[]
Definition: calc_par.c:114
getosa
float getosa(float wl, float sza, float wind, float chl, float fr, const luts_par *luts_data)
Definition: par_utils.c:1544
get_luts_data
void get_luts_data(l2str *l2rec, luts_par *luts_data)
Get the luts data object For now it includes the MERRA data. Should be moved after the merra data is ...
Definition: par_utils.c:244
interp_as_tauhigh
float interp_as_tauhigh(float csz)
Definition: par_utils.c:1349
omega
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
getcldalbe
void getcldalbe(float *TauCld, float *CF, float cosSZ, float t_obs, float *t_range, float *albe_obs, float *TauCld_obs, float *CF_obs, size_t t_step, float *wl, size_t bands_num)
//
Definition: par_utils.c:1384
func
subroutine func(x, conec, n, bconecno, bn, units, u, inno, i, outno, o, Input, Targ, p, sqerr)
Definition: ffnet.f:287
slope
float32 slope[]
Definition: l2lists.h:30
rc
float rc(float x, float y)
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
grid_scalar_inst_par
float * grid_scalar_inst_par[3]
Definition: calc_par.c:94
luts_data
luts_par luts_data
Definition: calc_par.c:89
PAR_0
#define PAR_0
Definition: par_utils.h:24
N_5nm
const size_t N_5nm
Definition: calc_par.c:106
triseset
void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset)
Definition: par_utils.c:1131
color_dtdb.gamma
float gamma
Definition: color_dtdb.py:447
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
EstimateWatVap
float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:1002
hico.proc_hico.sola
sola
Definition: proc_hico.py:242
B
const float B
Definition: calc_par.c:101
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
OCTS
#define OCTS
Definition: sensorDefs.h:14
EARTH_SURF_PRESSURE
#define EARTH_SURF_PRESSURE
Definition: par_utils.h:17
interp4d
float interp4d(const size_t *dims, const float *point, float **grid, const float *lut)
4-dimensional interpolation
Definition: par_utils.c:1749
get_solz
float get_solz(int32_t jday, float time, float lon, float lat)
Definition: par_utils.c:1172
beta
double beta
Definition: GEO_earth_location.c:14
geo_eval.cf
cf
Definition: geo_eval.py:29
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
F0_5nm
const float F0_5nm[]
Definition: calc_par.c:107
C
const float C
Definition: calc_par.c:102
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
Ps0
float Ps0
Definition: calc_par.c:103
splint
subroutine splint(xa, ya, y2a, n, x, y)
Definition: atrem_app_refl_f90_cubeio.f:3270
ANGSTROM_DEFAULT_VALUE
#define ANGSTROM_DEFAULT_VALUE
Definition: par_utils.h:19
EstimateDobson
float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:876
grid_watvap
float * grid_watvap[2]
Definition: calc_par.c:91
VIIRSJ2
#define VIIRSJ2
Definition: sensorDefs.h:44
c
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
interpnd
float interpnd(size_t n_dims, const size_t *dims, const float *point, float **grid, const float *lut)
N-dimensional interpolation.
Definition: par_utils.c:158
PAR_0_2023
#define PAR_0_2023
Definition: par_utils.h:25
i
int i
Definition: decode_rs.h:71
grid_tg
float * grid_tg[4]
Definition: calc_par.c:92
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
MODISA
#define MODISA
Definition: sensorDefs.h:19
TAU_550_LOW_BOUND
#define TAU_550_LOW_BOUND
Definition: par_utils.h:18
VIIRSJ1
#define VIIRSJ1
Definition: sensorDefs.h:37
TAUCONS_LOW
#define TAUCONS_LOW
Definition: par_utils.h:26
EINSTEIN
#define EINSTEIN
Definition: par_utils.h:22
par_planar_b_inst
float * par_planar_b_inst
Definition: calc_par.c:611
par_utils.h
interp_as_taulow
float interp_as_taulow(float csz, float tau)
Definition: par_utils.c:1258
calc_scalar_par_mean_cosine
void calc_scalar_par_mean_cosine(l2str *l2rec, int ip, float par_above, float par_c, float *scalar_par, float *mean_cosine)
Definition: calc_par.c:160


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster