OB.DAAC Logo
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, in integration
48  * * Added Abar 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 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 <stdio.h>
80 #include <stdlib.h>
81 #include <math.h>
82 #include "l2_struc.h"
83 #include "par_utils.h"
84 #include <timeutils.h>
85 #include "l12_proto.h"
86 
87 float calc_par(l2str *l2rec, int ip, int nbands, float *Lt, float taua,
88  float angstrom, float *wl, float *fo, float *ko3, float *taumolbar) {
89 
90  int i, band;
91  float r[nbands], csolz, transg;
92  float weight[nbands];
93  float rp[nbands], transa, sa, ps0, del;
94  float asymfac, beta, tau[nbands], taumol[nbands], tauaer[nbands];
95  float gamma, cos2x, pmol;
96  float csenz, transav;
97 
98  float par0, integtrans, delt, saavg;
99  float tmstep, sz, cossz, tg, td;
100  float sola;
101  float asavg, func[11], ra[nbands], rc[nbands];
102 
103  float sumSa, sumrc, sumtdir, sumtdif, denom, tatemp, ravgsat;
104  float tdir, tdif;
105  float q, q4, taucons, q4tau, f, asat, asbar, abar, acldmax;
106  float fcsolz, fcsenz, fcossz;
107  float tg_oz, tg_wv;
108  float trise, tset;
109  float par;
110 
111  int widx490, widx510;
112  float ang500, ta500, slope;
113 
114  int cloudy = FALSE;
115 
116  static float kavg_oz = 5.2e-05;
117  static float kavg_wv = 0.002;
118 
119  // need some memory allocation
120  float *aerphasefunc;
121  float *omega;
122  float *modelAngstrom;
123 
124  l1str *l1rec = l2rec->l1rec;
125 
126  int16_t year, month, mday, doy;
127  double sec;
128  unix2yds(l1rec->scantime, &year, &doy, &sec);
129  unix2ymds(l1rec->scantime, &year, &month, &mday, &sec);
130 
131  aerphasefunc = (float*) allocateMemory(nbands * sizeof (float),
132  "aerphasefunc");
133  omega = (float*) allocateMemory(nbands * sizeof (float), "omega");
134  modelAngstrom = (float*) allocateMemory(nbands * sizeof (float),
135  "modelAngstrom");
136 
137  int asat_is_max = 0;
138 
139  /****************************************************************************
140  * Perform data checks; assign defaults if necessary: *
141  * If data invalid and uncorrectable, PAR set to BAD_FLT missing value. *
142  ***************************************************************************/
143 
144  if (l1rec->solz[ip] < 0.0f || l1rec->solz[ip] > 90.0f) {
145  printf("Error: subroutine CalcPAR. computation failed.\n");
146  printf("SolZen (outside valid range 0-90) = %f\n", l1rec->solz[ip]);
147  par = BAD_FLT;
148  return par;
149  }
150 
151  if (l1rec->senz[ip] < 0.0f || l1rec->senz[ip] > 90.0f) {
152  printf("Error: subroutine CalcPAR. computation failed.\n");
153  printf("ViewZen (outside valid range 0-90) = %f\n", l1rec->senz[ip]);
154  par = BAD_FLT;
155  return par;
156  }
157 
158  if (l1rec->delphi[ip] < -180.0f || l1rec->delphi[ip] > 180.0f) {
159  printf("Error: subroutine CalcPAR. computation failed.\n");
160  printf("delphi (outside valid range -180:180) = %f\n",
161  l1rec->delphi[ip]);
162  par = BAD_FLT;
163  return par;
164  }
165 
166  float dobson = l1rec->oz[ip];
167 
168  csolz = l1rec->csolz[ip];
169  csenz = l1rec->csenz[ip];
170  cos2x = powf(cosf((l1rec->scattang[ip] * M_PI / 180.0)), 2.0);
171 
172  if (l1rec->scattang[ip] < 0.0 || l1rec->scattang[ip] > 180.0) {
173  printf("Error: subroutine CalcPAR. computation failed.\n");
174  printf("scatang (outside valid range 0-180) = %f\n",
175  l1rec->scattang[ip]);
176  printf("scatang = -(cossz*zosvz)+(sinsz*sinvz*cosra)\n");
177  par = BAD_FLT;
178  return par;
179  }
180 
181  if (taua <= 0.0f) {
182  taua = 0.1f;
183  }
184 
185  float surfpress = l1rec->pr[ip];
186  if (surfpress <= 0.0f) {
187  surfpress = 1013.15f;
188  }
189 
190  /* Need to use conventional Angstrom exponent.
191  if it's missing GetAerPhase should work and return a default value */
192 
193  GetAerPhase(l2rec, ip, nbands, angstrom, aerphasefunc, omega,
194  modelAngstrom);
195 
196  /****************************************************************************
197  * Angstrom exponents are conventionally > 0, but actually < 0 *
198  ****************************************************************************/
199 
200  angstrom = -1 * angstrom;
201 
202 
203  if (dobson < 0.0f) {
204  dobson = EstimateDobson(year, month, mday, l1rec->lat[ip]);
205  }
206 
207  float watvap = l1rec->wv[ip];
208  if (watvap < 0.0f) {
209  watvap = EstimateWatVap(year, month, mday, l1rec->lat[ip]);
210  }
211 
212  /*****************************************************************************
213  * Transform TOA radiances into reflectances: *
214  ****************************************************************************/
215  // TODO: consider passing in rhot instead of Lt
216  float Airmass = (1.f / csolz) + (1.f / csenz);
217 
218  for (band = 0; band < nbands; band++) {
219  r[band] = (M_PI * (Lt[band])) / (fo[band] * l1rec->fsol * csolz);
220  // Compute gaseous transmittance:
221  transg = expf((-ko3[band] * dobson * Airmass));
222  rp[band] = r[band] / transg;
223  }
224 
225  /****************************************************************************
226  * Compute reflectances of the cloud-surface subsystem: *
227  ***************************************************************************/
228 
229  /* Setup for transmittance equation:*/
230 
231  ps0 = 1013.5f;
232  del = 0.0095f;
233 
234  pmol = (2.0f * (1.0f - del) * 0.75f * (1.0f + cos2x) + 3.0f * del)
235  / (2.0f + del);
236 
237  asymfac = 0.6666667f;
238  beta = 0.5f * (1.0f + asymfac);
239  gamma = 1.0f - asymfac;
240 
241  /******************************************************************************
242  * Compute diffuse atmospheric transmittance: *
243  ******************************************************************************/
244  sumSa = 0.0f;
245  sumrc = 0.0f;
246  sumtdir = 0.0f;
247  sumtdif = 0.0f;
248  denom = 0.0f;
249 
250  weight[0] = (wl[0] - 400.0) + (wl[1] - wl[0]) / 2.0;
251  for (band = 1; band < nbands - 1; band++) {
252  weight[band] = (wl[band] - wl[band - 1]) / 2.0
253  + (wl[band + 1] - wl[band]) / 2.0;
254  }
255  weight[nbands - 1] = (700.0 - wl[nbands - 1])
256  + (wl[nbands - 1] - wl[nbands - 2]) / 2.0;
257 
258  for (band = 0; band < nbands; band++) {
259  taumol[band] = taumolbar[band] * (surfpress / ps0);
260 
261  tauaer[band] = taua
262  * powf(wl[band] / (float) input->aer_wave_long,
263  modelAngstrom[band]);
264 
265  tau[band] = taumol[band] + gamma * tauaer[band];
266 
267  ra[band] = ((taumol[band] * pmol)
268  + (omega[band] * tauaer[band] * aerphasefunc[band]))
269  / (4.0f * csolz * csenz);
270 
271  /* Compute diffuse atmospheric transmittance:*/
272 
273  transa = expf(-(taumol[band] + tauaer[band]) / csolz)
274  * expf((0.52f * taumol[band] + beta * tauaer[band]) / csolz);
275  transav = expf(-(taumol[band] + tauaer[band]) / csenz)
276  * expf((0.52f * taumol[band] + beta * tauaer[band]) / csenz);
277 
278  // Get direct and diffuse transmittances
279  tdir = expf(-tau[band] / csolz);
280 
281  tdif = (expf(-tau[band] / csolz))
282  * (expf((0.52f * taumol[band] + beta * tauaer[band]) / csolz)
283  - 1.f);
284  sumtdir += (tdir * fo[band]);
285  sumtdif += (tdif * fo[band]);
286 
287  /* Compute spherical albedo of the atmosphere:*/
288  sa = expf(-(taumol[band] + gamma * tauaer[band]))
289  * (0.92f * taumol[band] + gamma * tauaer[band]);
290 
291  sumSa += sa * fo[band] * weight[band];
292  denom += fo[band] * weight[band];
293 
294  rc[band] = (rp[band] - ra[band])
295  / ((transa * transav) + (sa * (rp[band] - ra[band])));
296  sumrc += (rc[band] * fo[band] * weight[band]);
297  }
298 
299  //Derive tauaer(500) from (469,555)
300  // indexes into modelAngstrom need to use bindex
301  switch (l1rec->l1file->sensorID) {
302  case MODIST:
303  case MODISA:
304  widx490 = windex(469., wl, nbands);
305  widx510 = windex(555., wl, nbands);
306  break;
307  case VIIRSN:
308  case VIIRSJ1:
309  widx490 = windex(490., wl, nbands);
310  widx510 = windex(555., wl, nbands);
311  break;
312  case OCTS:
313  widx490 = windex(490., wl, nbands);
314  widx510 = windex(520., wl, nbands);
315  break;
316  default:
317  widx490 = windex(490., wl, nbands);
318  widx510 = windex(510., wl, nbands);
319  break;
320  }
321 
322  slope = (modelAngstrom[widx510] - modelAngstrom[widx490])
323  / (wl[widx510] - wl[widx490]);
324  ang500 = modelAngstrom[widx490] + slope * (wl[widx510] - wl[widx490]);
325  ta500 = taua * powf((500.0 / (float) input->aer_wave_long), ang500);
326 
327  /******************************************************************************
328  * Compute daily average PAR: *
329  ******************************************************************************/
330 
331  /* Set up terms constant in the integration:*/
332  // Sa_bar:
333  saavg = sumSa / denom;
334  // R_bar:
335  ravgsat = sumrc / denom;
336 
337  /******************************************************************************
338  *Get simple bi-directional reflectance factor: *
339  ******************************************************************************/
340 
341  q = 1.0f / (3.0f * (1.0f - 0.853f));
342  q4 = 4.0f * q;
343 
344  taucons = 15.0f;
345  q4tau = q4 / (taucons + q4);
346 
347  fcsenz = (3.f / 7.f) * (1.0f + 2.0f * csenz);
348  fcsolz = (3.f / 7.f) * (1.0f + 2.0f * csolz);
349 
350  f = (1.0f - q4tau * fcsolz)
351  / (0.49f * ((1.f + 4.f * csolz * csenz) / (csolz + csenz))
352  - q4tau * fcsolz * fcsenz);
353 
354  /******************************************************************************
355  * Get As_bar: *
356  ******************************************************************************/
357  asbar = interp_as_taulow(csolz, ta500);
358 
359  /******************************************************************************
360  * Get Albedo from R_bar: *
361  ******************************************************************************/
362  asat = (ravgsat - asbar) * f + asbar;
363 
364  /******************************************************************************
365  * Test Asat to see if it's too large; adjust as necessary: *
366  ******************************************************************************/
367 
368  acldmax = 1.0f - (q4 / (300.0f + q4)) * fcsolz;
369 
370  if ((asat - asbar) >= acldmax) {
371  asat_is_max = 1;
372  taucons = 300.0f;
373  q4tau = q4 / (taucons + q4);
374  } else {
375  asat_is_max = 0;
376  taucons = q4 * ((fcsolz / (1.0f - (asat - asbar))) - 1.0f);
377 
378  /***************************************************************************
379  * Compute new Asat if needed: *
380  ***************************************************************************/
381 
382  if (taucons > 15.0f) {
383  q4tau = q4 / (taucons + q4);
384 
385  f = (1.0f - q4tau * fcsolz)
386  / (0.49f * ((1.f + 4.f * csolz * csenz) / (csolz + csenz))
387  - q4tau * fcsolz * fcsenz);
388  asat = (ravgsat - asbar) * f + asbar;
389 
390  /* New Test:*/
391 
392  if ((asat - asbar) >= acldmax) {
393  asat_is_max = 1;
394  taucons = 300.0f;
395  q4tau = q4 / (taucons + q4);
396  }
397  } else {
398  taucons = 15.0f;
399  q4tau = q4 / (taucons + q4);
400  }
401  }
402 
403  // Test to estimate if pixel is clear or cloudy:
404  cloudy = FALSE;
405  if (asat > asbar)
406  cloudy = TRUE;
407 
408  /*****************************************************************************
409  * Get parameters used in integration time-steps: *
410  *****************************************************************************/
411 
412  /* Get the rise and set time for this day (for integral):*/
413 
414  triseset(doy, l1rec->lon[ip], l1rec->lat[ip], &trise, &tset);
415 
416  /* Set up for 10 intervals*/
417  delt = (tset - trise) / 10.0f;
418 
419  /***************************************************************************
420  * Get parameters used in integration time-steps: *
421  ***************************************************************************/
422 
423  for (i = 0; i < 11; i++) {
424  tmstep = trise + i * delt;
425  int intyear, intday;
426  intyear = year;
427  intday = doy;
428  sunangs_(&intyear, &intday, &tmstep, l1rec->lon + ip,
429  l1rec->lat + ip, &sz, &sola);
430  cossz = cosf(sz * (float) M_PI / 180.0f);
431  fcossz = (3.f / 7.f) * (1.0f + 2.0f * cossz);
432 
433  /***************************************************************************
434  * Get gaseous transmittance term: *
435  ***************************************************************************/
436  if (sz >= 89.5f) {
437  tg_oz = 0.0f;
438  tg_wv = 0.0f;
439  } else {
440  //kavg_oz (5.2E-5) is valid for dobson in true dobson units, not the ozone units in l2gen
441  tg_oz = expf(-kavg_oz * powf((dobson * 1000. / cossz), 0.99f));
442  tg_wv = expf(-kavg_wv * powf((watvap / cossz), 0.87f));
443  }
444 
445  tg = tg_oz * tg_wv;
446  /***************************************************************************
447  * Get average diffuse atmospheric transmittance and As_bar terms: *
448  ***************************************************************************/
449  td = 0.0f;
450 
451  if (sz < 89.5f) {
452  for (band = 0; band < nbands; band++) {
453  tatemp = expf(-(taumol[band] + tauaer[band]) / cossz)
454  * expf(
455  (0.52f * taumol[band] + beta * tauaer[band])
456  / cossz);
457  td += tatemp * fo[band] * weight[band];
458  }
459  td /= denom;
460  }
461 
462  if (cossz < 0.05) {
463  cossz = 0.05;
464  }
465  if (cloudy) {
466  asavg = interp_as_tauhigh(cossz);
467  } else {
468  asavg = interp_as_taulow(cossz, ta500);
469  }
470 
471  /*****************************************************************************
472  *Get A_bar term: *
473  *****************************************************************************/
474 
475  if (asat_is_max) {
476  abar = 1.0f - q4tau * fcossz;
477  } else {
478  abar = asat * (1.0f - q4tau * fcossz) / (1.0f - q4tau * fcsolz);
479  }
480 
481  /*new test:*/
482  if (abar >= 1.0f) {
483  taucons = 300.0f;
484  q4tau = q4 / (taucons + q4);
485  abar = 1.0f - q4tau * fcossz;
486  }
487 
488  /****************************************************************************
489  * Get "function" integrand term: *
490  ****************************************************************************/
491 
492  if (abar <= asavg) {
493  func[i] = cossz * tg * td / (1.0 - abar * saavg);
494  } else {
495  func[i] = cossz * tg * td * (1.0f - abar) / (1.0f - asavg)
496  / (1.0f - abar * saavg);
497  }
498  }
499 
500  /* Use trapezoidal rule for integration:*/
501 
502  integtrans = 0.0f;
503 
504  for (i = 0; i < 10; i++) {
505  integtrans += delt * (func[i] + func[i + 1]) / 2.0f;
506  }
507 
508  /* Finally, determine the PAR:*/
509  // par0 determined as the mean value [400-700nm] for Thuillier 2003
510  // from run/data/common/Thuillier_F0.dat
511  par0 = 176.323f;
512  par = par0 * l1rec->fsol * integtrans / 24.0f;
513 
514  free(aerphasefunc);
515  free(omega);
516  free(modelAngstrom);
517 
518  return par;
519 }
data_t q
Definition: decode_rs.h:74
int r
Definition: decode_rs.h:73
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:99
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
#define FALSE
Definition: rice.h:164
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:87
read l1rec
#define VIIRSN
Definition: sensorDefs.h:23
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
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
#define TRUE
Definition: rice.h:165
#define MODIST
Definition: sensorDefs.h:18
void GetAerPhase(l2str *l2rec, int ip, int32_t nbands, float angstrom, float *aerphasefunc, float *omega, float *modelAngstrom)
Definition: par_utils.c:51
instr * input
void sunangs_(int *, int *, float *, float *, float *, float *, float *)
double precision function f(R1)
Definition: tmd.lp.f:1454
float interp_as_tauhigh(float csz)
Definition: par_utils.c:769
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
subroutine func(x, conec, n, bconecno, bn, units, u, inno, i, outno, o, Input, Targ, p, sqerr)
Definition: ffnet.f:287
float32 slope[]
Definition: l2lists.h:30
float rc(float x, float y)
#define M_PI
Definition: pml_iop.h:15
void unix2yds(double usec, short *year, short *day, double *secs)
void triseset(int32_t jday, float xlon, float xlat, float *trise, float *tset)
Definition: par_utils.c:529
float gamma
Definition: color_dtdb.py:447
float EstimateWatVap(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:381
#define BAD_FLT
Definition: jplaeriallib.h:19
#define OCTS
Definition: sensorDefs.h:14
double beta
int32_t nbands
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
float EstimateDobson(int32_t year, int32_t month, int32_t day, float lat)
Definition: par_utils.c:250
int i
Definition: decode_rs.h:71
#define MODISA
Definition: sensorDefs.h:19
#define VIIRSJ1
Definition: sensorDefs.h:37
float interp_as_taulow(float csz, float tau)
Definition: par_utils.c:661