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
get_ctht.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include "l12_proto.h"
4 #include "l2prod.h"
5 #include "l1.h"
6 #include "get_ctht.h"
7 #include "sensorInfo.h"
8 #define C_TO_K 273.15
9 
10 ctht_lins_str ctht_lins;
11 /*
12  file get_ctht.c contains most of the cloud top height source code
13 */
14 float *get_ctht_lin( l2str *l2rec, int prodnum )
15 /*
16 get_ctht_lin - get the cloud top height info for a current scan
17  Returns - float * - the record of the scan's product data
18 
19  Parameters: (in calling order)
20  Type Name I/O Description
21  ---- ---- --- -----------
22  l2str * l2rec I l2 record of data
23  int prodnum I Product ID to return
24  Note that if prodnum = -1
25  we won't return any product
26  This is when the code is called
27  to prepare for getting cmp
28  cloud values
29 
30  W. Robinson, SAIC, 22 Jun 2023
31 */
32  {
33  static int32_t cur_ctht_scan = -1;
34  // being managed: 0 - not filled, 1 filled
35  // (may need status for '-1 line or nlin + 1)
36  static int32_t firstin = 0;
37  static float *outbuf;
38  static int32_t *outint; // so ints will get treated right
39  int32_t cur_scan = l2rec->l1rec->iscan;
40  int32_t npix = l2rec->l1rec->npix;
41  int32_t clin, ipix, ntyp = 2, phase_v;
42  ctht_lins.nrec = 3;
43 
44  outint = (int32_t *)outbuf;
45 
46  if( firstin == 0 )
47  {
48  firstin = 1;
49  if( ( outbuf = malloc( npix * ntyp * sizeof( float ) ) ) == NULL )
50  {
51  printf( "%s, %d: E - Unable to allocate the ctht output buffer\n",
52  __FILE__, __LINE__ );
53  exit(1);
54  }
55  }
56 
57  clin = ctht_lins.nrec / 2;
58  if( cur_ctht_scan != cur_scan )
59  {
60  cur_ctht_scan = cur_scan;
61  if( comp_ctht( l2rec ) != 0 )
62  {
63  printf( "%s, %d: E - CTH computation failure\n", __FILE__, __LINE__ );
64  exit(1);
65  }
66  }
67  /* output the specified data */
68  // NOTE this should be in the center line of the 'queue'
69  switch( prodnum )
70  {
71  case CAT_Cld_p:
72  memcpy( outbuf, ctht_lins.ct[clin]->ctp, npix * sizeof( float ) );
73  return outbuf;
74  break;
75  case CAT_Cld_t:
76  memcpy( outbuf, ctht_lins.ct[clin]->ctt, npix * sizeof( float ) );
77  return outbuf;
78  break;
79  case CAT_Cld_h:
80  memcpy( outbuf, ctht_lins.ct[clin]->cth, npix * sizeof( float ) );
81  return outbuf;
82  break;
83  case CAT_cth_cod:
84  for( ipix = 0; ipix < npix; ipix++ )
85  {
86  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
87  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
88  ctht_lins.ct[clin]->cth_cod_all[ phase_v + ntyp * ipix ];
89  }
90  return outbuf;
91  break;
92  case CAT_cth_cod_all:
93  memcpy( outbuf, ctht_lins.ct[clin]->cth_cod_all,
94  npix * ntyp * sizeof( float ) );
95  return outbuf;
96  break;
97  case CAT_cth_alb_all:
98  memcpy( outbuf, ctht_lins.ct[clin]->cth_alb_all,
99  npix * ntyp * sizeof( float ) );
100  return outbuf;
101  break;
102  case CAT_cth_alb:
103  for( ipix = 0; ipix < npix; ipix++ )
104  {
105  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
106  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
107  ctht_lins.ct[clin]->cth_alb_all[ phase_v + ntyp * ipix ];
108  }
109  return outbuf;
110  break;
111  case CAT_cth_dcod:
112  for( ipix = 0; ipix < npix; ipix++ )
113  {
114  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
115  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
116  ctht_lins.ct[clin]->dcod[ phase_v + ntyp * ipix ];
117  }
118  return outbuf;
119  break;
120  case CAT_cth_dlcod:
121  for( ipix = 0; ipix < npix; ipix++ )
122  {
123  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
124  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
125  ctht_lins.ct[clin]->dlcod[ phase_v + ntyp * ipix ];
126  }
127  return outbuf;
128  break;
129  case CAT_cth_dcth:
130  for( ipix = 0; ipix < npix; ipix++ )
131  {
132  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
133  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
134  ctht_lins.ct[clin]->dcth[ phase_v + ntyp * ipix ];
135  }
136  return outbuf;
137  break;
138  case CAT_cth_dctp:
139  for( ipix = 0; ipix < npix; ipix++ )
140  {
141  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
142  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
143  ctht_lins.ct[clin]->dctp[ phase_v + ntyp * ipix ];
144  }
145  return outbuf;
146  break;
147  case CAT_cth_dctt:
148  for( ipix = 0; ipix < npix; ipix++ )
149  {
150  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
151  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
152  ctht_lins.ct[clin]->dctt[ phase_v + ntyp * ipix ];
153  }
154  return outbuf;
155  break;
156  case CAT_cth_dalb:
157  for( ipix = 0; ipix < npix; ipix++ )
158  {
159  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
160  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
161  ctht_lins.ct[clin]->dalb[ phase_v + ntyp * ipix ];
162  }
163  return outbuf;
164  break;
165  case CAT_cth_phase:
166  memcpy( outbuf, ctht_lins.ct[clin]->ct_phase,
167  npix * sizeof( unsigned char ) );
168  return outbuf;
169  break;
170  case CAT_cth_lcod_all:
171  memcpy( outbuf, ctht_lins.ct[clin]->cth_lcod_all,
172  npix * ntyp * sizeof( float ) );
173  return outbuf;
174  break;
175  case CAT_cth_cth_all:
176  memcpy( outbuf, ctht_lins.ct[clin]->cth_all,
177  npix * ntyp * sizeof( float ) );
178  return outbuf;
179  break;
180  case CAT_cth_ctp_all:
181  memcpy( outbuf, ctht_lins.ct[clin]->ctp_all,
182  npix * ntyp * sizeof( float ) );
183  return outbuf;
184  break;
185  case CAT_cth_cost_all:
186  memcpy( outbuf, ctht_lins.ct[clin]->cost_ss,
187  npix * ntyp * sizeof( float ) );
188  return outbuf;
189  break;
190  case CAT_cth_acost_all:
191  memcpy( outbuf, ctht_lins.ct[clin]->acost,
192  npix * ntyp * sizeof( float ) );
193  return outbuf;
194  break;
195  case CAT_cth_iter_all:
196  memcpy( outint, ctht_lins.ct[clin]->nitr,
197  npix * ntyp * sizeof( int32_t) );
198  return outbuf;
199  break;
200  case CAT_cth_cth_raw_all:
201  memcpy( outbuf, ctht_lins.ct[clin]->cth_raw,
202  npix * ntyp * sizeof( float ) );
203  return outbuf;
204  break;
205  case CAT_cth_ctt_all:
206  memcpy( outbuf, ctht_lins.ct[clin]->ctt_all,
207  npix * ntyp * sizeof( float ) );
208  return outbuf;
209  break;
210  // back to all 2-d arrays that get extracted from the 3-d
211  case CAT_cth_lcod:
212  for( ipix = 0; ipix < npix; ipix++ )
213  {
214  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
215  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
216  ctht_lins.ct[clin]->cth_lcod_all[ phase_v + ntyp * ipix ];
217  }
218  return outbuf;
219  break;
220  case CAT_cth_cost:
221  for( ipix = 0; ipix < npix; ipix++ )
222  {
223  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
224  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
225  ctht_lins.ct[clin]->cost_ss[ phase_v + ntyp * ipix ];
226  }
227  return outbuf;
228  break;
229  case CAT_cth_acost:
230  for( ipix = 0; ipix < npix; ipix++ )
231  {
232  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
233  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
234  ctht_lins.ct[clin]->acost[ phase_v + ntyp * ipix ];
235  }
236  return outbuf;
237  break;
238  case CAT_cth_iter:
239  for( ipix = 0; ipix < npix; ipix++ )
240  {
241  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
242  outint[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_INT :
243  ctht_lins.ct[clin]->nitr[ phase_v + ntyp * ipix ];
244  }
245  return outbuf;
246  break;
247  case CAT_cth_cth_raw:
248  for( ipix = 0; ipix < npix; ipix++ )
249  {
250  phase_v = ctht_lins.ct[clin]->ct_phase[ipix];
251  outbuf[ipix] = ( phase_v == BAD_UBYTE ) ? BAD_FLT :
252  ctht_lins.ct[clin]->cth_raw[ phase_v + ntyp * ipix ];
253  }
254  return outbuf;
255  break;
256  case CAT_cth_dcod_all:
257  memcpy( outbuf, ctht_lins.ct[clin]->dcod,
258  npix * ntyp * sizeof( float ) );
259  return outbuf;
260  break;
261  case CAT_cth_dlcod_all:
262  memcpy( outbuf, ctht_lins.ct[clin]->dlcod,
263  npix * ntyp * sizeof( float ) );
264  return outbuf;
265  break;
266  case CAT_cth_dcth_all:
267  memcpy( outbuf, ctht_lins.ct[clin]->dcth,
268  npix * ntyp * sizeof( float ) );
269  return outbuf;
270  break;
271  case CAT_cth_dctp_all:
272  memcpy( outbuf, ctht_lins.ct[clin]->dctp,
273  npix * ntyp * sizeof( float ) );
274  return outbuf;
275  break;
276  case CAT_cth_dctt_all:
277  memcpy( outbuf, ctht_lins.ct[clin]->dctt,
278  npix * ntyp * sizeof( float ) );
279  return outbuf;
280  break;
281  case CAT_cth_dalb_all:
282  memcpy( outbuf, ctht_lins.ct[clin]->dalb,
283  npix * ntyp * sizeof( float ) );
284  return outbuf;
285  break;
286  case -1:
287  return NULL;
288  }
289  /* if this is reached, we need another CAT */
290  printf( "%s, %d: E - CTH computation could not find the catalog # %d\n",
291  __FILE__, __LINE__, prodnum );
292  exit(1);
293  }
294 
295 int32_t comp_ctht( l2str *l2rec )
296  /*
297  comp_ctht will run the cloud top height algorithm on the lines of data in
298  the l1que that have not been processed yet and keep a local copy of the lines
299 
300  Returns - int32_t 0 if all good, no error conditions yet
301 
302  Parameters: (in calling order)
303  Type Name I/O Description
304  ---- ---- --- -----------
305  l2str * l2rec I l2 record of data
306 
307  W. Robinson, SAIC, 22 Jun 2023
308  */
309  {
310  int32_t ctht_nrec = 3; // # lines in our ctht records, matches those
311  // in the cmp, see get_cmp.c
312  static int32_t *ctht_stat;
313  static int32_t ctht_init = 0;
314  int32_t ilin, il2, lin_in_que, scn_in_que, ntyp = 2;
315  extern l1qstr l1que;
316  static ctht_lin_str *ctht_sav;
317  int32_t npix = l2rec->l1rec->npix;
318 
319  if( input->proc_cloud != 1 ) {
320  printf( "E, %s, %d: For cloud processing, proc_cloud must be set to 1\n",
321  __FILE__, __LINE__ );
322  exit(1);
323  }
324 
325  if( l2rec->l1rec->anc_add == NULL ) {
326  printf( "E, %s, %d: For cloud processing, anc_profile_1, _2, _3 must be specified\n", __FILE__, __LINE__ );
327  exit(1);
328  }
329 
330  // until climatology sfcp, sfct is cleared up, we can't use climatology
331  if( strstr( input->met1, "climatology" ) != NULL ) {
332  printf( "E, %s, %d: Ancillary meteorological data from a climatology file cannot be used for cloud processing\n", __FILE__, __LINE__ );
333  printf( "Use daily files in MET1, 2, 3 l2gen inputs\n" );
334  exit(1);
335  }
336 
337  /* initialize the ctht line storage and lines */
338  if( ctht_init == 0 )
339  {
340  ctht_init = 1;
341  /* set up the ctht record status */
342  if( ( ( ctht_lins.ct = (ctht_lin_str **)malloc(
343  ctht_nrec * sizeof( ctht_lin_str * ) ) ) == NULL ) ||
344  ( ( ctht_stat = (int32_t *)malloc(
345  ctht_nrec * sizeof( int32_t ) ) ) == NULL ) )
346  {
347  printf( "%s - %d: Allocation of ctht record space failed\n",
348  __FILE__, __LINE__ );
349  exit(1);
350  }
351  for( ilin = 0; ilin < ctht_nrec; ilin++ )
352  {
353  if( ( ctht_lins.ct[ilin] =
354  (ctht_lin_str *) malloc( sizeof( ctht_lin_str) ) ) == NULL )
355  {
356  printf( "%s - %d: Allocation of ctht record space failed\n",
357  __FILE__, __LINE__ );
358  exit(1);
359  }
360  ctht_lins.ct[ilin]->iscan = -1;
361  ctht_lins.ct[ilin]->npix = npix;
362  if( ( ( ctht_lins.ct[ilin]->cth =
363  (float *) malloc(npix * sizeof(float) ) ) == NULL ) ||
364  ( ( ctht_lins.ct[ilin]->ctp =
365  (float *) malloc(npix * sizeof(float) ) ) == NULL ) ||
366  ( ( ctht_lins.ct[ilin]->ctt =
367  (float *) malloc(npix * sizeof(float) ) ) == NULL ) ||
368  ( ( ctht_lins.ct[ilin]->ct_phase =
369  (unsigned char *) malloc(npix * sizeof(unsigned char) ) ) == NULL )
370  || ( ( ctht_lins.ct[ilin]->cth_all =
371  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
372  ( ( ctht_lins.ct[ilin]->cth_raw =
373  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
374  ( ( ctht_lins.ct[ilin]->ctp_all =
375  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
376  ( ( ctht_lins.ct[ilin]->ctt_all =
377  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
378  ( ( ctht_lins.ct[ilin]->cth_lcod_all =
379  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
380  ( ( ctht_lins.ct[ilin]->cth_alb_all =
381  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
382  ( ( ctht_lins.ct[ilin]->cost_ss =
383  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
384  ( ( ctht_lins.ct[ilin]->acost =
385  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
386  ( ( ctht_lins.ct[ilin]->nitr =
387  (int32_t *) malloc(npix * ntyp * sizeof(int32_t) ) ) == NULL ) )
388  {
389  printf( "%s - %d: Allocation of ctht record space failed\n",
390  __FILE__, __LINE__ );
391  exit(1);
392  }
393  if( ( ( ctht_lins.ct[ilin]->cth_cod_all =
394  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
395  ( ( ctht_lins.ct[ilin]->oe_akdiag =
396  (float *) malloc(npix * 3 * ntyp * sizeof(float) ) ) == NULL ) ||
397  ( ( ctht_lins.ct[ilin]->dlcod =
398  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
399  ( ( ctht_lins.ct[ilin]->dcod =
400  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
401  ( ( ctht_lins.ct[ilin]->dcth =
402  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
403  ( ( ctht_lins.ct[ilin]->dctp =
404  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
405  ( ( ctht_lins.ct[ilin]->dctt =
406  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) ||
407  ( ( ctht_lins.ct[ilin]->dalb =
408  (float *) malloc(npix * ntyp * sizeof(float) ) ) == NULL ) )
409  {
410  printf( "%s - %d: Allocation of ctht record space failed\n",
411  __FILE__, __LINE__ );
412  exit(1);
413  }
414  }
415  }
416  /* OK, arrange the existant ctht records that we have and process any
417  records we don't have in comp_ctht_lin
418  */
419 
420  // go thru ctht lines
421  for( ilin = 0; ilin < ctht_nrec; ilin++ )
422  {
423  ctht_stat[ilin] = 0;
424  }
425  for( ilin = 0; ilin < ctht_nrec; ilin++ )
426  {
427  lin_in_que = ilin + l1que.nq / 2 - ctht_nrec / 2;
428  // if lin_in_que < 0 we're in trouble...
429  scn_in_que = l1que.r[lin_in_que].iscan;
430  // search for the l1rec match line in ctht
431  for( il2 = 0; il2 < ctht_nrec; il2++ )
432  {
433  if( ctht_lins.ct[il2]->iscan == scn_in_que )
434  {
435  // switch addresses
436  ctht_sav = ctht_lins.ct[ilin];
437  ctht_lins.ct[ilin] = ctht_lins.ct[il2];
438  ctht_lins.ct[il2] = ctht_sav;
439  ctht_stat[ilin] = 1;
440  break;
441  }
442  }
443  }
444  // for any lines not transferred, process ctht
445  for( ilin = 0; ilin < ctht_nrec; ilin++ )
446  {
447  if( ctht_stat[ilin] == 0 )
448  {
449  lin_in_que = ilin + l1que.nq / 2 - ctht_nrec / 2;
450  scn_in_que = l1que.r[lin_in_que].iscan;
451  //if( comp_ctht_lin( l2rec->l1rec, ctht_lins.ct[ilin] ) != 0 )
452  if( comp_ctht_lin( &(l1que.r[lin_in_que]), ctht_lins.ct[ilin] ) != 0 )
453  {
454  printf( "%s - %d: Problem found running comp_ctht_lin\n",
455  __FILE__, __LINE__ );
456  exit(1);
457  }
458  }
459  }
460  return 0;
461  }
462 
463 int32_t comp_ctht_lin( l1str *l1rec, ctht_lin_str *ctht_lin )
464  /*
465  this will process the cloud top height for a line of satellite data
466 
467  Returns - int32_t - 0 for all good , only status for now
468 
469  Parameters: (in calling order)
470  Type Name I/O Description
471  ---- ---- --- -----------
472  l1str * l1rec input l1 record to use
473  ctht_lin_str *ctht_lin the cloud top height record to fill
474 
475  W. Robinson, SAIC, 22 Jun 2023, based on A. Sayer's IDL code
476  W. Robinson, SAIC, 3 May 2024, improve the merra 2 profile treatment
477  */
478  {
479  int32_t n_ctht_bnd = 8, n_prd, ny;
480  static int32_t ctht_init = 0;
481  static int32_t l1_bnd_ix[8];
482  static double *sy, *tmp_lut;
483  static ctht_lut_str full_lut_wtr, full_lut_ice;
484  static ctht_unc_str unc;
485  static oe_info_str oe_info;
486  static ctht_lut_str *lut;
487  static float *lcl_prof_t, *lcl_prof_h, *lcl_p_set;
488  static int32_t lcl_nlvl, lcl_lvl_ptr, fill_lvl;
489  float solz, senz, relaz, windsp, glint_ref;
490  int32_t ipix, ityp, npix, nbnd_ref, do_const;
491  int32_t ntyp = 2, foff, stype;
492  int32_t lut_dims[4];
493  float err1, err2, ev1, ev2;
494  double rhot[8], xa[3];
495  float *sxax[4], p_log_lo, p_log_hi;
496  int32_t snax[4];
497  int64_t n_block;
498  int32_t ibnd, ibnd2, icod, icth, ialb, icor, ib1, ib2, ilev;
499  // old:
500  //int32_t ctht_wav[] = { 755, 757, 760, 762, 765, 767, 770, 772 };
501  // new 14 May 2024
502  int32_t ctht_wav[] = { 754, 757, 759, 762, 764, 767, 769, 772 };
503  int32_t ctht_all_bnd = 0, close_ind;
504  int64_t ar_off[16], tmp_off, tmp_off2, off_cost;
505  float ar_wt[16], orig_wt[4];
506  float dp, found_hgt, lo_cost, spress_out;
507  int32_t bad_lt;
508  // for gsl vector, matrix work
509  double min_cost, yres[8], xres[8];
510  double fg_cost_t1, fg_cost_t2;
511  double int_pls, int_neg, int_ctr, coord, pctr;
512  int32_t min_cost_loc[3];
513  static double *fg_cost;
514  static gsl_matrix *gm_sa, *gm_sy, *gm_sy_inv, *gm_sa_inv;
515  static gsl_vector *gv_yres, *gv_sy_inv_yres, *gv_xres, *gv_sa_inv_xres;
516  // The final retrieved data to hang on to, yes, maybe we'll need both
517  // (or all) phases NOTE that ret_<name> will have all the product for
518  // phase 0 followed by phase 1
519  // phase 0 for ipix: ret_<name>[ ipix + npix * ityp ]
520  static float *ret_lcod;
521  // MERRA2 p levels
522  float merra_p_set[] = { 1000., 975., 950., 925., 900., 875., 850., 825.,
523  800., 775., 750., 725., 700., 650., 600., 550., 500., 450., 400.,
524  350., 300., 250., 200., 150., 100., 70., 50., 40., 30., 20., 10.,
525  7., 5., 4., 3., 2., 1., 0.7, 0.5, 0.4, 0.3, 0.1 };
526  int32_t nlvl_merra2 = 42, found_phase;
527 
528  // set some errors
529  unc.noise_err = 0.005;
530  unc.cal_err = 0.02;
531  float lut_dalb = 0.1; // Temporary till lut read in and albedo delta can be measured
532 
533  npix = l1rec->npix;
534  do_const = 2;
535  if( do_const != 2 )
536  {
537  comp_ctht_tst( l1rec, do_const, ctht_lin );
538  }
539  else
540  {
541  n_prd = 3;
542  ny = 8;
543  /*
544  * initialize the routine
545  */
546  if( ctht_init == 0 )
547  {
548  // get the location of ref_true (rhot) or bands to use for ctht
549  ctht_init = 1;
550  ctht_all_bnd = 1; // all bands there
551  for( ibnd = 0; ibnd < n_ctht_bnd; ibnd++ )
552  {
553  l1_bnd_ix[ibnd] = bindex_get( ctht_wav[ibnd] );
554  if( l1_bnd_ix[ibnd] < 0 )
555  {
556  ctht_all_bnd = 0;
557  printf( "%s, %d, I: CTHT wave: %d not found\n", __FILE__, __LINE__,
558  ctht_wav[ibnd] );
559  }
560  else
561  {
562  printf( "%s, %d, I: Got CTHT band wavelength %d as index %d\n",
563  __FILE__, __LINE__, ctht_wav[ibnd], l1_bnd_ix[ibnd] );
564  }
565  }
566  if( ctht_all_bnd == 0 )
567  {
568  printf( "%s, %d, E: Not all bands available for CTHT work\n",
569  __FILE__, __LINE__ );
570  exit(1);
571  }
572  // read in the cloud height tables uncertainty file parameters
573  if( ctht_tbl_init( n_ctht_bnd, l1rec, &unc, &full_lut_wtr, &full_lut_ice ) != 0 )
574  {
575  printf( "%s, %d, E: Failure reading unc or lut\n", __FILE__, __LINE__ );
576  exit(1);
577  }
578  // check the prods and bands (ny) vs expected and 1 table against
579  // the other
580  if( ( full_lut_wtr.lut_nwave != full_lut_ice.lut_nwave ) ||
581  ( full_lut_wtr.lut_nwave != ny ) )
582  {
583  printf( "%s, %d, E: expected # bands is not in the table\n",
584  __FILE__,__LINE__ );
585  exit(1);
586  }
587  if( ( full_lut_wtr.lut_nsolz != full_lut_ice.lut_nsolz ) ||
588  ( full_lut_wtr.lut_nsenz != full_lut_ice.lut_nsenz ) ||
589  ( full_lut_wtr.lut_nrelaz != full_lut_ice.lut_nrelaz ) ||
590  ( full_lut_wtr.lut_nlev != full_lut_ice.lut_nlev ) ||
591  ( full_lut_wtr.lut_npress != full_lut_ice.lut_npress ) ||
592  ( full_lut_wtr.lut_nalb != full_lut_ice.lut_nalb ) ||
593  ( full_lut_wtr.lut_ncod != full_lut_ice.lut_ncod ) ||
594  ( full_lut_wtr.lut_naod != full_lut_ice.lut_naod ) ||
595  ( full_lut_wtr.lut_ncth != full_lut_ice.lut_ncth ) )
596  {
597  printf( "%s, %d, E: LUTs (ice and water) have different array sizes\n",
598  __FILE__,__LINE__ );
599  exit(1);
600  }
601  // set some arrays used below
602  sy = (double *)malloc( n_ctht_bnd * n_ctht_bnd * sizeof(double) );
603 
604  if( ( fg_cost = (double *)
605  malloc( full_lut_wtr.lut_ncod * full_lut_wtr.lut_ncth *
606  full_lut_wtr.lut_nalb * sizeof(double) ) ) == NULL )
607  {
608  printf( "%s, %d, E: unable to allocate the cost array\n",
609  __FILE__, __LINE__ );
610  exit(1);
611  }
612  // set up the temp lut once
613  lut = &full_lut_wtr;
614  if( ( tmp_lut = (double *) malloc( lut->lut_nwave * lut->lut_ncod *
615  lut->lut_ncth * lut->lut_nalb * sizeof(double) ) ) == NULL )
616  {
617  printf( "%s - %d: E: Failure to allocate the pixel-specific LUT\n",
618  __FILE__, __LINE__ );
619  exit(1);
620  }
621  //
622  gm_sa = gsl_matrix_alloc( n_prd, n_prd );
623  gm_sy = gsl_matrix_alloc( ny, ny );
624  gv_yres = gsl_vector_alloc( ny );
625  gv_sy_inv_yres = gsl_vector_alloc( ny );
626  gv_xres = gsl_vector_alloc(n_prd);
627  gv_sa_inv_xres = gsl_vector_alloc(n_prd);
628 
629  // set up the outputs
630  // OK we could do an if( all mallocs == NULL err
631  if( ( ret_lcod = (float *) malloc( npix * ntyp * sizeof(float) ) ) == NULL )
632  {
633  printf(
634  "%s, %d, E: ctht was unable to allocate all the storage needed\n",
635  __FILE__, __LINE__ );
636  exit(1);
637  }
638  // set the local profiles that are specific to each pixel
639  if( ( ( lcl_prof_t = (float *) malloc( ( nlvl_merra2 + 1 ) * sizeof(float) ) ) == NULL )
640  || ( ( lcl_prof_h = (float *) malloc( ( nlvl_merra2 + 1 ) * sizeof(float) ) ) == NULL )
641  || ( ( lcl_p_set = (float *) malloc( ( nlvl_merra2 + 1 ) * sizeof(float) ) ) == NULL ) )
642  {
643  printf( "%s, %d, E: Unable to allocate local profile storage\n",
644  __FILE__, __LINE__ );
645  exit(1);
646  }
647  } // end general preparation
648 
649  npix = l1rec->npix;
650  nbnd_ref = l1rec->l1file->nbands;
651  /*
652  * loop through all the pixels
653  */
654  for( ipix = 0; ipix < npix; ipix++ )
655  {
656  /*
657  * do only the cloudy pixels, all others will be BAD_FLT
658  */
659  init_ctht_parms( ctht_lin, ipix );
660  // replace below with true eval of cloudy condition
661  if( l1rec->cloud[ipix] == 1 )
662 //printf( "%s,%d T: only testing pix 166\n", __FILE__, __LINE__ );
663 //if( ipix == 166 ) // WDR for test of odd point
664  {
665  solz = l1rec->solz[ipix];
666  senz = l1rec->senz[ipix];
667  relaz = l1rec->delphi[ipix];
668  windsp = l1rec->ws[ipix];
669 
670  // looks like relaz can be outside 0:180, so fix it
671  relaz = fmod( fabs( relaz ), 360. );
672  relaz = ( relaz <= 180. ) ? relaz : 360. - relaz;
673  // get toa reflectance for ctht bands
674  bad_lt = 0;
675  for( ibnd = 0; ibnd < n_ctht_bnd; ibnd++ )
676  {
677  foff = l1_bnd_ix[ibnd] + nbnd_ref * ipix;
678  if( !( isfinite( l1rec->Lt[foff] ) ) || l1rec->Lt[foff] == BAD_FLT )
679  bad_lt = 1;
680  rhot[ibnd] = l1rec->Lt[foff] * PI / ( cos(solz *
681  PI / 180. ) * l1rec->Fo[l1_bnd_ix[ibnd]] );
682  }
683  if( bad_lt ) break; // just skip this pixel
684  // get the pixel type: 0 water, 1 land, 2 snow/ice
685  stype = 0; // water
686  if( l1rec->land[ipix] ) stype = 1; // land
687  if( l1rec->ice[ipix] ) stype = 2; // snow, ice has precedence
688 
689 // Looks like the xa computation may be able to be done here and applied to both types
690 
691  // for this pixel, create the local profiles of p, t, h
692  // first level is the surface
693  lcl_prof_t[0] = l1rec->sfct[ipix];
694  lcl_prof_h[0] = l1rec->height[ipix];
695  lcl_p_set[0] = l1rec->sfcp[ipix];
696 
697  lcl_lvl_ptr = 0;
698  fill_lvl = 0;
699  // fill the profiles
700  for( int ilvl = 0; ilvl < nlvl_merra2; ilvl ++ )
701  {
702  // find the 1st level that is not fill and has p > the level p
703  // and the h < level h
704  if( lcl_lvl_ptr == 0 )
705  {
706  if( ( l1rec->anc_add->prof_height[ ipix * nlvl_merra2 + ilvl ]
707  != BAD_FLT ) &&
708  ( l1rec->anc_add->prof_temp[ ipix * nlvl_merra2 + ilvl ]
709  != BAD_FLT ) &&
710  ( merra_p_set[ilvl] < l1rec->sfcp[ipix] ) &&
711  ( l1rec->anc_add->prof_height[ ipix * nlvl_merra2 + ilvl ]
712  > l1rec->height[ipix] ) )
713  {
714  fill_lvl = 1;
715  lcl_lvl_ptr++;
716  }
717  }
718 
719  if( fill_lvl == 1 )
720  {
721  lcl_prof_t[lcl_lvl_ptr] = l1rec->anc_add->prof_temp[ ipix * nlvl_merra2 + ilvl ];
722  lcl_prof_h[lcl_lvl_ptr] = l1rec->anc_add->prof_height[ ipix * nlvl_merra2 + ilvl ];
723  lcl_p_set[lcl_lvl_ptr] = merra_p_set[ilvl];
724  lcl_lvl_ptr++;
725  }
726  }
727  lcl_nlvl = lcl_lvl_ptr;
728  // End local profile set-up
729  /*
730  * process for each cloud top particle type (phase)
731  */
732  for( ityp = 0; ityp < ntyp; ityp++ )
733  {
734  /*
735  * Set up the inputs for optimization
736  */
737  // Sy the uncerainty
738  for( ib1= 0; ib1 < n_ctht_bnd; ib1++ )
739  {
740  for( ib2 = 0; ib2 < n_ctht_bnd; ib2++ )
741  {
742  *( sy + ib1 + n_ctht_bnd * ib2 ) = 0.;
743  if( ib1 == ib2 )
744  {
745  // add random error
746  *( sy + ib1 + n_ctht_bnd * ib2 ) += unc.noise_err *
747  unc.noise_err * rhot[ib1] * rhot[ib2];
748  // add calib error
749  *( sy + ib1 + n_ctht_bnd * ib2 ) += unc.cal_err * unc.cal_err *
750  rhot[ib1] * rhot[ib2];
751  }
752  ev1 = unc.icpts[ ityp + ntyp * ( stype + unc.nsurf * ib1 ) ] +
753  unc.grads[ ityp + ntyp * ( stype + unc.nsurf * ib1 ) ] *
754  rhot[ib1];
755  ev2 = unc.fit_rms[ ityp + ntyp * ( stype + unc.nsurf * ib1 ) ];
756  err1 = ( ev1 > ev2 ) ? ev1 : ev2;
757 
758  ev1 = unc.icpts[ ityp + ntyp * ( stype + unc.nsurf * ib2 ) ] +
759  unc.grads[ ityp + ntyp * ( stype + unc.nsurf * ib2 ) ] *
760  rhot[ib2];
761  ev2 = unc.fit_rms[ ityp + ntyp * ( stype + unc.nsurf * ib2 ) ];
762  err2 = ( ev1 > ev2 ) ? ev1 : ev2;
763 
764  *( sy + ib1 + n_ctht_bnd * ib2 ) += err1 * err2 *
765  unc.err_corrs_category[ityp + ntyp *
766  ( stype + unc.nsurf * ( ib1 + n_ctht_bnd * ib2 ) ) ];
767  }
768  }
769  /*
770  * get the glint (Andy's version of Cox, Munk)
771  */
772  glint_ref = ctht_glint( solz, senz, relaz, windsp );
773  /*
774  * compute the xa - a-priori state vector [ hgt, p, albedo ]
775  */
776  double max_glint_for_table;
777  max_glint_for_table = lut->lut_alb[lut->lut_nalb-1];
778  xa[0] = 3.; xa[1] = 3.;
779  if( stype == 0 ) // water
780  {
781  xa[2] = ( 0.02 + glint_ref ) / lut_dalb;
782  }
783  else if( stype == 1 ) // land
784  {
785  xa[2] = l1rec->cld_dat->cth_alb_init[ipix] / lut_dalb;
786  }
787  else // snow/ice
788  {
789  xa[2] = 0.9 / lut_dalb;
790  }
791  if( xa[2] >= max_glint_for_table / lut_dalb )
792  {
793  //WDR avoid these as inst can view direct sun glint
794  //printf(
795  // "%s, %d: W: a-priori albedo exceeds table range, set to %f\n",
796  // __FILE__, __LINE__, max_glint_for_table );
797  xa[2] = max_glint_for_table / lut_dalb;
798  }
799  if( xa[2] < 0 )
800  {
801  //printf(
802  // "%s, %d: W: a-priori albedo below table range, set to 0\n",
803  // __FILE__, __LINE__ );
804  xa[2] = 0.;
805  }
806  /*
807  * make the covariance matrix sa
808  */
809  gsl_matrix_set_zero( gm_sa );
810  // Very weak prior on COD
811  gsl_matrix_set( gm_sa, 0, 0, pow( 100., 2 ) );
812  // Very weak prior on CTH
813  gsl_matrix_set( gm_sa, 1, 1, pow( 100., 2 ) );
814 
815  if( stype == 0 ) // water
816  {
817  // Stronger prior on albedo for water. This is +/-0.01. +
818  // 20% of the glint strength
819  gsl_matrix_set( gm_sa, 2, 2,
820  pow( ( 0.01 + 0.2 * glint_ref ) / lut_dalb, 2. ) );
821  }
822  else if( stype == 1 ) // land
823  {
824  gsl_matrix_set( gm_sa, 2, 2,
825  pow( ( l1rec->cld_dat->cth_alb_unc_init[ipix] / lut_dalb ), 2. ) );
826  // Prior comes from climatology
827  }
828  else
829  {
830  gsl_matrix_set( gm_sa, 2, 2, pow( 0.05 / lut_dalb, 2. ) );
831  // Overwrite if snow/ice cover > 50%,
832  // assume 0.05 uncertainty on it
833  }
834 
835  /*
836  * interpolate the lut to the current senz, solz, relaz, sfc press
837  */
838  lut = ( ityp == 0 ) ? &full_lut_wtr : &full_lut_ice;
839 
840  sxax[0] = lut->lut_press;
841  snax[0] = lut->lut_npress;
842  sxax[1] = lut->lut_solz;
843  snax[1] = lut->lut_nsolz;
844  sxax[2] = lut->lut_senz;
845  snax[2] = lut->lut_nsenz;
846  sxax[3] = lut->lut_relaz;
847  snax[3] = lut->lut_nrelaz;
848  n_block = lut->lut_ncod * lut->lut_ncth * lut->lut_nalb *
849  lut->lut_nwave;
850 
851  // For the idl match up work, I've been using l1rec->pr[ipix]
852  // as sfc pressure - seems to match what Andy has? but Andy
853  // did not want pr as it uses height and sea lvl p to get the sfc p
854  // and so wanted to use sfcp - ? a bit of an odd match. Just the
855  // same, the surface p is used in at least 3 parts below. I'll
856  // just declare spress_out here as one or the other
857  // spress_out = l1rec->pr[ipix]; // used for idl match work
858  spress_out = l1rec->sfcp[ipix];
859  float xpt[4] = { spress_out, solz, senz, relaz };
860 
861  // looks like the delphi can be outside range 0 - 180, so:
862  xpt[3] = fabs( fmod( xpt[3], 180. ) );
863  // This is a quick fix to get things running, the interpolation
864  // that got this bad p should get checked
865  xpt[0] = ( xpt[0] > 1050. ) ? 1050. : xpt[0];
866  xpt[0] = ( xpt[0] < 600. ) ? 600. : xpt[0];
867 
868  // set up. interpolate the LUT to the current point's 3 view
869  // angles and the surface p
870  // get the weights and offsets first
871 // **** this can be moved outside the 2 LUT types
872  if( int_4d( sxax, snax, n_block, xpt, ar_off, ar_wt, orig_wt ) != 0 )
873  {
874  //printf( "%s - %d: E: Failure in int_4d\n", __FILE__, __LINE__ );
875  //for( ibnd = 0; ibnd < 4; ibnd++ )
876  // printf( "x index: %d, value: %f, lo: %f, hi: %f\n", ibnd,
877  // xpt[ibnd], sxax[ibnd][0], sxax[ibnd][snax[ibnd]-1] );
878  break; // leave the prod arrays as bad for this point
879  }
880  // apply the weights to the LUT and reduce it from 9 to 4 dimensions
881  // [ band, cod, cth, alb ]
882  // apply the weights (and re-organize to have bands slowest)
883  for( ibnd = 0; ibnd < lut->lut_nwave; ibnd++ )
884  {
885  for( icod = 0; icod < lut->lut_ncod; icod++ )
886  {
887  for( icth = 0; icth < lut->lut_ncth; icth++ )
888  {
889  for( ialb = 0; ialb < lut->lut_nalb; ialb++ )
890  {
891  tmp_off = icod + lut->lut_ncod * ( icth + lut->lut_ncth *
892  ( ialb + lut->lut_nalb * ibnd ) );
893  tmp_off2 = ialb + lut->lut_nalb * ( icth + lut->lut_ncth *
894  ( icod + lut->lut_ncod * ibnd ) );
895  *( tmp_lut + tmp_off ) = 0;
896  for( icor = 0; icor < 16; icor++ )
897  {
898  *( tmp_lut + tmp_off ) += ar_wt[icor] *
899  *(lut->lut + ar_off[icor] + tmp_off2 );
900  }
901  }
902  }
903  }
904  }
905  /*
906  * get the first guess - the last piece in the puzzle
907  */
908  min_cost = 1.e11; // This should go down in checks below
909  for( icod = 0; icod < lut->lut_ncod; icod++ )
910  {
911  for( icth = 0; icth < lut->lut_ncth; icth++ )
912  {
913  // set the fg_cost to 1.e10, may shift to BAD later
914  for( ialb = 0; ialb < lut->lut_nalb; ialb++ )
915  {
916  // As orig: off_cost = ialb + lut->lut_nalb *
917  // ( icth + lut->lut_ncth * icod );
918  off_cost = icod + lut->lut_ncod *
919  ( icth + lut->lut_ncth * ialb );
920  *( fg_cost + off_cost ) = 1.e10;
921  }
922  ialb = (int32_t) ( xa[2] + .5 ); // a-priori albedo
923  // I orig thought off_cost = ialb + lut->lut_nalb *
924  // ( icth + lut->lut_ncth * icod );
925  // but below is correct
926  off_cost = icod + lut->lut_ncod * ( icth + lut->lut_ncth * ialb );
927  for( ibnd = 0; ibnd < lut->lut_nwave; ibnd++ )
928  {
929  tmp_off = icod + lut->lut_ncod * ( icth + lut->lut_ncth *
930  ( ialb + lut->lut_nalb * ibnd ) );
931  yres[ibnd] = *( tmp_lut + tmp_off ) - rhot[ibnd];
932  // ABOVE looks a lot like it could go in loop further up!!!
933  }
934  xres[0] = icod - xa[0];
935  xres[1] = icth - xa[1];
936  xres[2] = ialb - xa[2];
937 
938  // to replace Andy's code line:
939  // fg_cost[aa,bb,cc]=transpose(yres)#invert(Sy)#yres + $
940  // xres#invert(Sa)#xres
941  // make the matrix for sy and the vector for yres
942  for( ibnd = 0; ibnd < lut->lut_nwave; ibnd++ )
943  {
944  gsl_vector_set( gv_yres, ibnd, yres[ibnd] );
945  for( ibnd2= 0; ibnd2 < lut->lut_nwave; ibnd2++ )
946  {
947  gsl_matrix_set( gm_sy, ibnd, ibnd2,
948  *( sy + ibnd + 8 * ibnd2 ) );
949  }
950  }
951 
952  // term 1 of fg_cost: transpose(yres)#invert(Sy)#yres portion
953  gm_sy_inv = invert_a_matrix( gm_sy );
954  gsl_blas_dgemv( CblasNoTrans, 1., gm_sy_inv, gv_yres, 0., gv_sy_inv_yres );
955  gsl_matrix_free( gm_sy_inv );
956  gsl_blas_ddot( gv_sy_inv_yres, gv_yres, &fg_cost_t1 );
957  //
958  for( ibnd = 0; ibnd < 3; ibnd++ )
959  {
960  gsl_vector_set( gv_xres, ibnd, xres[ibnd] );
961  }
962 
963  // term 2 of fg_cost: xres#invert(Sa)#xres portion
964  gm_sa_inv = invert_a_matrix( gm_sa );
965  gsl_blas_dgemv( CblasNoTrans, 1., gm_sa_inv, gv_xres, 0.,
966  gv_sa_inv_xres );
967  gsl_matrix_free( gm_sa_inv );
968  gsl_blas_ddot( gv_sa_inv_xres, gv_xres, &fg_cost_t2 );
969  // sum the 2 values
970  // fg_cost[ialb,icth,icod] = fg_cost_t1 + fg_cost_t2;
971  *( fg_cost + off_cost ) = fg_cost_t1 + fg_cost_t2;
972 
973  // end Andy's code line:
974  /*
975  * find the locatuion and value of the fg_cost minimum
976  */
977  if( *( fg_cost + off_cost ) < min_cost )
978  {
979  min_cost = *( fg_cost + off_cost );
980  min_cost_loc[0] = icod;
981  min_cost_loc[1] = icth;
982  min_cost_loc[2] = ialb;
983  }
984  //
985  }
986  }
987  /*
988  * check the min cost so it is < 1.e10
989  */
990  if( min_cost >= 1.e10 )
991  {
992  printf( "%s, %d, I: the minimum cost is 1.e10. or higher\n",
993  __FILE__, __LINE__ );
994  }
995  /*
996  * perform the optimization, (ams_oe_inversion.pro part, maybe
997  * supplied by Laughlin
998  */
999  // oe_out=ams_oe_inversion(ref_true,Sy,xa,Sa,tmp_lut,fg=inds_min)
1000  // TO
1001  lut_dims[0] = lut->lut_ncod;
1002  lut_dims[1] = lut->lut_ncth;
1003  lut_dims[2] = lut->lut_nalb;
1004  lut_dims[3] = lut->lut_nwave;
1005  ams_oe_inversion( rhot, gm_sy, xa, gm_sa, tmp_lut,
1006  lut_dims, min_cost_loc, &oe_info );
1007  /*
1008  * post-optimization work
1009  */
1010  // get the retrieved parameters - xlate from index to value
1011  // WE'LL try a linear interp to emulate interpolate as Andy does
1012  // 26 calls to interpolate. Try on the 1st
1013  // of his
1014  // Do: ret_lcod[out_inds[0],out_inds[1],t] =
1015  // interpolate( lut_lcod, oe_out.x[0] )
1016 
1017  // Retrieved parameters
1018  // ret_lcod
1019  ctht_lin->cth_lcod_all[ ityp + ntyp * ipix ] =
1020  axis_interp( lut->lut_lcod, lut->lut_ncod, oe_info.x_prd_state[0] );
1021  // ret_cod
1022  ctht_lin->cth_cod_all[ ityp + ntyp * ipix ] =
1023  pow( 10., ctht_lin->cth_lcod_all[ ityp + ntyp * ipix ] );
1024  // for height, Use l1rec 'height'
1025  // ret_cth
1026  ctht_lin->cth_all[ ityp + ntyp * ipix ] =
1027  axis_interp( lut->lut_cth, lut->lut_ncth, oe_info.x_prd_state[1] )
1028  + l1rec->height[ipix] / 1000.;
1029  // ret_ctp
1030  ctht_lin->ctp_all[ ityp + ntyp * ipix ] =
1031  axis_interp( lut->lut_pp0, lut->lut_ncth, oe_info.x_prd_state[1] ) *
1032  spress_out;
1033  // ret_alb
1034  ctht_lin->cth_alb_all[ ityp + ntyp * ipix ] =
1035  axis_interp( lut->lut_alb, lut->lut_nalb, oe_info.x_prd_state[2] );
1036  // Diagnostics
1037  // ret_ss[npix,nlin,nphase] - cost
1038  ctht_lin->cost_ss[ ityp + ntyp * ipix ] = oe_info.cost;
1039  // ret_acost[npix,nlin,nphase]
1040  ctht_lin->acost[ ityp + ntyp * ipix ] = oe_info.acost;
1041  // ret_iter[npix,nlin,nphase]
1042  ctht_lin->nitr[ ityp + ntyp * ipix ] = oe_info.nitr;
1043 
1044  // Avg kernel - HOLD outputting this, need yet another 3rd dim size
1045  // oe_akdiag[npix,nlin,nphase,3]
1046  int nxx = 3;
1047  for( int xx = 0; xx < nxx; xx++ )
1048  ctht_lin->oe_akdiag[ xx + nxx * ( ityp + ntyp * ipix ) ] =
1049  gsl_matrix_get( oe_info.gm_ak, xx, xx );
1050 
1051  // Uncertainty estimates
1052  // ret_dlcod[npix,nlin,nphase]
1053  /* initial code:
1054  ul=10.^interpolate(lut_lcod,oe_out.x[0]+sqrt(oe_out.sx[0,0]) < $
1055  (lut_ncod-1)) - 10.^interpolate(lut_lcod,oe_out.x[0])
1056  ll=10.^interpolate(lut_lcod,oe_out.x[0]-sqrt(oe_out.sx[0,0]) > 0 ) - $
1057  10.^interpolate(lut_lcod,oe_out.x[0])
1058  ret_dcod[out_inds[0],out_inds[1],t]=max(abs([ul,ll]))
1059 
1060  ul=interpolate(lut_lcod,oe_out.x[0]+sqrt(oe_out.sx[0,0]) < $
1061  (lut_ncod-1)) - interpolate(lut_lcod,oe_out.x[0])
1062  ll=interpolate(lut_lcod,oe_out.x[0]-sqrt(oe_out.sx[0,0]) > 0 ) - $
1063  interpolate(lut_lcod,oe_out.x[0])
1064  ret_dlcod[out_inds[0],out_inds[1],t]=max(abs([ul,ll]))
1065  */
1066  // ret_dcod[npix,nlin,nphase]
1067  coord = oe_info.x_prd_state[0] +
1068  sqrt( gsl_matrix_get( oe_info.gm_sx, 0, 0 ) );
1069  if( lut->lut_ncod - 1 < coord ) coord = lut->lut_ncod - 1;
1070  int_ctr = axis_interp( lut->lut_lcod, lut->lut_ncod,
1071  oe_info.x_prd_state[0] );
1072  int_pls = axis_interp( lut->lut_lcod, lut->lut_ncod, coord );
1073 
1074  coord = oe_info.x_prd_state[0] -
1075  sqrt( gsl_matrix_get( oe_info.gm_sx, 0, 0 ) );
1076  if( coord < 0 ) coord = 0;
1077  int_neg = axis_interp( lut->lut_lcod, lut->lut_ncod, coord );
1078 
1079  ctht_lin->dlcod[ ityp + ntyp * ipix ] =
1080  fmax( fabs( int_pls - int_ctr ), fabs( int_neg -int_ctr ) );
1081  // ret_dcod[npix,nlin,nphase]
1082  pctr = pow( 10., int_ctr );
1083  int_pls = pow( 10., int_pls );
1084  int_neg = pow( 10., int_neg );
1085  ctht_lin->dcod[ ityp + ntyp * ipix ] =
1086  fmax( fabs( int_pls - pctr ), fabs( int_neg - pctr ) );
1087 
1088  // ret_dcth[npix,nlin,nphase]
1089  coord = oe_info.x_prd_state[1] +
1090  sqrt( gsl_matrix_get( oe_info.gm_sx, 1, 1 ) );
1091  if( lut->lut_ncth - 1 < coord ) coord = lut->lut_ncth - 1;
1092  int_ctr = axis_interp( lut->lut_cth, lut->lut_ncth,
1093  oe_info.x_prd_state[1] );
1094  int_pls = axis_interp( lut->lut_cth, lut->lut_ncth, coord );
1095 
1096  coord = oe_info.x_prd_state[1] -
1097  sqrt( gsl_matrix_get( oe_info.gm_sx, 1, 1 ) );
1098  if( coord < 0 ) coord = 0;
1099  int_neg = axis_interp( lut->lut_cth, lut->lut_ncth, coord );
1100 
1101  ctht_lin->dcth[ ityp + ntyp * ipix ] =
1102  fmax( fabs( int_pls - int_ctr ), fabs( int_neg - int_ctr ) );
1103 
1104  // Retrieved parameters (more which require MERRA2 profiles)
1105  // ret_cth_raw[nphase, npix,nlin]
1106  // the cth_raw will be the non-MERRA2 -derived height, in case
1107  // first, save the cth_all
1108  ctht_lin->cth_raw[ ityp + ntyp * ipix ] =
1109  ctht_lin->cth_all[ ityp + ntyp * ipix ];
1110  // find start of interval in MERRA2 levels containg the ctp
1111  close_ind = -1;
1112  for( ilev = 0; ilev < lcl_nlvl; ilev++ )
1113  {
1114  if( lcl_p_set[ilev] < ctht_lin->ctp_all[ ityp + ntyp * ipix ] )
1115  {
1116  close_ind = ( ilev == 0 ) ? 0 : ilev -1;
1117  break;
1118  }
1119  }
1120  if( ( close_ind == -1 ) || ( close_ind == lcl_nlvl-1 ) )
1121  {
1122  p_log_lo = log10( lcl_p_set[lcl_nlvl-2] );
1123  p_log_hi = log10( lcl_p_set[lcl_nlvl-1] );
1124  }
1125  else
1126  {
1127  p_log_lo = log10( lcl_p_set[close_ind] );
1128  p_log_hi = log10( lcl_p_set[ close_ind + 1 ] );
1129  }
1130  dp = ( log10( ctht_lin->ctp_all[ ityp + ntyp * ipix ] ) - p_log_lo )
1131  / ( p_log_hi - p_log_lo );
1132 
1133  // ret_cth (another calc using MERRA2)[npix,nlin,nphase] above
1134  // Note that we make the height into km as Andy has
1135  found_hgt = axis_interp( lcl_prof_h, lcl_nlvl, close_ind + dp );
1136 
1137  if( found_hgt < l1rec->height[ipix] ) found_hgt =
1138  l1rec->height[ipix] + 10.;
1139  ctht_lin->cth_all[ ityp + ntyp * ipix ] = found_hgt / 1000.;
1140 
1141  // ret_ctt (calc using MERRA2)[npix,nlin,nphase]
1142  ctht_lin->ctt_all[ ityp + ntyp * ipix ] = C_TO_K +
1143  axis_interp( lcl_prof_t, lcl_nlvl, close_ind + dp );
1144  // <a cth sanity check>
1145  if( ctht_lin->cth_all[ ityp + ntyp * ipix ] > 100. )
1146  {
1147  printf( "%s, %d, W: crazy CTH, profile issue, will use non MERRA\n",
1148  __FILE__, __LINE__ );
1149  ctht_lin->cth_all[ ityp + ntyp * ipix ] =
1150  ctht_lin->cth_raw[ ityp + ntyp * ipix ];
1151  }
1152  // <check the ctt, if > 350, use sfc t (10 m T )>
1153  if( ctht_lin->ctt_all[ ityp + ntyp * ipix ] >= 350. )
1154  ctht_lin->ctt_all[ ityp + ntyp * ipix ] =
1155  l1rec->sfct[ipix] + C_TO_K;
1156 
1157  // Uncertainty estimates, round 2
1158  // ret_dctp[npix,nlin,nphase]
1159  // coord: oe_out.x[1]+sqrt(oe_out.sx[1,1])
1160  // or lut_ncth-1, which ever lower
1161  coord = oe_info.x_prd_state[1] +
1162  sqrt( gsl_matrix_get( oe_info.gm_sx, 1, 1 ) );
1163  if( lut->lut_ncth - 1 < coord ) coord = lut->lut_ncth - 1;
1164  int_ctr = axis_interp( lut->lut_pp0, lut->lut_ncth,
1165  oe_info.x_prd_state[1] ) * spress_out;
1166  int_pls = axis_interp( lut->lut_pp0, lut->lut_ncth, coord ) *
1167  spress_out;
1168 
1169  coord = oe_info.x_prd_state[1] -
1170  sqrt( gsl_matrix_get( oe_info.gm_sx, 1, 1 ) );
1171  if( coord < 0 ) coord = 0;
1172  int_neg = axis_interp( lut->lut_pp0, lut->lut_ncth, coord ) *
1173  spress_out;
1174 
1175  ctht_lin->dctp[ ityp + ntyp * ipix ] =
1176  fmax( fabs( int_pls - int_ctr ), fabs( int_neg - int_ctr ) );
1177 
1178  // ret_dctt[npix,nlin,nphase] a bit easier, Likewise no super
1179  // easy way to make an equivalent for CTT uncertainty.
1180  // Since average lapse rate in troposphere is 6 K per km,
1181  // assume that holds, and scale based on CTH uncertainty.
1182  ctht_lin->dctt[ ityp + ntyp * ipix ] =
1183  6. * ctht_lin->dcth[ ityp + ntyp * ipix ];
1184 
1185  // ret_dalb[npix,nlin,nphase]
1186  coord = oe_info.x_prd_state[2] +
1187  sqrt( gsl_matrix_get( oe_info.gm_sx, 2, 2 ) );
1188  if( lut->lut_nalb - 1 < coord ) coord = lut->lut_nalb - 1;
1189  int_ctr = axis_interp( lut->lut_alb, lut->lut_nalb,
1190  oe_info.x_prd_state[2] );
1191  int_pls = axis_interp( lut->lut_alb, lut->lut_nalb, coord );
1192 
1193  coord = oe_info.x_prd_state[2] -
1194  sqrt( gsl_matrix_get( oe_info.gm_sx, 2, 2 ) );
1195  int_neg = axis_interp( lut->lut_alb, lut->lut_nalb, coord );
1196 
1197  ctht_lin->dalb[ ityp + ntyp * ipix ] =
1198  fmax( fabs( int_pls - int_ctr ), fabs( int_neg - int_ctr ) );
1199 
1200  } // end cloud top particle type
1201  // determine the 'best' phase - the phase will control what is
1202  // output
1203  lo_cost = 1.e20;
1204  found_phase = BAD_UBYTE;
1205  for( ityp = 0; ityp < ntyp; ityp++ )
1206  {
1207  if( ctht_lin->cost_ss[ ityp + ntyp * ipix ] == BAD_FLT )
1208  {
1209  found_phase = BAD_UBYTE;
1210  break;
1211  }
1212  else
1213  {
1214  if( ctht_lin->cost_ss[ ityp + ntyp * ipix ] < lo_cost )
1215  {
1216  found_phase = ityp;
1217  lo_cost = ctht_lin->cost_ss[ ityp + ntyp * ipix ];
1218  }
1219  }
1220  }
1221  ctht_lin->ct_phase[ipix] = found_phase;
1222  // set the main 3 products that go into the CMP code
1223  if( found_phase != BAD_UBYTE )
1224  {
1225  ctht_lin->cth[ipix] = ctht_lin->cth_all[ found_phase + ntyp * ipix ];
1226  ctht_lin->ctp[ipix] = ctht_lin->ctp_all[ found_phase + ntyp * ipix ];
1227  ctht_lin->ctt[ipix] = ctht_lin->ctt_all[ found_phase + ntyp * ipix ];
1228  }
1229  else
1230  {
1231  ctht_lin->cth[ipix] = BAD_FLT;
1232  ctht_lin->ctp[ipix] = BAD_FLT;
1233  ctht_lin->ctt[ipix] = BAD_FLT;
1234  }
1235  } // end cloudy pixel work, may need to set some BADs here too
1236  } // end pixel loop
1237 // WDR temp, stop after 1st line
1238 // printf( "WDR - get_ctht early end\n" );
1239 // exit(8);
1240 
1241  // set scan and # pix
1242  ctht_lin->iscan = l1rec->iscan;
1243  ctht_lin->npix = l1rec->npix;
1244  }
1245  return 0;
1246  }
1247 int32_t init_ctht_parms( ctht_lin_str *ctht_lin, int32_t ipix )
1248  /*
1249  * just initialize the cloud top height values
1250  *
1251  * Return 0 - all well
1252  * ctht_lin_str * ctht_lin structure of all ctht info
1253  * int32_t ipix pixell to set values at
1254  */
1255  {
1256  int32_t ntyp = 2;
1257 
1258  // basic values needed by chimaera code
1259  ctht_lin->cth[ipix] = BAD_FLT;
1260  ctht_lin->ctp[ipix] = BAD_FLT;
1261  ctht_lin->ctt[ipix] = BAD_FLT;
1262  // other cloud height values
1263  ctht_lin->ct_phase[ipix] = BAD_UBYTE;
1264 
1265  // 'all' from all phases
1266  for( int i = 0; i < ntyp; i++ )
1267  {
1268  ctht_lin->cth_all[ i + ntyp * ipix ] = BAD_FLT;
1269  ctht_lin->ctp_all[ i + ntyp * ipix ] = BAD_FLT;
1270  ctht_lin->ctt_all[ i + ntyp * ipix ] = BAD_FLT;
1271  ctht_lin->cth_cod_all[ i + ntyp * ipix ] = BAD_FLT;
1272  ctht_lin->cth_lcod_all[ i + ntyp * ipix ] = BAD_FLT;
1273  ctht_lin->cth_raw[ i + ntyp * ipix ] = BAD_FLT;
1274  ctht_lin->cth_alb_all[ i + ntyp * ipix ] = BAD_FLT;
1275  ctht_lin->cost_ss[ i + ntyp * ipix ] = BAD_FLT;
1276  ctht_lin->acost[ i + ntyp * ipix ] = BAD_FLT;
1277  ctht_lin->nitr[ i + ntyp * ipix ] = BAD_INT;
1278  for( int inx = 0; inx < 3; inx++ )
1279  ctht_lin->oe_akdiag[ inx * ( i + ntyp * ipix ) ] = BAD_FLT;
1280  ctht_lin->dlcod[ i + ntyp * ipix ] = BAD_FLT;
1281  ctht_lin->dcod[ i + ntyp * ipix ] = BAD_FLT;
1282  ctht_lin->dcth[ i + ntyp * ipix ] = BAD_FLT;
1283  ctht_lin->dctp[ i + ntyp * ipix ] = BAD_FLT;
1284  ctht_lin->dctt[ i + ntyp * ipix ] = BAD_FLT;
1285  ctht_lin->dalb[ i + ntyp * ipix ] = BAD_FLT;
1286  }
1287  return 0;
1288  }
1289 int32_t comp_ctht_tst(l1str *l1rec, int32_t proc_flg, ctht_lin_str *ctht_lin )
1290  /*
1291  comp_ctht_tst will supply test ctht data if called for a line of
1292  satellite data
1293 
1294  Returns - int32_t - 0 for all good , only status for now
1295 
1296  Parameters: (in calling order)
1297  Type Name I/O Description
1298  ---- ---- --- -----------
1299  l1str * l1rec input l1 record to use
1300  int32_t proc_flg either 0 to return constants, or 1 to return
1301  values based on scan mod 3
1302  ctht_lin_str *ctht_lin the cloud top height record to fill
1303 
1304  W. Robinson, SAIC, 22 Jun 2023
1305  */
1306  {
1307  int32_t ipix, npix;
1308  float cth, ctp, ctt;
1309 
1310  npix = l1rec->npix;
1311  if( proc_flg == 0 )
1312  {
1313  cth = 3.;
1314  ctp = 480.;
1315  ctt = 250.;
1316  }
1317  else if( proc_flg== 1 )
1318  {
1319  switch ( l1rec->iscan % 3 )
1320  {
1321  case 0:
1322  cth = 2.;
1323  ctp = 600.;
1324  ctt = 270.;
1325  break;
1326  case 1:
1327  cth = 3.;
1328  ctp = 480.;
1329  ctt = 250.;
1330  break;
1331  case 2:
1332  cth = 4.;
1333  ctp = 400.;
1334  ctt = 230.;
1335  break;
1336  }
1337  }
1338  for( ipix = 0; ipix < npix; ipix++ )
1339  {
1340  ctht_lin->cth[ipix] = cth;
1341  ctht_lin->ctp[ipix] = ctp;
1342  ctht_lin->ctt[ipix] = ctt;
1343  }
1344  ctht_lin->iscan = l1rec->iscan;
1345  ctht_lin->npix = l1rec->npix;
1346  return 0;
1347  }
1348 
1349 int32_t ctht_tbl_init( int32_t nband, l1str *l1rec, ctht_unc_str *unc,
1350  ctht_lut_str *full_lut_wtr, ctht_lut_str *full_lut_ice )
1351  /*
1352  * ctht_tbl_init - initialize the cloud height tables
1353  *
1354  * returns int32_t status: 0 good else bad
1355  * int32_t nband input expected # bands
1356  * l1str * l1rec, the L1 record
1357  * ctht_unc_str *unc output filled uncertainty information
1358  * ctht_lut_str *full_lut_wtr output filled LUT information for water
1359  * ctht_lut_str *full_lut_ice output filled LUT information for ice
1360  */
1361  {
1362  int ncid, dim_loc, varid;
1363  char *env_arr;
1364  char unc_file[FILENAME_MAX], lut_wtr_file[FILENAME_MAX], lut_ice_file[FILENAME_MAX];
1365  char tit_unc[] = { "Forward model uncertainty coefficients for cloud top pressure retrieval from oci" };
1366  char tit_lut[] = { "LUT for cloud top pressure retrieval from PACE OCI" };
1367  ctht_lut_str *lut_ptr;
1368  char *lut_fil;
1369  int32_t ilut, ilev, ihgt, close_ind;
1370  float pstep, frac_ind;
1371 
1372  /* set the path to the tables and table names */
1373  if ((env_arr = getenv("OCDATAROOT")) == NULL) {
1374  printf("-E- %s, %d: Error looking up environmental variable OCDATAROOT\n", __FILE__, __LINE__ );
1375  exit(1);
1376  }
1377  strcpy( unc_file, env_arr );
1378  strcat( unc_file, "/cloud/" );
1379  strcat( unc_file, sensorId2SensorDir(l1rec->l1file->sensorID));
1380  strcat( unc_file, "/" );
1381  strcpy( lut_wtr_file, unc_file );
1382  strcpy( lut_ice_file, unc_file );
1383  // make individual file names here
1384  strcat( unc_file, "oci_fm_err_fits_lut_v17_sim_v5.nc" );
1385  strcat( lut_wtr_file, "oci_o2_cth_lut_v17_water.nc" );
1386  strcat( lut_ice_file, "oci_o2_cth_lut_v17_column_8elements.nc" );
1387 
1388  /* open the uncertainty file */
1389  printf( "%s, %d: I: Reading cloud top uncertainty file: %s\n",
1390  __FILE__, __LINE__, unc_file );
1391  DPTB( nc_open( unc_file, 0, &ncid ) );
1392  /* check the title for a match */
1393  /* must do the below for a string attr, not text */
1394  size_t attlen = 0, dim_len;
1395  DPTB( nc_inq_attlen(ncid, NC_GLOBAL, "description", &attlen) );
1396  char **string_attr = (char**)malloc(attlen * sizeof(char*));
1397  memset(string_attr, 0, attlen * sizeof(char*));
1398  DPTB( nc_get_att_string(ncid, NC_GLOBAL, "description", string_attr) );
1399 
1400  if( strncmp( string_attr[0], tit_unc, strlen(tit_unc) ) != 0 )
1401  {
1402  printf( "%s, %d: E - Description of ctht uncertainty file is bad\n",
1403  __FILE__, __LINE__ );
1404  exit(1);
1405  }
1406  /* read sizes */
1407  DPTB( nc_inq_dimid( ncid, "nphase", &dim_loc ) );
1408  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1409  unc->nphase = dim_len;
1410 
1411  DPTB( nc_inq_dimid( ncid, "nsurf", &dim_loc ) );
1412  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1413  unc->nsurf = dim_len;
1414 
1415  DPTB( nc_inq_dimid( ncid, "nbands", &dim_loc ) );
1416  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1417  unc->nbands = dim_len;
1418  /* read axis coordinates */
1419  /* in this case the wavelengths 'better' match ctht_wav, we have 2 phases
1420  water and ice and we have 3 surfaces, 0 water, 1 land, 2 snowice */
1421  /* read the uncertainty arrays */
1422  DPTB( nc_inq_varid( ncid, "icpts", &varid ) );
1423  unc->icpts = (float *)malloc( unc->nphase * unc->nsurf * unc->nbands *
1424  sizeof(float) );
1425  DPTB( nc_get_var_float( ncid, varid, unc->icpts ) );
1426 
1427  DPTB( nc_inq_varid( ncid, "grads", &varid ) );
1428  unc->grads = (float *)malloc( unc->nphase * unc->nsurf * unc->nbands *
1429  sizeof(float) );
1430  DPTB( nc_get_var_float( ncid, varid, unc->grads ) );
1431 
1432  DPTB( nc_inq_varid( ncid, "fit_rms", &varid ) );
1433  unc->fit_rms = (float *)malloc( unc->nphase * unc->nsurf * unc->nbands *
1434  sizeof(float) );
1435  DPTB( nc_get_var_float( ncid, varid, unc->fit_rms ) );
1436 
1437  DPTB( nc_inq_varid( ncid, "err_corrs_category", &varid ) );
1438  unc->err_corrs_category = (float *)malloc( unc->nphase * unc->nsurf *
1439  unc->nbands * unc->nbands * sizeof(float) );
1440  DPTB( nc_get_var_float( ncid, varid, unc->err_corrs_category ) );
1441  DPTB( nc_close( ncid ) );
1442  /*
1443  * open the LUT file for water and ice
1444  */
1445  for( ilut = 0; ilut < 2; ilut++ )
1446  {
1447  lut_ptr = ( ilut == 0 ) ? full_lut_wtr : full_lut_ice;
1448  lut_fil = ( ilut == 0 ) ? lut_wtr_file : lut_ice_file;
1449 
1450  printf( "%s, %d: I: Reading cloud top phase LUT file: %s\n",
1451  __FILE__, __LINE__, lut_fil );
1452  DPTB( nc_open( lut_fil, 0, &ncid ) );
1453  /* check the title for a match */
1454  size_t attlen = 0, dim_len;
1455  DPTB( nc_inq_attlen(ncid, NC_GLOBAL, "description", &attlen) );
1456  char **string_attr = (char**)malloc(attlen * sizeof(char*));
1457  memset(string_attr, 0, attlen * sizeof(char*));
1458  DPTB( nc_get_att_string(ncid, NC_GLOBAL, "description", string_attr) );
1459 
1460  if( strncmp( string_attr[0], tit_lut, strlen(tit_lut) ) != 0 )
1461  {
1462  printf( "%s, %d: E - Description of ctht lut file is bad\n",
1463  __FILE__, __LINE__ );
1464  printf( "For file: %s\n", lut_fil );
1465  exit(1);
1466  }
1467  /* read sizes */
1468  DPTB( nc_inq_dimid( ncid, "nsza", &dim_loc ) );
1469  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1470  lut_ptr->lut_nsolz = dim_len;
1471 
1472  DPTB( nc_inq_dimid( ncid, "nvza", &dim_loc ) );
1473  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1474  lut_ptr->lut_nsenz = dim_len;
1475 
1476  DPTB( nc_inq_dimid( ncid, "nazi", &dim_loc ) );
1477  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1478  lut_ptr->lut_nrelaz = dim_len;
1479 
1480  DPTB( nc_inq_dimid( ncid, "nlevs", &dim_loc ) );
1481  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1482  lut_ptr->lut_nlev = dim_len;
1483 
1484  DPTB( nc_inq_dimid( ncid, "npress", &dim_loc ) );
1485  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1486  lut_ptr->lut_npress = dim_len;
1487 
1488  DPTB( nc_inq_dimid( ncid, "nalb", &dim_loc ) );
1489  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1490  lut_ptr->lut_nalb = dim_len;
1491 
1492  DPTB( nc_inq_dimid( ncid, "ncod", &dim_loc ) );
1493  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1494  lut_ptr->lut_ncod = dim_len;
1495 
1496  DPTB( nc_inq_dimid( ncid, "naod", &dim_loc ) );
1497  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1498  lut_ptr->lut_naod = dim_len;
1499 
1500  DPTB( nc_inq_dimid( ncid, "nbands", &dim_loc ) );
1501  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1502  lut_ptr->lut_nwave = dim_len;
1503 
1504  DPTB( nc_inq_dimid( ncid, "ncth", &dim_loc ) );
1505  DPTB( nc_inq_dimlen( ncid, dim_loc, &dim_len ) );
1506  lut_ptr->lut_ncth = dim_len;
1507 
1508  //printf( "%s, %d: TEMP: nsolz: %d, nsenz: %d, nrelaz: %d\n", __FILE__, __LINE__, lut_ptr->lut_nsolz, lut_ptr->lut_nsenz, lut_ptr->lut_nrelaz );
1509  //printf( "nlev: %d, npress: %d, nalb: %d, ncod: %d\n", lut_ptr->lut_nlev, lut_ptr->lut_npress, lut_ptr->lut_nalb, lut_ptr->lut_ncod );
1510  //printf( "naod: %d, nbands: %d, ncth: %d\n", lut_ptr->lut_naod, lut_ptr->lut_nwave, lut_ptr->lut_ncth );
1511 
1512  if( lut_ptr->lut_nwave != nband )
1513  {
1514  printf( "%s, %d: E: CTHT LUT # bands != expected #: %d\n", __FILE__, __LINE__, nband );
1515  exit(1);
1516  }
1517  /* read axis coordinates */
1518  lut_ptr->lut_solz = (float *) malloc( lut_ptr->lut_nsolz * sizeof( float) );
1519  lut_ptr->lut_senz = (float *) malloc( lut_ptr->lut_nsenz * sizeof( float) );
1520  lut_ptr->lut_relaz = (float *) malloc( lut_ptr->lut_nrelaz * sizeof( float) );
1521  if( ( lut_ptr->lut_solz == NULL ) || ( lut_ptr->lut_senz == NULL ) || ( lut_ptr->lut_relaz == NULL ) )
1522  {
1523  printf( "%s, %d: E: unable to allocate lut_solz, lut_senz, or lut_relaz\n", __FILE__, __LINE__ );
1524  exit(1);
1525  }
1526  DPTB( nc_inq_varid( ncid, "sza", &varid ) );
1527  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_solz ) );
1528 
1529  DPTB( nc_inq_varid( ncid, "vza", &varid ) );
1530  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_senz) );
1531 
1532  DPTB( nc_inq_varid( ncid, "azi", &varid ) );
1533  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_relaz ) );
1534 
1535  lut_ptr->lut_press = (float *) malloc( lut_ptr->lut_npress * sizeof( float) );
1536  lut_ptr->lut_alb = (float *) malloc( lut_ptr->lut_nalb * sizeof( float) );
1537  lut_ptr->lut_cod = (float *) malloc( lut_ptr->lut_ncod * sizeof( float) );
1538  lut_ptr->lut_lcod = (float *) malloc( lut_ptr->lut_ncod * sizeof( float) );
1539  if( ( lut_ptr->lut_press == NULL ) || ( lut_ptr->lut_alb == NULL ) || ( lut_ptr->lut_cod == NULL ) )
1540  {
1541  printf( "%s, %d: E: unable to allocate lut_press, lut_alb, or lut_cod\n", __FILE__, __LINE__ );
1542  exit(1);
1543  }
1544  DPTB( nc_inq_varid( ncid, "spress", &varid ) );
1545  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_press ) );
1546 
1547  DPTB( nc_inq_varid( ncid, "cod", &varid ) );
1548  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_cod ) );
1549 
1550  DPTB( nc_inq_varid( ncid, "lcod", &varid ) );
1551  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_lcod ) );
1552 
1553  DPTB( nc_inq_varid( ncid, "salb", &varid ) );
1554  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_alb ) );
1555 
1556  lut_ptr->lut_aod = (float *) malloc( lut_ptr->lut_naod * sizeof( float) );
1557  lut_ptr->lut_wave = (float *) malloc( lut_ptr->lut_nwave * sizeof( float) );
1558  lut_ptr->lut_cth = (float *) malloc( lut_ptr->lut_ncth * sizeof( float) );
1559  if( ( lut_ptr->lut_aod == NULL ) || ( lut_ptr->lut_wave == NULL ) || ( lut_ptr->lut_cth == NULL ) )
1560  {
1561  printf( "%s, %d: E: unable to allocate lut_aod, lut_wave, or lut_cth\n", __FILE__, __LINE__ );
1562  exit(1);
1563  }
1564  DPTB( nc_inq_varid( ncid, "aod", &varid ) );
1565  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_aod ) );
1566 
1567  DPTB( nc_inq_varid( ncid, "band_wavelengths", &varid ) );
1568  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_wave ) );
1569 
1570  DPTB( nc_inq_varid( ncid, "cth", &varid ) );
1571  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut_cth ) );
1572 
1573  /* there is a state_inds array, but avoid this for now */
1574  /* read the LUT array and height rel press */
1575  lut_ptr->z_and_pp0 = (float *) malloc( lut_ptr->lut_nlev * 2 * sizeof( float) );
1576  lut_ptr->lut_pp0 = (float *) malloc( lut_ptr->lut_ncth * sizeof( float ) );
1577  lut_ptr->lut = (float *) malloc( lut_ptr->lut_nrelaz * lut_ptr->lut_nsenz
1578  * lut_ptr->lut_nsolz * lut_ptr->lut_npress * lut_ptr->lut_nwave *
1579  lut_ptr->lut_ncod * lut_ptr->lut_ncth * lut_ptr->lut_nalb *
1580  sizeof( float) );
1581  if( ( lut_ptr->z_and_pp0 == NULL ) || ( lut_ptr->lut == NULL ) ||
1582  ( lut_ptr->lut_pp0 == NULL ) )
1583  {
1584  printf( "%s, %d: E: unable to allocate lut_z_and_pp0, lut_pp0, or lut\n",
1585  __FILE__, __LINE__ );
1586  exit(1);
1587  }
1588  DPTB( nc_inq_varid( ncid, "z_and_pp0", &varid ) );
1589  DPTB( nc_get_var_float( ncid, varid, lut_ptr->z_and_pp0 ) );
1590 
1591  DPTB( nc_inq_varid( ncid, "lut", &varid ) );
1592  DPTB( nc_get_var_float( ncid, varid, lut_ptr->lut ) );
1593  // that should do it
1594  DPTB( nc_close( ncid ) );
1595  /*
1596  * The lut_pp0 is created from the z_and_pp0 - do that here
1597  */
1598  float *and_pp0 = (float *) malloc( lut_ptr->lut_nlev * sizeof(float) );
1599  for( ihgt = 0; ihgt < lut_ptr->lut_nlev; ihgt++ )
1600  and_pp0[ihgt] = lut_ptr->z_and_pp0[ 1 + 2 * ihgt ];
1601  for( ilev = 0; ilev < lut_ptr->lut_ncth; ilev++ )
1602  {
1603  // get the pressure for all heights in the height axis
1604  lut_ptr->lut_pp0[ilev] = BAD_FLT;
1605  close_ind = -1;
1606  for( ihgt = 0; ihgt < lut_ptr->lut_nlev; ihgt++ )
1607  {
1608  if( lut_ptr->z_and_pp0[ ihgt * 2 ] < lut_ptr->lut_cth[ilev] )
1609  {
1610  close_ind = ( ihgt == 0 ) ? 0 : ihgt - 1;
1611  break;
1612  }
1613  }
1614  if( close_ind < lut_ptr->lut_nlev - 1 )
1615  pstep = lut_ptr->z_and_pp0[ ( close_ind + 1 ) * 2 ] -
1616  lut_ptr->z_and_pp0[ close_ind * 2 ];
1617  else
1618  pstep = lut_ptr->z_and_pp0[ close_ind * 2 ] -
1619  lut_ptr->z_and_pp0[ ( close_ind - 1 ) * 2 ];
1620  frac_ind = close_ind + ( lut_ptr->lut_cth[ilev] -
1621  lut_ptr->z_and_pp0[ close_ind * 2 ] ) / pstep;
1622  lut_ptr->lut_pp0[ilev] = axis_interp( and_pp0, lut_ptr->lut_nlev,
1623  frac_ind );
1624  }
1625  free(and_pp0);
1626  }
1627  /*
1628  * Both tables read
1629  */
1630  return 0;
1631  }
1632 
1633 #define DTOR PI / 180.
1634 
1635 float ctht_glint( float solz, float senz, float relaz, float windsp )
1636  /*
1637  * ctht_glint will vget the glint for a pixel and wind speed, adapted from
1638  * calc_sunglint.pro f A. Sayer
1639  *
1640  * return float glint value as reflectance
1641  * float solz solar zenith, in deg
1642  * float senz sensor zenith, degrees
1643  * float relaz relative azimuth
1644  * float windsp wind speed, m/s
1645  *
1646  * follows the Cox and Munk 1954 paper
1647  *
1648  * 24 Jan 2024, WDR make more variables as double and avoid NaN in
1649  * R calculation
1650  *
1651  */
1652  {
1653  double cos_two_omega, cos_beta;
1654  double a1, b1, c1, d1, rgl, wprime, ref_coef;
1655  float kai, nr;
1656  double zxprime, zyprime, sigx, sigy, w;
1657  double zeta, eta, p_gram, zx, zy;
1658  double dsenz, dsolz, drelaz;
1659 
1660  int32_t do_gordon = 1;
1661 
1662  dsenz = (double) senz * DTOR;
1663  dsolz = (double) solz * DTOR;
1664  drelaz = ( do_gordon ) ? ( 180. - (double)relaz ) * DTOR :
1665  (double)relaz * DTOR;
1666 
1667  nr = 1.341;
1668 
1669  // Specular reflectance from ruffled sea
1670  kai = 0.;
1671  rgl = 0.0;
1672  kai = kai * DTOR;
1673 
1674  // Defines surface slopes
1675  zx = ( -1.0 * sin(dsenz) * sin(drelaz) ) / ( cos(dsolz) + cos(dsenz) );
1676  zy = ( sin(dsolz) + sin(dsenz) * cos( drelaz ) ) /
1677  ( cos(dsolz) + cos(dsenz) );
1678 
1679  // Use these lines to make independent of wind direction
1680  zxprime = zx;
1681  zyprime = zy;
1682 
1683  // Coefficients from Cox and Munk, 1954 (Statistics of...).
1684  // Wind isotropic slope as function of find speed.
1685  sigx = sqrt( 0.003 + 0.00192 * windsp );
1686  sigy = sqrt( 0.00316 * windsp );
1687 
1688  zeta = zxprime / sigx;
1689  eta = zyprime / sigy;
1690 
1691  p_gram = ( 1. / ( 2. * PI * sigx * sigy ) ) *
1692  exp( -0.5 * ( pow( zeta, 2. ) + pow( eta, 2. ) ) );
1693 
1694  // Cox and Munk (1954) geometry (Measurements of...)
1695  // 2 omega = angle between incident ray (sun) and instrument
1696  // with the surface as the sloping sea
1697 
1698  cos_two_omega = cos(dsenz) * cos(dsolz) +
1699  sin(dsenz) * sin( dsolz ) * cos(drelaz);
1700  cos_beta = ( cos(dsolz) + cos(dsenz) ) / ( sqrt( 2 + 2 * cos_two_omega ) );
1701  w = ( cos_two_omega >= 1. ) ? 0. : 0.5 * acos( cos_two_omega );
1702 
1703  // Fresnel component
1704  wprime = asin( 1.00029 * sin(w) / nr );
1705  a1 = sin( w - wprime );
1706  b1 = sin( w + wprime );
1707  c1 = tan( w - wprime );
1708  d1 = tan( w + wprime );
1709  // make glint refl 0 if b1 or d1 is 0
1710  if( ( b1 == 0 ) || ( d1 == 0 ) )
1711  rgl = 0.;
1712  else {
1713  ref_coef = 0.5 * ( ( a1 * a1 ) / ( b1 * b1 ) + ( c1 * c1 ) / ( d1 * d1 ) );
1714  // compute the glint reflectance
1715  rgl = PI * p_gram * ref_coef / ( 4. * cos( dsolz ) * cos( dsenz ) *
1716  ( pow( cos_beta, 4. ) ) );
1717  }
1718  return rgl;
1719  }
1720 // the matrix invert:
1721 
1722 gsl_matrix *invert_a_matrix(gsl_matrix *matrix)
1723  /**************************************************
1724  * invert_a_matrix was copied as a way to do an inverse of a gsl matrix
1725  * The only thing was that it affects the input matris. So, we copy it
1726  * and work on that
1727  *************************************************/
1728 {
1729  size_t size = matrix->size1;
1730  gsl_permutation *p = gsl_permutation_alloc(size);
1731  gsl_matrix *gm_work;
1732  int s;
1733 
1734  // copy the input matris to work on
1735  gm_work = gsl_matrix_alloc( size, size );
1736  gsl_matrix_memcpy( gm_work, matrix );
1737 
1738  // Compute the LU decomposition of this matrix
1739  gsl_linalg_LU_decomp(gm_work, p, &s);
1740 
1741  // Compute the inverse of the LU decomposition
1742  gsl_matrix *inv = gsl_matrix_alloc(size, size);
1743  gsl_linalg_LU_invert(gm_work, p, inv);
1744 
1745  gsl_permutation_free(p);
1746  gsl_matrix_free( gm_work );
1747 
1748  return inv;
1749 }
1750 
1751 // a diagnostic routines to allow printing matricise and vectors
1752 // (the gsl_matrix_fprintf is bit vague)
1753 // the gsl_vector_fprintf has no such confusion
1754 void
1755 print_mat_contents(gsl_matrix *matrix, int nrow, int ncol )
1756 {
1757  double element;
1758  int i, j;
1759  printf( "%d rows, %d columns\n", nrow, ncol );
1760  for (i = 0; i < nrow; ++i) {
1761  for (j = 0; j < ncol; ++j) {
1762  element = gsl_matrix_get(matrix, i, j);
1763  printf("%f ", element);
1764  }
1765  printf("\n");
1766  }
1767 }
int_4d
int32_t int_4d(float *sxax[4], int32_t snax[4], int64_t n_block, float x[4], int64_t[16], float[16], float[4])
Definition: int_4d.c:12
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
ctht_tbl_init
int32_t ctht_tbl_init(int32_t nband, l1str *l1rec, ctht_unc_str *unc, ctht_lut_str *full_lut_wtr, ctht_lut_str *full_lut_ice)
Definition: get_ctht.c:1349
CAT_cth_cost
#define CAT_cth_cost
Definition: l2prod.h:522
j
int j
Definition: decode_rs.h:73
CAT_cth_cost_all
#define CAT_cth_cost_all
Definition: l2prod.h:504
CAT_cth_dctp_all
#define CAT_cth_dctp_all
Definition: l2prod.h:515
comp_ctht_lin
int32_t comp_ctht_lin(l1str *l1rec, ctht_lin_str *ctht_lin)
Definition: get_ctht.c:463
get_ctht_lin
float * get_ctht_lin(l2str *l2rec, int prodnum)
Definition: get_ctht.c:14
NULL
#define NULL
Definition: decode_rs.h:63
axis_interp
float axis_interp(float *axis_arr_vals, int32_t n_ax, float ind_frac)
Definition: axis_interp.c:19
CAT_Cld_h
#define CAT_Cld_h
Definition: l2prod.h:484
coord
subroutine coord(X, GM, Y)
Definition: coord.f:2
ams_oe_inversion
int32_t ams_oe_inversion(double *rhot, gsl_matrix *gm_sy, double *xa, gsl_matrix *gm_sa, double *lut, int32_t *lut_dims, int32_t *min_cost_loc, oe_info_str *oe_info)
Definition: ams_oe_inversion.c:9
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
print_mat_contents
void print_mat_contents(gsl_matrix *matrix, int nrow, int ncol)
Definition: get_ctht.c:1755
CAT_cth_dlcod_all
#define CAT_cth_dlcod_all
Definition: l2prod.h:510
CAT_cth_dctt
#define CAT_cth_dctt
Definition: l2prod.h:536
CAT_cth_lcod_all
#define CAT_cth_lcod_all
Definition: l2prod.h:500
BAD_UBYTE
#define BAD_UBYTE
Definition: genutils.h:26
matrix
float ** matrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:60
CAT_Cld_t
#define CAT_Cld_t
Definition: l2prod.h:425
get_ctht.h
CAT_cth_ctp_all
#define CAT_cth_ctp_all
Definition: l2prod.h:502
DTOR
#define DTOR
Definition: get_ctht.c:1633
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
CAT_cth_dctp
#define CAT_cth_dctp
Definition: l2prod.h:535
input
instr * input
Definition: main_l2binmatch.cpp:27
CAT_cth_dalb
#define CAT_cth_dalb
Definition: l2prod.h:537
l12_proto.h
PI
#define PI
Definition: l3_get_org.c:6
CAT_cth_ctt_all
#define CAT_cth_ctt_all
Definition: l2prod.h:514
DPTB
#define DPTB(function)
Definition: passthebuck.h:24
l1que
l1qstr l1que
Definition: getl1rec.c:7
ctht_glint
float ctht_glint(float solz, float senz, float relaz, float windsp)
Definition: get_ctht.c:1635
CAT_cth_alb
#define CAT_cth_alb
Definition: l2prod.h:496
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
CAT_cth_cth_raw
#define CAT_cth_cth_raw
Definition: l2prod.h:521
comp_ctht
int32_t comp_ctht(l2str *l2rec)
Definition: get_ctht.c:295
CAT_cth_dalb_all
#define CAT_cth_dalb_all
Definition: l2prod.h:519
CAT_cth_dcod
#define CAT_cth_dcod
Definition: l2prod.h:539
CAT_cth_dcod_all
#define CAT_cth_dcod_all
Definition: l2prod.h:509
dp
float dp[MODELMAX]
Definition: atrem_corl1v2.h:248
CAT_cth_cth_all
#define CAT_cth_cth_all
Definition: l2prod.h:501
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
CAT_cth_iter
#define CAT_cth_iter
Definition: l2prod.h:523
b1
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor b1
Definition: HISTORY.txt:576
CAT_cth_lcod
#define CAT_cth_lcod
Definition: l2prod.h:520
nband
int32_t nband
Definition: pml_iop_tables.c:25
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
C_TO_K
#define C_TO_K
Definition: get_ctht.c:8
fabs
#define fabs(a)
Definition: misc.h:93
CAT_cth_dlcod
#define CAT_cth_dlcod
Definition: l2prod.h:540
CAT_cth_acost_all
#define CAT_cth_acost_all
Definition: l2prod.h:506
CAT_cth_dcth
#define CAT_cth_dcth
Definition: l2prod.h:529
BAD_INT
#define BAD_INT
Definition: genutils.h:23
CAT_cth_dctt_all
#define CAT_cth_dctt_all
Definition: l2prod.h:518
mapgen_overlay.size
size
Definition: mapgen_overlay.py:242
comp_ctht_tst
int32_t comp_ctht_tst(l1str *l1rec, int32_t proc_flg, ctht_lin_str *ctht_lin)
Definition: get_ctht.c:1289
s
data_t s[NROOTS]
Definition: decode_rs.h:75
CAT_cth_acost
#define CAT_cth_acost
Definition: l2prod.h:524
coord
Definition: dataday.h:40
invert_a_matrix
gsl_matrix * invert_a_matrix(gsl_matrix *matrix)
Definition: get_ctht.c:1722
CAT_cth_phase
#define CAT_cth_phase
Definition: l2prod.h:497
CAT_cth_cod_all
#define CAT_cth_cod_all
Definition: l2prod.h:498
i
int i
Definition: decode_rs.h:71
CAT_Cld_p
#define CAT_Cld_p
Definition: l2prod.h:424
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
CAT_cth_cth_raw_all
#define CAT_cth_cth_raw_all
Definition: l2prod.h:513
CAT_cth_dcth_all
#define CAT_cth_dcth_all
Definition: l2prod.h:512
ctht_lins
ctht_lins_str ctht_lins
Definition: get_ctht.c:10
npix
int npix
Definition: get_cmp.c:28
CAT_cth_iter_all
#define CAT_cth_iter_all
Definition: l2prod.h:507
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
init_ctht_parms
int32_t init_ctht_parms(ctht_lin_str *ctht_lin, int32_t ipix)
Definition: get_ctht.c:1247
CAT_cth_alb_all
#define CAT_cth_alb_all
Definition: l2prod.h:503
CAT_cth_cod
#define CAT_cth_cod
Definition: l2prod.h:495


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