ocssw V2020
fluorescence.c
Go to the documentation of this file.
1 //---------------------------------------------------------------------
2 // fluorescence.c - functions to support chlorophyll fluorescence.
3 //
4 // Algorithm Reference:
5 // Behrenfeld, M.J., T.K. Westberry, E.S. Boss, R.T. O'Malley, D.A.
6 // Siegel, J.D. Wiggert, B.A. Franz, C.R. McClain, G.C. Feldman, S.C.
7 // Doney, J.K. Moore, G. Dall'Olmo, A. J. Milligan, I. Lima, and N.
8 // Mahowald (2009). Satellite-detected fluorescence reveals global
9 // physiology of ocean phytoplankton, Biogeosci., 6, 779-794.
10 //
11 // Implementation: B. Franz, Oct 2009.
12 //---------------------------------------------------------------------
13 
14 #include "l12_proto.h"
15 
16 #define PARW1 400
17 #define PARW2 700
18 #define PARWN (PARW2-PARW1)+1
20 static float fqymin = 0.0;
21 static float fqymax = 0.3;
22 static float flhmin = 0.0;
23 
24 /*---------------------------------------------------------------------*/
25 /* fsat_modis - normalized fluorescence line height for each pixel */
26 /* fsat output is in radiance units (mW/cm^2/um/sr) */
27 /*---------------------------------------------------------------------*/
28 
35 void fsat_modis(l2str *l2rec, float flh[]) {
36  static int32_t ib667, ib678, ib748;
37  static int firstCall = 1;
38 
39  int32_t ip, ipb;
40  float nLw1;
41  float nLw2;
42  float nLw3;
43  float bias = input->flh_offset;
44 
45  l1str *l1rec = l2rec->l1rec;
46  filehandle *l1file = l1rec->l1file;
47 
48  if (firstCall) {
49  firstCall = 0;
50  ib667 = windex(667., l1file->fwave, l1file->nbands);
51  ib678 = windex(678., l1file->fwave, l1file->nbands);
52  ib748 = windex(748., l1file->fwave, l1file->nbands);
53  }
54 
55  for (ip = 0; ip < l1rec->npix; ip++) {
56 
57  flh[ip] = BAD_FLT;
58 
59  ipb = l1file->nbands * ip;
63  nLw1 = l2rec->nLw[ipb + ib667];
64  nLw2 = l2rec->nLw[ipb + ib678];
65  nLw3 = l2rec->nLw[ipb + ib748];
66 
71  if (l1rec->mask[ip] || nLw1 < -0.01 || nLw2 < -0.01 || nLw3 < -0.01) {
72  l1rec->flags[ip] |= PRODFAIL;
73  continue;
74 
75  } else {
76 
80  flh[ip] = nLw2 - (70. / 81.) * nLw1 - (11. / 81.) * nLw3;
81 
82  // base = nLw3 + (nLw1 - nLw3) * ((748.0 - 678.0) / (748.0 - 667.0));
83  // flh[ip] = nLw2 - base - input->flh_offset;
84 
88  flh[ip] -= bias;
89  }
90  }
91 }
92 
93 
94 /*---------------------------------------------------------------------*/
95 /* flh_modis - fluorescence line height for each pixel in rec */
96 /* flh output is in radiance units (mW/cm^2/um/sr) */
97 /* OLD - SHOULD GO AWAY */
98 /*---------------------------------------------------------------------*/
99 
107 void flh_modis(l2str *l2rec, float flh[]) {
108  static int32_t ib667, ib678, ib748;
109  static int firstCall = 1;
110 
111  int32_t ip, ipb;
112  float Lw1;
113  float Lw2;
114  float Lw3;
115  float base;
116  float bias = input->flh_offset;
117 
118  l1str *l1rec = l2rec->l1rec;
119  filehandle *l1file = l1rec->l1file;
120 
121  if (firstCall) {
122  firstCall = 0;
123  ib667 = windex(667., l1file->fwave, l1file->nbands);
124  ib678 = windex(678., l1file->fwave, l1file->nbands);
125  ib748 = windex(748., l1file->fwave, l1file->nbands);
126  }
127 
128  for (ip = 0; ip < l1rec->npix; ip++) {
129 
130  flh[ip] = BAD_FLT;
131 
132  ipb = l1file->nbands * ip;
136  Lw1 = l2rec->Lw[ipb + ib667];
137  Lw2 = l2rec->Lw[ipb + ib678];
138  Lw3 = l2rec->Lw[ipb + ib748];
139 
140  // scattering-angle-dependent spectral difference in aerosol model phase
141  // functions are introducing geometric artifacts
142  // replace aerosol difference between red bands with simple model
143  // rho_a[678] = rho_a[667]*(678/667)^-0.65
144 
145  //La = l2rec->La[ipb+ib678]/l2rec->t_sen[ipb+ib678]; // original aerosol
146  //Lap = 0.99*l2rec->La[ipb+ib667]/l2rec->Fo[ib667]*l2rec->Fo[ib678]/l2rec->t_sen[ipb+ib678];
147  //Lw2 = Lw2 + La - Lap;
152  if (l1rec->mask[ip] || Lw1 < -0.01 || Lw2 < -0.01 || Lw3 < -0.01) {
153  l1rec->flags[ip] |= PRODFAIL;
154  continue;
155 
156  } else {
157 
162  base = Lw3 + (Lw1 - Lw3) * ((748.0 - 678.0) / (748.0 - 667.0));
163  flh[ip] = Lw2 - base;
164 
168  flh[ip] -= bias;
169 
170  if (flh[ip] < flhmin) {
171  flh[ip] = 0.0;
172  }
173  }
174  }
175 }
176 
177 /*---------------------------------------------------------------------*/
178 /* fqy - fluorescence quantum yield */
179 /*---------------------------------------------------------------------*/
180 
188 void fqy_modis(l2str *l2rec, float fqy[]) {
189  static float badval = BAD_FLT;
190  static int firstCall = 1;
191  static float *ipar;
192  static float *fsat;
193 
194  int32_t ip;
195 
196  l1str *l1rec = l2rec->l1rec;
197 
198  if (firstCall) {
199 
200  firstCall = 0;
201 
202  if ((ipar = (float *) calloc(l1rec->npix, sizeof (float))) == NULL) {
203  printf("-E- %s line %d: Unable to allocate space for ipar.\n",
204  __FILE__, __LINE__);
205  exit(1);
206  }
207  if ((fsat = (float *) calloc(l1rec->npix, sizeof (float))) == NULL) {
208  printf("-E- %s line %d: Unable to allocate space for fsat.\n",
209  __FILE__, __LINE__);
210  exit(1);
211  }
212  }
213 
214  // Compute ipar and fsat at all pixels
215 
216  get_ipar(l2rec, ipar);
217  fsat_modis(l2rec, fsat);
218 
219  // Compute fluorescence quantum yield for each pixel
220 
221  for (ip = 0; ip < l1rec->npix; ip++) {
222 
223  fqy[ip] = badval;
224 
225  if (!l1rec->mask[ip] && l2rec->chl[ip] > 0.0 && fsat[ip] >= 0.0) {
226 
227  // simple form from Behrenfeld, A12
228  fqy[ip] = 0.01 * fsat[ip] / (0.0302 * l2rec->chl[ip]) * ipar[ip]
229  * 1e6 / 1590;
230 
231  if (!isfinite(fqy[ip])) {
232  fqy[ip] = badval;
233  l1rec->flags[ip] |= PRODFAIL;
234  } else if (fqy[ip] < fqymin || fqy[ip] > fqymax) {
235  l1rec->flags[ip] |= PRODWARN;
236  }
237 
238  } else {
239  fqy[ip] = badval;
240  l1rec->flags[ip] |= PRODFAIL;
241  }
242  }
243 }
244 
251 void get_fqy(l2str *l2rec, float fqy[]) {
252  switch (l2rec->l1rec->l1file->sensorID) {
253  case MODISA:
254  case MODIST:
255  fqy_modis(l2rec, fqy);
256  break;
257  default:
258  printf("No fluorescence algorithm available for this sensor.\n");
259  exit(1);
260  break;
261  }
262 }
263 
269 void get_flh(l2str *l2rec, float flh[]) {
270  switch (l2rec->l1rec->l1file->sensorID) {
271  case MODISA:
272  case MODIST:
273  flh_modis(l2rec, flh);
274  break;
275  default:
276  printf("No fluorescence algorithm available for this sensor.\n");
277  exit(1);
278  break;
279  }
280 }
281 
289 void get_fsat(l2str *l2rec, float flh[]) {
290  switch (l2rec->l1rec->l1file->sensorID) {
291  case MODISA:
292  case MODIST:
293  fsat_modis(l2rec, flh);
294  break;
295  default:
296  printf("No fluorescence algorithm available for this sensor.\n");
297  exit(1);
298  break;
299  }
300 }
301 
302 // *************** old stuff *******************
303 
304 /* we want flhq in same units as arp (Ein/m^2/s) */
305 /* flh is in mW/cm^2/um/sr / 100 = W/m^2/nm/sr */
306 /* hc = 1986.5e-19 J*nm/photon, Lambda 676.7 nm */
307 /* radiance->iradiance = 4*pi */
308 /* line height -> area = 43.38 nm */
309 /* 6.023e23 per mole Quanta */
310 
311 /* Notes from modcol (anly8dbl.f90)
312  ! 4*$PI - convert radiance to scalar irradiance
313  ! energy of 1 photon = hc/lambda
314  ! h = 6.6261e-34 [Js]; c = 299,792,458 [m/s] (*10^9 for nm/s)
315  ! energy of one photon = (1986.5*10^-19)/Lambda [W s]
316  ! where lambda is [nm]
317  !
318  ! FLHQ units are quanta/m^2/s
319  ! to convert FLH into total fluorescence under the fluorescence
320  ! curve :
321  ! The fluorescence curve is described as a gaussian curve with
322  ! center at 683 nm and half maximum emission of 25 nm (from
323  ! Collins et al. 1985
324  ! If L683 (1 nm width) = 1 then the area under the curve= 26.61
325  ! Band14 center= 677 bandwidth= 10 Fluorescence signal read= 8.61
326  ! Band13 center= 665 bandwidth= 10 Fluorescence signal read= 2.86
327  ! Band15 center= 747 bandwidth= 10 Fluorescence signal read= 1.55e-6
328  ! All band readings correspond to a 10 nm bandwidth signal. For this
329  ! reason the signals must be divided by 10 to get the reading per nm
330  !
331  ! Because band 13 is affected by the fluorescence signal then there
332  ! is a contribution of this signal equal to 0.2478 to the baseline
333  ! correction.
334  !
335  ! The conversion factor of FLH into area under the curve can be
336  ! computed: area under the curve / (Band14 read - baseline contrib)
337  ! factor = 26.61 /(0.861 - 0.2478) = 43.38
338  ! rescale: 43.38 (nm) -> 0.04338 (um)
339  !
340  ! FLHQ = FLH * Lambda/(hc) * (radiance->irradiance)*(FLH->area)
341  ! FLHQ is in quanta m-2 s-1 and ARP is in Einstein m-2 s-1
342  ! 1 Einstein = 1 mol quanta = 6.023e23 quanta
343  */
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
void get_fqy(l2str *l2rec, float fqy[])
Definition: fluorescence.c:251
#define NULL
Definition: decode_rs.h:63
read l1rec
#define PRODWARN
Definition: l2_flags.h:11
float32 base
Definition: maplists.h:106
#define MODIST
Definition: sensorDefs.h:18
void fqy_modis(l2str *l2rec, float fqy[])
Definition: fluorescence.c:188
#define PRODFAIL
Definition: l2_flags.h:39
instr * input
void flh_modis(l2str *l2rec, float flh[])
Definition: fluorescence.c:107
#define BAD_FLT
Definition: jplaeriallib.h:19
void get_fsat(l2str *l2rec, float flh[])
Definition: fluorescence.c:289
void get_ipar(l2str *l2rec, float ipar[])
Definition: ipar.c:18
void get_flh(l2str *l2rec, float flh[])
Definition: fluorescence.c:269
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:71
#define MODISA
Definition: sensorDefs.h:19
void fsat_modis(l2str *l2rec, float flh[])
Definition: fluorescence.c:35