1 # include <stdlib.h>

2 # include <stdio.h>

3 #include "l12_proto.h"

4 #include "l2prod.h"

5 #include <netcdf.h>

6 #include "met_cvt.h"

7 #include "scene_meta.h"

8 

9 extern void ch_cld_sci_( float *, int *, unsigned char *, int *,

10  int32_t *, int *, int *, char * );

11 /*

12  routine get_cmp.c - Made from tst_cld_sci.c - this will be a stub to

13  read the small chimaera data set and call the chimaera code to run

14  a small set of pixels through the science code.

15 

16  This will first demonstrate the union of l2gen and the chimaera code

17 

18  This is the c portion, which reads in the data from a netcdf file

19  wc_file.nc and set up the data to pass to the f90 code: ch_cld_sci.f90

20 

21  I have also put in 2 netcdf convenience routines based on the ncio_grab_f_ds

22  in the l2gen src area file ncio.c

23 */

24 

25 /* this must be global or the f90 won't see it */

26 struct

27  {

28  int npix; /* array sizes */

29  int nlin;

30  int nbnd;

31  int nbnd_albedo;

32  int nlvl_model;

33  int scan; /* mainly used for reading same data */

34  int st_samp; /* from MODIS L2 cloud file */

35  int g_year; /* the granule year */

36  int g_day; /* the granule day-of-year */

37  } dim_ctl_;

38 

39  static int32_t cur_cmp_rec = -1;

40 

41 int get_cmp( l2str *l2rec, int prodnum, float prod[])

42 /*

43 get_cmp will just initiate the derivation of cloud microphysical properties

44 (CMP) for the line if we haven't done so already. Then, it will pass back the

45 proper property array for that line

46 

47  Parameters: (in calling order)

48  Type Name I/O Description

49  ---- ---- --- -----------

50  l2str * l2rec I l2 record of data

51  int prodnum I Product ID to return

52  float [] prod O returned data - this storage is already allocated

53 

54 W. Robinson, SAIC, 4 Jan 2019

55 

56  W. Robinson, 14 Jun 2019 Add the cirrus reflectance data to the radiance set

57 */

58  {

59  int compute_cmp( l2str * );

60  void get_cmp_prod_( int *, float *, int * );

61  int ipix, n_prd;

62  /* if proc_cloud not set, error out */

63  if( input->proc_cloud == 0 ) {

64  printf( "%s, %d: Cloud products require the proc_cloud=1 set\n", __FILE__,

65  __LINE__ );

66  exit(1);

67  }

68  /* WDR this is just to force the call all the time */

69  if( l2rec->l1rec->iscan != cur_cmp_rec )

70  {

71  cur_cmp_rec = l2rec->l1rec->iscan;

72  if( compute_cmp( l2rec ) != 0 )

73  {

74  printf( "%s, %d: E - CMP computation failure\n", __FILE__, __LINE__ );

75  exit(1);

76  }

77  }

78  /*

79  * call get_cmp_prod to get the cmp desired

80  */

81  /* printf( "%s, %d: calling the get_cmp_prod\n", __FILE__, __LINE__ ); */

82  /* WDR this is awkward */

83 

84  n_prd = 1;

85  if( ( prodnum >= 480 ) && ( prodnum <= 482 ) ) n_prd = 10;

86  if( prodnum == 492 ) n_prd = 10;

87  get_cmp_prod_( &prodnum, prod, &n_prd );

88  /*

89  * if any product value is < -900., set prodfail

90  * For now, limit to XXX_2100 prods - some better treatment for the

91  * different cmp products is needed going forward though

92  */

93  if( ( prodnum == CAT_CER_2100 ) || ( prodnum == CAT_COT_2100 ) ||

94  ( prodnum == CAT_COT_2100 ) )

95  {

96  for( ipix = 0; ipix < l2rec->l1rec->npix; ipix++ )

97  if( *( prod + ipix ) < -900 )

98  l2rec->l1rec->flags[ipix] |= PRODFAIL;

99  }

100  /*

101  * do the same for 1600 -> PRODWARN and 1621 -> NAVFAIL

102  * 2200 -> ATMFAIL

103  */

104  if( ( prodnum == CAT_CER_1600 ) || ( prodnum == CAT_COT_1600 ) ||

105  ( prodnum == CAT_COT_1600 ) )

106  {

107  for( ipix = 0; ipix < l2rec->l1rec->npix; ipix++ )

108  if( *( prod + ipix ) < -900 )

109  l2rec->l1rec->flags[ipix] |= PRODWARN;

110  }

111  if( ( prodnum == CAT_CER_1621 ) || ( prodnum == CAT_COT_1621 ) ||

112  ( prodnum == CAT_COT_1621 ) )

113  {

114  for( ipix = 0; ipix < l2rec->l1rec->npix; ipix++ )

115  if( *( prod + ipix ) < -900 )

116  l2rec->l1rec->flags[ipix] |= NAVFAIL;

117  }

118  if( ( prodnum == CAT_CER_2200 ) || ( prodnum == CAT_COT_2200 ) ||

119  ( prodnum == CAT_COT_2200 ) )

120  {

121  for( ipix = 0; ipix < l2rec->l1rec->npix; ipix++ )

122  if( *( prod + ipix ) < -900 )

123  l2rec->l1rec->flags[ipix] |= ATMFAIL;

124  }

125  return(0);

126  }

127 

128 unsigned char *get_cmp_byt( l2str *l2rec, int prodnum )

129 /*

130 get_cmp_byt is like get_cmp but will provide the byte flag values

131 

132  Returns unsigned char array of the flag valueson the line

133 

134  Parameters: (in calling order)

135  Type Name I/O Description

136  ---- ---- --- -----------

137  l2str * l2rec I l2 record of data

138  int prodnum I Product ID to return

139 

140  W. Robinson, SAIC, 5 Aug 2019

141 

142 */

143  {

144  int compute_cmp( l2str * );

145  static unsigned char *bprod = NULL;

146  void get_cmp_byt_( int *, unsigned char * );

147  /* WDR this is just to force the call all the time */

148  if( l2rec->l1rec->iscan != cur_cmp_rec )

149  {

150  cur_cmp_rec = l2rec->l1rec->iscan;

151  if( compute_cmp( l2rec ) != 0 )

152  {

153  printf( "%s, %d: E - CMP computation failure\n", __FILE__, __LINE__ );

154  exit(1);

155  }

156  }

157 

158  if( bprod == NULL )

159  {

160  if( ( bprod = (unsigned char *)

161  malloc( l2rec->l1rec->npix * sizeof(unsigned char) ) ) == NULL )

162  {

163  printf( "%s, %d: E - unable to alocate byte storage\n",

164  __FILE__, __LINE__ );

165  exit(1);

166  }

167  }

168  /*

169  * call get_cmp_prod to get the cmp desired

170  */

171  /* printf( "%s, %d: calling the get_cmp_prod\n", __FILE__, __LINE__ ); */

172  get_cmp_byt_( &prodnum, bprod );

173 

174  return( bprod );

175  }

176 

177 int compute_cmp( l2str *l2rec )

178 /*

179  compute_cmp will derive the cloud microphysical properties

180 

181 */

182  {

183  /*

184  * set the inputs based on a netcdf file for the test

185  */

186  int set_cmp( l2str *, float **, int32_t *, int32_t **, int32_t *,

187  unsigned char **, int32_t * );

188  void get_cmp_byt_( int *, unsigned char * );

189  unsigned char *lcl_prd;

190  int32_t nfloat, nint32, nubyte;

191  static float *tdat;

192  static unsigned char *ubdat;

193  static int32_t *i32dat;

194  int sensor_id = l2rec->l1rec->l1file->sensorID;

195  int ip, npix = l2rec->l1rec->npix;

196  int cat_val = CAT_Cld_Phase_2100;

197 

198  /*

199  * assure that the cloud top height parmeters are made for this line

200  */

201  get_ctht_lin( l2rec, -1 );

202  //printf( "%s, %d: Just created the ctht prods for this scan\n", __FILE__, __LINE__ );

203  /*

204  * Fill with data from the l1que

205  */

206  if( set_cmp( l2rec, &tdat, &nfloat, &i32dat, &nint32, &ubdat, &nubyte ) != 0 )

207  return(-1);

208  /*

209  * call the chimaera code

210  */

211  /* OK, call the f90 routine */

212  ch_cld_sci_( tdat, &nfloat, ubdat, &nubyte, i32dat, &nint32, &sensor_id,

213  input->cloud_hgt_file );

214 /*

215  * OK, the cloud phase in 2.1 will be used to set the l2_flags as such:

216  * TURBIDW if phase is 2 or 4 (water) and COCCOLITH if phase is 3 (ice)

217  */

218  if( ( lcl_prd = (unsigned char *)malloc( npix * sizeof( char ) ) ) == NULL )

219  {

220  printf( "%s, %d: Unable to allocate cloud phase buffer\n", __FILE__, __LINE__ );

221  exit(1);

222  }

223  /*

224  * lastly, to set the phase at the pixels, do the following. The cloud

225  * phase 2100 applies as cloud phase for all the absorbing bands

226  */

227  get_cmp_byt_( &cat_val, lcl_prd );

228 

229  for( ip = 0; ip < npix; ip++ )

230  {

231  if( *( lcl_prd + ip ) == 3 )

232  l2rec->l1rec->flags[ip] |= COCCOLITH;

233  if( ( *( lcl_prd + ip ) == 2 ) || ( *( lcl_prd + ip ) == 4 ) )

234  l2rec->l1rec->flags[ip] |= TURBIDW;

235  }

236 

237  free(lcl_prd);

238 

239  return(0);

240  }

241 

242 int set_cmp( l2str *l2rec, float **tdat, int32_t *nfloat, int32_t **i32dat,

243  int32_t *nint32, unsigned char **ubdat, int32_t *nubyte )

244 /*

245  set_cmp will set up all the data required by the chimaera code from

246  the contents of the l1rec

247  30 Jun 2023, W Robinson, add the cloud top height arrays in (ctht)

248 

249 */

250  {

251  extern l1qstr l1que;

252  extern ctht_lins_str ctht_lins;

253  int mk_cmp_prof( float *, float *, float *, float *, float, float, float,

254  float, unsigned char *, unsigned char *, double *, double *, double *,

255  double * );

256 /* WDR **** temp to dummy the profile interp */

257  int mk_cmp_prof_dum( float *, float *, float *, float *, float, float, float,

258  float, unsigned char *, unsigned char *, double *, double *, double *,

259  double * );

260  int make_profile_101_( double * );

261  int32_t npix, nlin, nbnd_albedo, ctr_que,ibnd, iln, qln, ipx;

262  l1str *l1rec = l2rec->l1rec;

263  int32_t nbnd_ref, nbnd_emis, foff, foff2, uoff, uoff2, nbnd, sensor_id;

264  float *F0, rad, solz, out;

265  int32_t dbg_print, ix_bnd, force_oci, cld_missed;

266  int32_t ncmp_bnd = 15;

267  /* WDR a temp print routine for debug */

268  int profile_print( float *, unsigned char *, int32_t *, int32_t, int32_t,

269  int32_t, int32_t );

270 

271  /* some look-ups to id the type / location of the l2 rads relative to the

272  chimaera rads */

273  /* source for the chimaera rads: 1 reflective (Lt), 0 emmissive (Ltir)

274  and cirrus (rho_cirrus) */

275  static int32_t ref_for_cmb[] = { 1, 1, 1, 1, 1, 1, 0, 2, 0, 0, 1, 0, 0, 1, 1 };

276  /* cloud band is in l1rec if 1, else 0 */

277  static char cmp_there[15];

278  /* index in Lt (ref) or Ltir (emis) (after # vis bands subtracted) */

279  static int32_t l1ix[15];

280  /* cloud band ideal wavelengths for MODIS (and any other inst to try) */

281  int32_t cld_wav_mod[] = { 645, 859, 1240, 1640, 2130, 930, 3750, 1375,

282  8550, 11000, 443, 12000, 13600, 469, 555 };

283  /* cloud bands for OCI (only) */

284  // int32_t cld_wav_oci[] = { 645, 859, 1250, 1615, 2130, 930, 3750, 1375,

285  // 8550, 11000, 443, 12000, 13600, 469, 555 };

286  /* 10 Sep 21 to match the waves for the Kerry Meyer table */

287  /* 29 aug 22 the 8.8 um not in OCI, replace with the 2.26 um */

288  /* WDR pre-2.2 band use in OCI

289  int32_t cld_wav_oci[] = { 665, 865, 1250, 1616, 2130, 930, 3750, 1375,

290  2260, 11000, 443, 12000, 13600, 469, 555 };

291  */

292  /* WDR 18 Sep 23 set cld_wav_oci as cld_wav_ocis and make cld_wav_oci

293  the actual instrument wavelengths

294  int32_t cld_wav_oci[] = { 665, 865, 1250, 1616, 2130, 930, 2260, 1375,

295  8550, 11000, 443, 12000, 13600, 469, 555 };

296  */

297  // for inst OCIS

298  int32_t cld_wav_ocis[] = { 665, 865, 1250, 1616, 2130, 930, 2260, 1375,

299  8550, 11000, 443, 12000, 13600, 469, 555 };

300  // for inst OCI

301  int32_t cld_wav_oci[] = { 665, 865, 1249, 1618, 2131, 940, 2258, 1378,

302  8550, 11000, 442, 12000, 13600, 470, 555 };

303  int32_t cld_wav[15];

304  /* OCI bands available (for simulation) */

305  int32_t oci_bnd_avail[] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1 };

306  /* minimum bands required for cloud processing */

307  static int32_t cld_min_bnd[] = { 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1 };

308  /* The following emperically derived factors will correct the l2gen radiances

309  and reflectances to those that chimaera expects - probably the difference

310  between OBPG and EOSDIS calibration */

311  // *** I have temporarily removed the can change for a test

312  //float chim_rad_corr[] = { 0.965407, 0.973168, 1.00678, 1.00244, 0.963511,

313  float chim_rad_corr[] = { 1.0, 1.0, 1.0, 1.0, 1.0,

314  // WDR Jan 2021 try 5% low in 1st 5 bands

315  //float chim_rad_corr[] = { 0.5, 0.5, 0.5, 0.5, 0.5,

316  /* bands 645 859 1240 1640 2130 */

317  //-1., 0.999840, -1., 0.999979, 0.999794, 1.01814, 0.999760,

318  -1., 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,

319  /* 930 approx 3750 1375 approx 8550 11000 412 12000 */

320  //-1., 0.986288, 0.994909 };

321  -1., 1.0, 1.0 };

322  /* 13600 approx 469 555 */

323 

324  /* Guard the transfer array creation to only 1 per run */

325  static int32_t firstcall = 0;

326  /*

327  * for the met data,

328  * set the merra profile heights in reverse order to match the order in the

329  * profiles used by chimaera

330  */

331  int32_t nlvl = 42, nlvl_cmp = 101, ilvl, is_land, cld_flg;

332  int32_t bad_rad, ice_flg, glint_flg, n_sfc_albedo;

333  int32_t npix_scan, lst, len, pst, pen, apx, aln, itot, nland;

334  unsigned char sfc_lvl, trop_lvl;

335  float merra_p_set[] = { 0.1, 0.3, 0.4, 0.5, 0.7, 1., 2., 3., 4., 5., 7., 10.,

336  20., 30., 40., 50., 70., 100., 150., 200., 250., 300., 350., 400.,

337  450., 500., 550., 600., 650., 700., 725., 750., 775., 800., 825.,

338  850., 875., 900., 925., 950., 975., 1000. };

339  float merra_t[nlvl], merra_p[nlvl], merra_q[nlvl], merra_h[nlvl];

340  double cmp_t[nlvl_cmp], cmp_mixr[nlvl_cmp], cmp_h[nlvl_cmp];

341  double cmp_p[nlvl_cmp], cmp_p_set[nlvl_cmp];

342  scnstr *meta;

343  /*

344  get the # pixels

345  and for now, make sure they match the # coming in

346  */

347  ctr_que = l1que.nq / 2; /* the center of the queue = our primary line */

348 

349  npix = l1rec->npix;

350  // WDR test nlin = 3;

351  nlin = 3;

352  dbg_print = 0;

353 

354  dim_ctl_.npix = npix;

355  dim_ctl_.nbnd_albedo = 6; /* the MODIS # 'albedo bands' */

356  nbnd_albedo = dim_ctl_.nbnd_albedo;

357  dim_ctl_.nlin = nlin;

358  dim_ctl_.nbnd = ncmp_bnd; /* a bnd # used from MODIS */

359  nbnd = ncmp_bnd;

360  dim_ctl_.nlvl_model = nlvl_cmp; /* the set # interpolated to for chimaera */

361  dim_ctl_.scan = l1que.r[ctr_que].iscan;

362  dim_ctl_.st_samp = 1; /* probably not used, but ftn 1-origin start samp */

363  meta = scene_meta_get();

364  /* get the time info down to the day - all that is used */

365  dim_ctl_.g_year = meta->start_year;

366  dim_ctl_.g_day = meta->start_day;

367 

368  sensor_id = l1rec->l1file->sensorID;

369  n_sfc_albedo = 5;

370  if( ( sensor_id == OCIS ) || ( sensor_id == OCI ) ) n_sfc_albedo = 6;

371 

372  /* Allocate arrays to pass to the f90 routine with all the information

373  together: all real arrays into a real array etc */

374  *nfloat = npix * nlin * nbnd + /* reflectance */

375  npix * nlin * nbnd_albedo + /* band uncertainty */

376  2 * nbnd_albedo + /* uncertainty info */

377  7 * npix * nlin + /* geom info */

378  4 * npix * nlin * nlvl_cmp + /* met profiles */

379  8 * npix * nlin + /* 2D anc data */

380  n_sfc_albedo * npix * nlin + /* surface albedo for 5/6 bands */

381  3 * npix * nlin; /* for the 3 arrays of cloud top height, pressure, and

382  temperature */

383 

384  nbnd_ref = l1rec->l1file->nbands;

385  nbnd_emis = l1rec->l1file->nbandsir;

386  if( firstcall == 0 )

387  {

388  firstcall = 1;

389  if( ( *tdat = (float *) malloc( *nfloat * sizeof(float) ) ) == NULL )

390  {

391  printf( "%s - %d: Allocation of real storage failed\n",

392  __FILE__, __LINE__ );

393  exit(3);

394  }

395  /*

396  * : all byte arrays

397  */

398  *nubyte = 2 * npix * nlin + /* for the ancillary byte stuff */

399  23 * npix * nlin + /* for the cloud mask stuff */

400  ncmp_bnd; /* for the band availability array */

401  if( ( *ubdat = (unsigned char *)

402  malloc( *nubyte * sizeof(unsigned char) ) ) == NULL )

403  {

404  printf( "%s - %d: Allocation of byte storage failed\n",

405  __FILE__, __LINE__ );

406  exit(3);

407  }

408  /*

409  * : all int32 arrays

410  */

411  *nint32 = npix * nlin;

412  if( ( *i32dat = (int32_t *) malloc( *nint32 * sizeof(int32_t) ) ) == NULL )

413  {

414  printf( "%s - %d: Allocation of i32 storage failed\n",

415  __FILE__, __LINE__ );

416  exit(3);

417  }

418  /*

419  * Set up the cloud band presence and L1 index here

420  */

421  for( ibnd = 0; ibnd < ncmp_bnd; ibnd++ )

422  {

423  if( sensor_id == OCIS )

424  cld_wav[ibnd] = cld_wav_ocis[ibnd];

425  else if( sensor_id == OCI )

426  cld_wav[ibnd] = cld_wav_oci[ibnd];

427  else

428  cld_wav[ibnd] = cld_wav_mod[ibnd];

429  }

430  /* if OCI, make the 2.2 reflective band */

431  if( ( sensor_id == OCIS ) || ( sensor_id == OCI ) )

432  {

433  ref_for_cmb[6] = 1;

434  cld_min_bnd[6] = 1;

435  }

436  for( ibnd = 0; ibnd < ncmp_bnd; ibnd++ )

437  {

438  ix_bnd = bindex_get( cld_wav[ibnd] );

439  if( ( ix_bnd < 0 ) && ( ref_for_cmb[ibnd] != 2 ) ) // if no matching,

440  // non-cirrus band was found

441  {

442  cmp_there[ibnd] = 0;

443  l1ix[ibnd] = -1;

444  }

445  else

446  {

447  cmp_there[ibnd] = 1;

448  if( ref_for_cmb[ibnd] == 0 ) // emissive

449  l1ix[ibnd] = ix_bnd;

450  else if( ref_for_cmb[ibnd] == 1 ) // reflective

451  l1ix[ibnd] = ix_bnd;

452  else //cirrus

453  l1ix[ibnd] = 0;

454  }

455  }

456  /*

457  * This code will optionally force the available bands down to the OCI bands

458  */

459  force_oci = 0;

460  printf( "OCI band forcing is set to %d\n", force_oci );

461  if( force_oci == 1 )

462  {

463  for( ibnd = 0; ibnd < ncmp_bnd; ibnd++ )

464  cmp_there[ibnd] *= oci_bnd_avail[ibnd];

465  }

466  /*

467  * If the minimum band set for 2.1, (2.2, OCI) and 1.6 nm cloud products

468  * (cld_min_bnd) is not in this instrument (cmp_there), report it and

469  * stop processing here

470  */

471  printf( "Talley of current instrument's wavelengths\n" );

472  printf( "Cloud wave is wave is\n" );

473  printf( "wavelength available required\n" );

474  cld_missed = 0;

475  for( ibnd = 0; ibnd < ncmp_bnd; ibnd++ )

476  {

477  if( ( cld_min_bnd[ibnd] == 1 ) && ( cmp_there[ibnd] == 0 ) )

478  cld_missed = 1;

479  printf( "%10d %9d %8d\n", cld_wav[ibnd], cmp_there[ibnd],

480  cld_min_bnd[ibnd] );

481  }

482  if( cld_missed == 1 )

483  {

484  printf( "%s, %d: The required wavelengths for cloud processing were not found, Exiting\n", __FILE__, __LINE__ );

485  exit(FATAL_ERROR);

486  }

487  }

488 

489  F0 = l1rec->Fo;

490 

491  /* go through the chimaera refl (refl) rad (emiss) bands */

492  for( ibnd = 0; ibnd < ncmp_bnd; ibnd++ )

493  {

494  for( iln = 0; iln < nlin; iln++ )

495  {

496  qln = iln + ctr_que - nlin / 2;

497  for( ipx = 0; ipx < npix; ipx++ )

498  {

499  /* only write the center line value */

500  if( ( iln == 1) && ( dbg_print == 1 ) )

501  {

502  printf( "%s, %d: bnd %d, lin %d, pix %d, rad %f\n",

503  __FILE__, __LINE__, ibnd, iln, ipx,

504  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] );

505  }

506  /* only substitute if the band exists (for now) */

507  if( cmp_there[ibnd] == 1 )

508  {

509  /* for the reflective */

510  if( ref_for_cmb[ibnd] == 1 )

511  {

512  rad = l1que.r[qln].Lt[ l1ix[ ibnd ] + nbnd_ref * ipx ];

513  solz = l1que.r[qln].solz[ipx ];

514  out = rad * PI / ( cos( solz * PI / 180. ) * F0[ l1ix[ ibnd ] ] );

515  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] =

516  chim_rad_corr[ibnd] * out;

517  }

518  else if( ref_for_cmb[ibnd] == 2 ) /* for the cirrus */

519  {

520  solz = l1que.r[qln].solz[ipx ];

521  rad = l1que.r[qln].rho_cirrus[ipx];

522  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] =

523  chim_rad_corr[ibnd] * rad * cos( solz * PI / 180. );

524  }

525  else /* for emissive */

526  {

527  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] =

528  chim_rad_corr[ibnd] *

529  10. * l1que.r[qln].Ltir[ ( l1ix[ ibnd ] - nbnd_ref ) +

530  nbnd_emis * ipx ];

531  }

532  }

533  /* for the missing... try -999.0 */

534  else

535  {

536  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] = -999.0;

537  }

538  /* repeat write of the center line value */

539  if( ( iln == 1 ) && (dbg_print == 1 ) )

540  {

541  printf( "%s, %d: bnd %d, lin %d, pix %d, rad %f\n",

542  __FILE__, __LINE__, ibnd, iln, ipx,

543  (*tdat)[ ipx + npix * ( ibnd + ncmp_bnd * iln ) ] );

544  }

545  }

546  }

547  }

548  /*

549  * the uncertainty with 6 band_albedo - I have a little on this:

550  * There is a equation in modis_science_module.f90 that shows how the

551  * info is used:

552  * refl_unc = spec_uncertain * exp(band_uncertainty) / uncertain_sf

553  * This is a kind of odd exponential relation.

554  * It also shows that I might be better using 1. for the dummy uncertain_sf

555  * value

556  *

557  * so 0 values are OK for now except the uncertain_sf (used in denominator)

558  *

559  * *** NOTE *** that the band_uncertainty in the chimaera code is

560  * actually a unsigned byte or i*1. Now, assigning the real data to that

561  * in the cloud code should work just fine - just some overkill in the

562  * value description. However, making this data a char array would be

563  * work for very little gain and could make problems if done wrong. It

564  * could also hide the problem we're looking for too. So, I'm leaving

565  * this as-is for now.

566  */

567  foff = npix * nlin * ncmp_bnd; /* offset to the output data */

568  for( ibnd = 0; ibnd < nbnd_albedo; ibnd++ )

569  {

570  for( iln = 0; iln < nlin; iln++ )

571  {

572  qln = iln + ctr_que - nlin / 2;

573  for( ipx = 0; ipx < npix; ipx++ )

574  {

575  /* only write the center line value */

576  if( ( iln == 1 ) && (dbg_print == 1 ) )

577  {

578  printf( "%s, %d: bnd %d, lin %d, pix %d, UNC %f\n",

579  __FILE__, __LINE__, ibnd, iln, ipx,

580  (*tdat)[ foff + ipx + npix * ( ibnd + nbnd_albedo * iln ) ] );

581  }

582  /* substitute a 0 for all */

583  (*tdat)[ foff + ipx + npix * ( ibnd + nbnd_albedo * iln ) ] = 0.;

584 

585  /* repeat write of the center line value */

586  if( ( iln == 1 ) && (dbg_print == 1 ) )

587  {

588  printf( "%s, %d: bnd %d, lin %d, pix %d, UNC %f\n",

589  __FILE__, __LINE__, ibnd, iln, ipx,

590  (*tdat)[ foff + ipx + npix * ( ibnd + nbnd_albedo * iln ) ] );

591  }

592  }

593  }

594  }

595  /*

596  * also the nbnd_albedo long arrays of spec_uncertain and uncertain_sf see

597  * above

598  */

599  foff += npix * nlin * nbnd_albedo;

600  foff2 = foff + nbnd_albedo;

601  for( ibnd = 0; ibnd < nbnd_albedo; ibnd++ )

602  {

603  (*tdat)[ foff + ibnd ] = 0.;

604  (*tdat)[ foff2 + ibnd ] = 1.;

605  }

606  /*

607  * transfer the geolocation and view angle information over

608  */

609  foff = foff2 + nbnd_albedo;

610  foff2 = npix * nlin;

611  /* lat and lon at the l1rec values */

612  for( iln = 0; iln < nlin; iln++ )

613  {

614  qln = iln + ctr_que - nlin / 2;

615  for( ipx = 0; ipx < npix; ipx++ )

616  {

617  if( ( iln == 1 ) && ( dbg_print == 1 ) )

618  {

619  printf( "%s, %d: lin %d, pix %d, INITIAL lon %f, lat %f, senz %f, sena %f, solz %f, sola %f, relaz %f\n",

620  __FILE__, __LINE__, iln, ipx,

621  (*tdat)[ foff + ipx + npix * iln ], (*tdat)[ foff + foff2 + ipx + npix * iln ],

622  (*tdat)[ foff + 2 * foff2 + ipx + npix * iln ],

623  (*tdat)[ foff + 3 * foff2 + ipx + npix * iln ],

624  (*tdat)[ foff + 4 * foff2 + ipx + npix * iln ],

625  (*tdat)[ foff + 5 * foff2 + ipx + npix * iln ],

626  (*tdat)[ foff + 6 * foff2 + ipx + npix * iln ] );

627  }

628  /* lat and lon */

629  out = l1que.r[qln].lat[ipx];

630  (*tdat)[ foff + ipx + npix * iln ] = out;

631  out = l1que.r[qln].lon[ipx];

632  (*tdat)[ foff + foff2 + ipx + npix * iln ] = out;

633  /* senz, sena, sola, solz, relaz */

634  out = l1que.r[qln].senz[ipx];

635  (*tdat)[ foff + 2 * foff2 + ipx + npix * iln ] = out;

636  out = l1que.r[qln].sena[ipx];

637  (*tdat)[ foff + 3 * foff2 + ipx + npix * iln ] = out;

638  out = l1que.r[qln].sola[ipx];

639  (*tdat)[ foff + 4 * foff2 + ipx + npix * iln ] = out;

640  out = l1que.r[qln].solz[ipx];

641  (*tdat)[ foff + 5 * foff2 + ipx + npix * iln ] = out;

642  out = l1que.r[qln].delphi[ipx];

643  (*tdat)[ foff + 6 * foff2 + ipx + npix * iln ] = -out; /* their relaz is (-) ours */

644 

645  if( ( iln == 1 ) && ( dbg_print == 1 ) )

646  {

647  printf( "%s, %d: lin %d, pix %d, FINAL lon %f, lat %f, senz %f, sena %f, solz %f, sola %f, relaz %f\n",

648  __FILE__, __LINE__, iln, ipx,

649  (*tdat)[ foff + ipx + npix * iln ], (*tdat)[ foff + foff2 + ipx + npix * iln ],

650  (*tdat)[ foff + 2 * foff2 + ipx + npix * iln ],

651  (*tdat)[ foff + 3 * foff2 + ipx + npix * iln ],

652  (*tdat)[ foff + 4 * foff2 + ipx + npix * iln ],

653  (*tdat)[ foff + 5 * foff2 + ipx + npix * iln ],

654  (*tdat)[ foff + 6 * foff2 + ipx + npix * iln ] );

655  }

656  }

657  }

658  /*

659  * next is the met info

660  */

661  foff += 7 * foff2; /* start of the profiles */

662  /* set up the clean hi res profile */

663  make_profile_101_( cmp_p_set );

664 

665  /* check that profiles were read in */

666  if( l1que.r[ctr_que].anc_add == NULL )

667  {

668  fprintf( stderr,

669  "-E- %s %d: Ancillary profile required for cloud processing is missing\n",

670  __FILE__, __LINE__);

671  exit(1);

672  }

673  uoff = 0;

674  for( iln = 0; iln < nlin; iln++ )

675  {

676  qln = iln + ctr_que - nlin / 2;

677  for( ipx = 0; ipx < npix; ipx++ )

678  {

679  memcpy( merra_p, merra_p_set, nlvl * sizeof(float) );

680  for( ilvl = 0; ilvl < nlvl; ilvl++ )

681  {

682  /* note that the profiles for chimaera are reverse of the MERRA,

683  so in making the merra arrays from the l1rec stuff, we reverse it */

684  merra_t[ nlvl - 1 - ilvl ] = l1que.r[qln].anc_add->prof_temp[ ilvl + nlvl * ipx ];

685  merra_q[ nlvl - 1 - ilvl ] = l1que.r[qln].anc_add->prof_q[ ilvl + nlvl * ipx ];

686  merra_h[ nlvl - 1 - ilvl ] = l1que.r[qln].anc_add->prof_height[ ilvl + nlvl * ipx ];

687  }

688 /* WDR trad test point tpix = 2, tlin = 1; */

689 int tpix = 2, tlin = 1;

690 if( (ipx == tpix) && ( iln == tlin ) && ( dbg_print == 1 ) )

691  {

692  printf( "\n\n%s, %d: merra_p for pix: %d, lin: %d\n", __FILE__, __LINE__, ipx, iln );

693  for( ilvl = 0; ilvl < nlvl; ilvl++ )

694  printf( "%f ", merra_p[ilvl] );

695  printf( "\n" );

696 

697  printf( "%s, %d: merra_t for pix: %d, lin: %d\n", __FILE__, __LINE__, ipx, iln );

698  for( ilvl = 0; ilvl < nlvl; ilvl++ )

699  printf( "%f ", merra_t[ilvl] );

700  printf( "\n" );

701 

702  printf( "%s, %d: merra_q for pix: %d, lin: %d\n", __FILE__, __LINE__, ipx, iln );

703  for( ilvl = 0; ilvl < nlvl; ilvl++ )

704  printf( "%f ", merra_q[ilvl] );

705  printf( "\n" );

706 

707  printf( "%s, %d: merra_h for pix: %d, lin: %d\n", __FILE__, __LINE__, ipx, iln );

708  for( ilvl = 0; ilvl < nlvl; ilvl++ )

709  printf( "%f ", merra_h[ilvl] );

710  printf( "\n" );

711  }

712 

713  /*

714  * for this profile, make the 101 level chimaera profile

715  * and put the chimaera profiles into the storage

716  */

717  memcpy( cmp_p, cmp_p_set, nlvl_cmp * sizeof( double) );

718 /* WDR **** test branch to get around interpolation of merra */

719 /* ilvl = -899; WDR use a pre-computed fixed profile */

720 ilvl = 0; /* WDR use the real profiles */

721 if( ilvl == -899 )

722  {

723  mk_cmp_prof_dum( merra_p, merra_t, merra_q, merra_h, l1que.r[qln].sfcp[ipx],

724  l1que.r[qln].sfct[ipx], l1que.r[qln].sfcrh[ipx],

725  l1que.r[qln].lat[ipx], &sfc_lvl, &trop_lvl, cmp_p, cmp_t,

726  cmp_mixr, cmp_h );

727  }

728 else

729  {

730  mk_cmp_prof( merra_p, merra_t, merra_q, merra_h, l1que.r[qln].sfcp[ipx],

731  l1que.r[qln].sfct[ipx], l1que.r[qln].sfcrh[ipx],

732  l1que.r[qln].lat[ipx], &sfc_lvl, &trop_lvl, cmp_p, cmp_t,

733  cmp_mixr, cmp_h );

734  }

735  /* as sfc_lvl, trop_lvl go to ftn, make them 1-origin */

736  sfc_lvl = ( sfc_lvl == nlvl_cmp ) ? sfc_lvl : sfc_lvl + 1;

737  trop_lvl = ( trop_lvl == nlvl_cmp ) ? trop_lvl : trop_lvl + 1;

738 

739  /*

740  * WDR Look through the profiles and see if all is OK with them.

741  * out, no need for debug

742  *** 5 Mar 2024 I may check this below in mk_cmp_prof

743  *** but comment out for now

744  */

745 /*

746  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

747  {

748  if( ( cmp_t[ilvl] > 350. ) || ( cmp_t[ilvl] < 100. ) )

749  {

750  printf( "PROFILE ODD T: %f, iln: %d, ipx: %d, ilvl: %d\n",

751  cmp_t[ilvl], iln, ipx, ilvl );

752  }

753  if( ( cmp_mixr[ilvl] > 100. ) || ( cmp_mixr[ilvl] < 0 ) )

754  {

755  printf( "PROFILE ODD MIXR: %f, iln: %d, ipx: %d, ilvl: %d\n",

756  cmp_mixr[ilvl], iln, ipx, ilvl );

757  }

758  if( ( cmp_h[ilvl] > 1000000. ) || ( cmp_h[ilvl] < -20000. ) )

759  {

760  printf( "PROFILE ODD height: %f, iln: %d, ipx: %d, ilvl: %d\n",

761  cmp_h[ilvl], iln, ipx, ilvl );

762  }

763  if( ( cmp_p[ilvl] > 1200. ) || ( cmp_p[ilvl] < 0. ) )

764  {

765  printf( "PROFILE ODD P: %f, iln: %d, ipx: %d, ilvl: %d\n",

766  cmp_p[ilvl], iln, ipx, ilvl );

767  }

768  }

769  end commenting */

770 /* the tdat at foff starts with mixr

771  it runs all pixels fastest, then lines, themn levels

772  so t is at foff + npix * nlin * nlvl_cmp

773  p is at foff + 2 * nlin * nlvl_cmp

774  h is at foff + 3 * nlin * nlvl_cmp

775 

776  the offset for pix, lin level is the base above plus

777  ipx * npix * ( iln + nlin * ilvl )

778 */

779 /* I need to print the profiles before they get re-done and then again after

780  they are updated. Maybe a routine for that would be usefull */

781 if( ( ipx == tpix ) && ( iln == tlin ) && ( dbg_print == 1 ) )

782  {

783  printf( "/n%s, %d: mix, t, p, h profiles BEFORE, pix %d, lin %d\n", __FILE__, __LINE__, ipx, iln );

784  profile_print( (*tdat) + foff, (*ubdat), (*i32dat), ipx, iln, npix, nlin );

785  printf( "OUR sfcp: %f, sfct: %f\n", l1que.r[qln].sfcp[ipx], l1que.r[qln].sfct[ipx] );

786  printf( "mk_cmp_prof sfc_lvl: %d, trop_lvl: %d\n", sfc_lvl, trop_lvl );

787  }

788  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

789  {

790  (*tdat)[foff + ipx + npix * ( iln + nlin * ilvl ) ] = cmp_mixr[ilvl];

791  foff2 = foff + npix * nlin * nlvl_cmp;

792  (*tdat)[foff2 + ipx + npix * ( iln + nlin * ilvl ) ] = cmp_t[ilvl];

793  foff2 = foff + 2 * npix * nlin * nlvl_cmp;

794  (*tdat)[foff2 + ipx + npix * ( iln + nlin * ilvl ) ] = cmp_p[ilvl];

795  foff2 = foff + 3 * npix * nlin * nlvl_cmp;

796  (*tdat)[foff2 + ipx + npix * ( iln + nlin * ilvl ) ] = cmp_h[ilvl];

797  }

798  /*

799  * add the single-level met parameters

800  * sfc T

801  */

802  foff2 = foff + 4 * npix * nlin * nlvl_cmp;

803  (*tdat)[foff2 + ipx + npix * iln ] = 273.15 + l1que.r[qln].sfct[ipx];

804  /* sfc P */

805  foff2 += npix * nlin;

806  (*tdat)[foff2 + ipx + npix * iln ] = l1que.r[qln].sfcp[ipx];

807  /* wind speed */

808  foff2 += npix * nlin;

809  (*tdat)[foff2 + ipx + npix * iln ] = l1que.r[qln].ws[ipx];

810  /* O3 */

811  foff2 += npix * nlin;

812  (*tdat)[foff2 + ipx + npix * iln ] = 1.e3 * l1que.r[qln].oz[ipx];

813  /* ice frac and snow frac */

814  is_land = l1que.r[qln].land[ipx];

815 /*

816 if( ( ipx == 2 ) && ( iln == 1 ) )

817 printf( "flag: %d\n", l1que.r[qln].flags[ipx]);

818 */

819  foff2 += npix * nlin;

820  (*tdat)[foff2 + ipx + npix * iln ] = ( is_land ) ? 0. :

821  l1que.r[qln].icefr[ipx];

822 

823  foff2 += npix * nlin;

824  (*tdat)[foff2 + ipx + npix * iln ] = ( is_land ) ?

825  l1que.r[qln].icefr[ipx] : 0;

826  /* alt o3 - originally an alt o3 source but... */

827  foff2 += npix * nlin;

828  (*tdat)[foff2 + ipx + npix * iln ] = 1.e3 * l1que.r[qln].oz[ipx];

829  /* alt ice conc */

830  foff2 += npix * nlin;

831  /* WDR to globalize this value

832  (*tdat)[foff2 + ipx + npix * iln ] = ( is_land ) ? 0. :

833  l1que.r[qln].icefr[ipx];

834  */

835  (*tdat)[foff2 + ipx + npix * iln ] = l1que.r[qln].icefr[ipx];

836  /* We will get the surface albedo from l2gen now */

837  foff2 += npix * nlin;

838  (*tdat)[foff2 + ipx + npix * iln ] =

839  l1que.r[qln].cld_dat->sfc_albedo_659[ipx];

840  foff2 += npix * nlin;

841  (*tdat)[foff2 + ipx + npix * iln ] =

842  l1que.r[qln].cld_dat->sfc_albedo_858[ipx];

843  foff2 += npix * nlin;

844  (*tdat)[foff2 + ipx + npix * iln ] =

845  l1que.r[qln].cld_dat->sfc_albedo_1240[ipx];

846  foff2 += npix * nlin;

847  (*tdat)[foff2 + ipx + npix * iln ] =

848  l1que.r[qln].cld_dat->sfc_albedo_1640[ipx];

849  foff2 += npix * nlin;

850  (*tdat)[foff2 + ipx + npix * iln ] =

851  l1que.r[qln].cld_dat->sfc_albedo_2130[ipx];

852  /* WDR for 2.2 as we don't have 2.2, duplicate the 2.1 here */

853  if( ( sensor_id == OCIS ) || ( sensor_id == OCI ) )

854  {

855  foff2 += npix * nlin;

856  (*tdat)[foff2 + ipx + npix * iln ] =

857  l1que.r[qln].cld_dat->sfc_albedo_2130[ipx];

858  }

859  /* add the cloud top height, pressure, and temperature */

860  foff2 += npix * nlin;

861  (*tdat)[foff2 + ipx + npix * iln ] =

862  ctht_lins.ct[qln]->cth[ipx];

863  foff2 += npix * nlin;

864  (*tdat)[foff2 + ipx + npix * iln ] =

865  ctht_lins.ct[qln]->ctp[ipx];

866  foff2 += npix * nlin;

867  (*tdat)[foff2 + ipx + npix * iln ] =

868  ctht_lins.ct[qln]->ctt[ipx];

869  /* end of float data set up */

870 

871  /* sfc_lvl and trop_lvl gotten with the profile above */

872  /* note that this is unsigned char data now */

873  /* also note that the ftn expects 1-origin levels so add 1 */

874  (*ubdat)[ uoff + ipx + npix * iln ] = sfc_lvl;

875  uoff2 = uoff + npix * nlin;

876  (*ubdat)[ uoff2 + ipx + npix * iln ] = trop_lvl;

877 

878  /* the I32 met info is the next */

879  /* just try inserting the ice fraction as percent */

880  (*i32dat)[ ipx + npix * iln ] = 100 * l1que.r[qln].icefr[ipx];

881 

882  /* the remaining byte data is all from the cloud mask that is

883  a MODIS product. We have one thing: the cloud mask

884  and we'll see how that does... */

885  /* some set-up and CLD_DET */

886  uoff2 += npix * nlin;

887  bad_rad = ( l1que.r[qln].Lt[l1ix[0]+nbnd_ref * ipx] == BAD_INT ) ? 1 : 0;

888  cld_flg = l1que.r[qln].cloud[ipx];

889  ice_flg = l1que.r[qln].ice[ipx];

890  glint_flg = l1que.r[qln].glint[ipx];

891  npix_scan = l1que.r[qln].npix;

892  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( bad_rad == 0 ) ? 1 : 0;

893  // We will use the cloud mask file categories if the file was specified

894  if( l1_input->cld_msk_file[0]) {

895  char cld_category;

896  get_sdps_cld_mask( &(l1que.r[qln]), ipx, &cld_category );

897  /* for missing cloud, set bad_rad, no cloud detected and make clear */

898  if( cld_category == BAD_BYTE ){

899  bad_rad = 1;

900  cld_category = 3;

901  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

902  }

903  /* CONF_CLD, CLR_66, CLR_95, CLR_99 */

904  uoff2 += npix * nlin;

905  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_category == 0 ) ? 1 : 0;

906  uoff2 += npix * nlin;

907  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_category == 1 ) ? 1 : 0;

908  uoff2 += npix * nlin;

909  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_category == 2 ) ? 1 : 0;

910  uoff2 += npix * nlin;

911  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_category >= 3 ) ? 1 : 0;

912  } else {

913  // use the binary cloud flag

914  /* CONF_CLD, CLR_66, CLR_95, CLR_99 */

915  uoff2 += npix * nlin;

916  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_flg != 0 ) ? 1 : 0;

917  uoff2 += npix * nlin;

918  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_flg != 0 ) ? 1 : 0;

919  uoff2 += npix * nlin;

920  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( cld_flg != 0 ) ? 1 : 0;

921  uoff2 += npix * nlin;

922  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

923  }

924  /* SNO_SFC */

925  uoff2 += npix * nlin;

926  (*ubdat)[ uoff2 + ipx + npix * iln ] =

927  ( ( is_land == 1 ) && ( ice_flg == 1 ) ) ? 1 : 0;

928  /* WTR_SFC */

929  uoff2 += npix * nlin;

930  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( is_land == 0 ) ? 1 : 0;

931  /* COAST_SFC */

932  uoff2 += npix * nlin;

933  lst = ( iln < 1 ) ? iln : iln - 1;

934  len = ( iln >= ( nlin - 1 ) ) ? iln : iln + 1;

935  pst = ( ipx <= 0 ) ? 0 : ipx -1;

936  pen = ( ipx >= ( npix_scan - 1 ) ) ? ipx : ipx + 1;

937  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

938  itot = 0;

939  nland = 0;

940  for( aln = lst; aln <= len; aln++ )

941  for( apx = pst; apx <= pen; apx++ )

942  {

943  itot++;

944  if( l1que.r[aln].land[apx] ) nland++;

945  }

946  if( ( nland != 0 ) && ( nland != itot ) )

947  (*ubdat)[ uoff2 + ipx + npix * iln ] = 1;

948  /* DESERT_SFC */

949  uoff2 += npix * nlin;

950  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

951  /* LND_SFC */

952  uoff2 += npix * nlin;

953  (*ubdat)[ uoff2 + ipx + npix * iln ] = ( is_land ) ? 1 : 0;

954  /* NIGHT - solz > 90 */

955  uoff2 += npix * nlin;

956  (*ubdat)[ uoff2 + ipx + npix * iln ] =

957  ( l1que.r[iln].solz[ipx] > 90. ) ? 1 : 0;

958  /* GLINT */

959  uoff2 += npix * nlin;

960  (*ubdat)[ uoff2 + ipx + npix * iln ] =

961  ( glint_flg ) ? 1 : 0;

962  /* OCEAN_NO_SNOW */

963  uoff2 += npix * nlin;

964  (*ubdat)[ uoff2 + ipx + npix * iln ] =

965  ( ( is_land == 0 ) && ( ice_flg == 0 ) ) ? 1 : 0;

966  /* OCEAN_SNOW */

967  uoff2 += npix * nlin;

968  (*ubdat)[ uoff2 + ipx + npix * iln ] =

969  ( ( is_land == 0 ) && ( ice_flg == 1 ) ) ? 1 : 0;

970  /* LND_NO_SNOW */

971  uoff2 += npix * nlin;

972  (*ubdat)[ uoff2 + ipx + npix * iln ] =

973  ( ( is_land == 1 ) && ( ice_flg == 0 ) ) ? 1 : 0;

974  /* LND_SNOW */

975  uoff2 += npix * nlin;

976  (*ubdat)[ uoff2 + ipx + npix * iln ] =

977  ( ( is_land == 1 ) && ( ice_flg == 1 ) ) ? 1 : 0;

978  /* TST_H_138 */

979  uoff2 += npix * nlin;

980  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

981  /* TST_VIS_REFL */

982  uoff2 += npix * nlin;

983  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

984  /* TST_VIS_RATIO */

985  uoff2 += npix * nlin;

986  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

987  /* VIS_CLD_250 */

988  uoff2 += npix * nlin;

989  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

990  /* APPL_HCLD_138 */

991  uoff2 += npix * nlin;

992  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

993  /* APPL_VIS_REFL */

994  uoff2 += npix * nlin;

995  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

996  /* APPL_VIS_NIR_RATIO */

997  uoff2 += npix * nlin;

998  (*ubdat)[ uoff2 + ipx + npix * iln ] = 0;

999  /* end all inputs */

1000 

1001 if( ( ipx == tpix ) && ( iln == tlin ) && ( dbg_print == 1 ) )

1002  {

1003  printf( "/n%s, %d: mix, t, p, h profiles AFTER, pix %d, lin %d\n", __FILE__, __LINE__, ipx, iln );

1004  profile_print( (*tdat) + foff , (*ubdat), (*i32dat), ipx, iln, npix, nlin );

1005  }

1006  /*

1007  * note that we also get the trop_lvl, sfc_lvl which go into another

1008  * storage area further on in the work

1009  */

1010  }

1011  }

1012  /*

1013  * add the cloud band presence array to the end of the byte array

1014  */

1015  memcpy( ( *ubdat + 25 * npix * nlin ), cmp_there, ncmp_bnd );

1016 

1017  return(0);

1018  }

1019 

1020 int profile_print( float *prof, unsigned char *ubdat, int32_t *i32dat,

1021  int32_t ipx, int32_t iln, int32_t npix, int32_t nlin )

1022  /*******************************************************************

1023 

1024  profile_print

1025 

1026  purpose: print the met profile and 2d data at a point

1027 

1028  Returns type: int - 0 if good

1029 

1030  Parameters: (in calling order)

1031  Type Name I/O Description

1032  ---- ---- --- -----------

1033  float * prof I pointer to the start of the

1034  met data

1035  unsigned char * ubdat I unsigned char data

1036  int32_t ipx I pixel to look at

1037  int32_t iln I line to look at

1038  int32_t npix I # pixels

1039  int32_t nlin I # lines

1040  Modification history:

1041  Programmer Date Description of change

1042  ---------- ---- ---------------------

1043  W. Robinson 28-Jan-2019 Original development

1044 

1045 *******************************************************************/

1046  {

1047  int nlvl_cmp = 101;

1048  int32_t off, ilvl;

1049  /* just print each profile */

1050  printf( "mixr profile\n" );

1051  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

1052  {

1053  printf( "%f ", prof[ ipx + npix * ( iln + nlin * ilvl ) ] );

1054  }

1055  printf( "\nT profile\n" );

1056  off = npix * nlin * nlvl_cmp;

1057  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

1058  {

1059  printf( "%f ", prof[ off + ipx + npix * ( iln + nlin * ilvl ) ] );

1060  }

1061  printf( "\nPRESS profile\n" );

1062  off += npix * nlin * nlvl_cmp;

1063  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

1064  {

1065  printf( "%f ", prof[ off + ipx + npix * ( iln + nlin * ilvl ) ] );

1066  }

1067  printf( "\nHEIGHT profile\n" );

1068  off += npix * nlin * nlvl_cmp;

1069  for( ilvl = 0; ilvl < nlvl_cmp; ilvl++ )

1070  {

1071  printf( "%f ", prof[ off + ipx + npix * ( iln + nlin * ilvl ) ] );

1072  }

1073  printf( "\n" );

1074  /* on to the single-level params */

1075  off += npix * nlin * nlvl_cmp;

1076  printf( "TSFC: %f, ", prof[ off + ipx + npix * iln] );

1077  off += npix * nlin;

1078  printf( "PSFC: %f, ", prof[ off + ipx + npix * iln] );

1079  off += npix * nlin;

1080  printf( "WIND SP: %f, ", prof[ off + ipx + npix * iln] );

1081  off += npix * nlin;

1082  printf( "O3: %f\n", prof[ off + ipx + npix * iln] );

1083 

1084  off += npix * nlin;

1085  printf( "ICE FR: %f, ", prof[ off + ipx + npix * iln] );

1086  off += npix * nlin;

1087  printf( "SNO FR: %f, ", prof[ off + ipx + npix * iln] );

1088  off += npix * nlin;

1089  printf( "ALT O3: %f, ", prof[ off + ipx + npix * iln] );

1090  off += npix * nlin;

1091  printf( "ALT ICEC: %f\n", prof[ off + ipx + npix * iln] );

1092  /* the first 2 met byte quantities - the 101 level profile surface

1093  and tropopause levels */

1094  off = 0;

1095  printf( "SFC_LVL: %d, ", ubdat[ off + ipx + npix * iln ] );

1096  off += npix * nlin;

1097  printf( "TROP_LVL: %d\n", ubdat[ off + ipx + npix * iln ] );

1098  /* the 'alternate snow/ice type (from NSIDC) */

1099  printf( "ALT ICE: %d\n", i32dat[ ipx + npix * iln ] );

1100  /* the remaining byte values are the cloud mask stuff */

1101  off = 2 * npix * nlin;

1102  printf( "---- Cloud Mask values ----\n" );

1103  printf( "CLD_DET: %d, ", ubdat[ off + ipx + npix * iln ] );

1104  off += npix * nlin;

1105  printf( "CONF_CLD: %d, ", ubdat[ off + ipx + npix * iln ] );

1106  off += npix * nlin;

1107  printf( "CLR_66: %d, ", ubdat[ off + ipx + npix * iln ] );

1108  off += npix * nlin;

1109  printf( "CLR_95: %d\n", ubdat[ off + ipx + npix * iln ] );

1110  off += npix * nlin;

1111  printf( "CLR_99: %d, ", ubdat[ off + ipx + npix * iln ] );

1112  off += npix * nlin;

1113  printf( "SNO_SFC: %d, ", ubdat[ off + ipx + npix * iln ] );

1114  off += npix * nlin;

1115  printf( "WTR_SFC: %d, ", ubdat[ off + ipx + npix * iln ] );

1116  off += npix * nlin;

1117  printf( "COAST_SFC: %d\n", ubdat[ off + ipx + npix * iln ] );

1118  off += npix * nlin;

1119  printf( "DESERT_SFC: %d, ", ubdat[ off + ipx + npix * iln ] );

1120  off += npix * nlin;

1121  printf( "LND_SFC: %d, ", ubdat[ off + ipx + npix * iln ] );

1122  off += npix * nlin;

1123  printf( "NIGHT: %d, ", ubdat[ off + ipx + npix * iln ] );

1124  off += npix * nlin;

1125  printf( "GLINT: %d\n", ubdat[ off + ipx + npix * iln ] );

1126  off += npix * nlin;

1127  printf( "OCEAN_NO_SNOW: %d, ", ubdat[ off + ipx + npix * iln ] );

1128  off += npix * nlin;

1129  printf( "OCEAN_SNOW: %d, ", ubdat[ off + ipx + npix * iln ] );

1130  off += npix * nlin;

1131  printf( "LND_NO_SNOW: %d, ", ubdat[ off + ipx + npix * iln ] );

1132  off += npix * nlin;

1133  printf( "LND_SNOW: %d\n", ubdat[ off + ipx + npix * iln ] );

1134  off += npix * nlin;

1135  printf( "TST_H_138: %d, ", ubdat[ off + ipx + npix * iln ] );

1136  off += npix * nlin;

1137  printf( "TST_VIS_REFL: %d, ", ubdat[ off + ipx + npix * iln ] );

1138  off += npix * nlin;

1139  printf( "TST_VIS_RATIO: %d, ", ubdat[ off + ipx + npix * iln ] );

1140  off += npix * nlin;

1141  printf( "VIS_CLD_250: %d\n", ubdat[ off + ipx + npix * iln ] );

1142  off += npix * nlin;

1143  printf( "APPL_HCLD_138: %d, ", ubdat[ off + ipx + npix * iln ] );

1144  off += npix * nlin;

1145  printf( "APPL_VIS_REFL: %d, ", ubdat[ off + ipx + npix * iln ] );

1146  off += npix * nlin;

1147  printf( "APPL_VIS_NIR_RATIO: %d\n", ubdat[ off + ipx + npix * iln ] );

1148  off += npix * nlin;

1149  /* end? */

1150  return 0;

1151  }

1152 /* WDR **** test end-around to mk_cmp_prof */

1153 int mk_cmp_prof_dum( float *merra_p, float *merra_t, float *merra_q,

1154  float *merra_h, float sfcp, float sfct, float sfcrh, float lat,

1155  unsigned char *sfc_lvl, unsigned char *trop_lvl, double *cmp_p,

1156  double *cmp_t, double *cmp_mixr, double *cmp_h )

1157  {

1158 int i;

1159 double cmp_t_sav[] = {169.25176566420146, 197.35942572651811, 229.18349451149692, 256.22238210609601, 236.47864107966123, 244.04297143175555, 250.64605396062848, 255.49783596464286, 259.17682096634508, 260.50430632229126, 257.86788008325004, 255.09753268656993, 251.71866532458688, 246.82264824496451, 241.36301970269031, 235.5349625222656, 231.00564672010188, 230.31934481415178, 230.00717490122784, 230.93984650147675, 231.87556038742846, 230.052070835427, 226.94875351064999, 223.97554681909486, 221.12246382681033, 218.38062746953614, 216.05483908715939, 214.30433557126059, 212.6156777870857, 210.98488863951465, 209.57034993410883, 208.26534656680437, 207.00140487891071, 205.76477951859837, 204.56229701519121, 202.84595260841155, 200.1528837094385, 197.5345872033561, 194.98734605240188, 193.13906753041528, 192.9482175288388, 192.76222219178675, 192.58086243374066, 192.4039337826197, 193.53690804615815, 196.34328186967463, 199.08486131023724, 201.76426680177968, 204.38396142841776, 206.94626326403974, 209.4240557622999, 211.67158053966773, 213.87221009140157, 216.02762928535265, 218.13943261959307, 220.23999924343335, 222.60373354556071, 224.92166874653532, 227.195309078017, 229.42608466977327, 232.27719286025973, 235.2490670940424, 238.1670467211311, 241.03276062011722, 243.91780941249539, 246.75240714620259, 249.53800835678237, 251.86502860055904, 253.8587763082154, 255.81919333430849, 257.63971457966039, 259.29891772708743, 260.9312384604865, 262.72494342240276, 264.64753610551918, 266.53988144929855, 268.11359331445857, 269.60080666913444, 271.01286592906087, 272.25331134687139, 273.4751007120372, 274.88788770642338, 276.32402654921054, 277.67617255366838, 278.94743905007527, 280.29684643297759, 281.87085668556011, 283.47072360439103, 285.24727818913414, 286.9975584364218, 288.51282458375113, 289.86087429979796, 290.9986136887851, 292.01521247752953, 292.98262145756013, 294.3041990659554, 296.79569609374886, 298.98862037321277, 299.76233992810609, 300.50664057118684, 301.24118800223761};

1160 double cmp_mixr_sav[] = {0.0030000000000000001, 0.0030000000000000001, 0.0030000000000000001, 0.0030000000000000001, 0.0039293894442529974, 0.0041105013907719468, 0.0042181900744832808, 0.0042143496342834938, 0.0041817655664085482, 0.0040206101090032517, 0.0038733849643674614, 0.0037377625227534746, 0.0036014648248259636, 0.0034516475997961819, 0.0033167026788735537, 0.0031947714963073749, 0.0030799596646610928, 0.0030107622372746654, 0.0029488582454394335, 0.0028982059410255015, 0.0028502657213637195, 0.0027944201318237261, 0.002735964452701974, 0.0026799596042977963, 0.0026262174675935851, 0.0025745708296473057, 0.0025405554258138494, 0.0025323804296621087, 0.0025244942581909551, 0.0025168783380373483, 0.0025258323028017112, 0.0025409217199495382, 0.0025555363514953174, 0.0025643709192777427, 0.0025716416839911219, 0.0025546327064356832, 0.0024931393172969081, 0.002433353276920958, 0.0023751897106692147, 0.0023510380165387293, 0.0024102799348228623, 0.0024680149123855986, 0.0025243109557570816, 0.0025792315352972882, 0.0028074932565275365, 0.0032578001900602956, 0.0036977103048750353, 0.004127644072243016, 0.004547996715297346, 0.0049591401890269129, 0.0066045793017774166, 0.015756209535939904, 0.024716888476006896, 0.033493476685974542, 0.042092466755785318, 0.053504463790107065, 0.094126872282568774, 0.13396219155187214, 0.17303627285925954, 0.21137369397371381, 0.2419968812989492, 0.27022626866094007, 0.2979437192418275, 0.32516470551490928, 0.47386770138225009, 0.6199703171142984, 0.76354752329508391, 1.0340410272474989, 1.3924930017834096, 1.7449525187114237, 1.9710464191297217, 2.0459308727278538, 2.1196020481996749, 2.5849697692954705, 3.3725913913975889, 4.1478216949005233, 4.7278626505249006, 5.2595586368380491, 5.6780053593354278, 5.7887852404941018, 5.897899014248777, 6.0126629733222074, 6.1272686176990323, 6.2612715902934664, 6.4138250668383652, 6.5300489624532512, 6.6239450570569449, 7.0237148492647679, 7.815759173341303, 8.9338062006928265, 10.017175527482264, 10.798574208009653, 11.482688439271284, 12.855657707404751, 14.06097496976534, 15.451677990360071, 8.9087798522868766, 0.078892832367705168, 0.0013115547243546515, 0.0021613451297737749, 0.0030000000000000001};

1161 double cmp_h_sav[] = {80980.567272357701, 75343.876190403083, 70225.982411112418, 65480.628319062816, 61428.924238065607, 58024.80202191414, 54950.191567802227, 52150.37234583289, 49587.069685576469, 47234.584895480606, 45084.856633888245, 43122.927977172032, 41324.952076504, 39676.318479632821, 38164.924587601548, 36777.388481093178, 35497.403939587413, 34300.298902152295, 33167.424793334678, 32088.91680780418, 31057.094888806638, 30073.813431125316, 29143.272510847557, 28263.483356793786, 27430.050001579777, 26639.098949354469, 25886.659507298184, 25168.400851046954, 24481.003799506794, 23822.279970437601, 23190.014293749664, 22582.119678830768, 21996.942353189086, 21433.083441024752, 20889.269993475449, 20365.028820377422, 19861.183794284312, 19377.767922145697, 18913.565376345625, 18466.731479797825, 18033.779798131261, 17612.242600176087, 17201.593473696255, 16801.341547642598, 16409.703419308389, 16023.364531449324, 15640.568137969527, 15261.317752839021, 14885.61129854926, 14513.441856962587, 14144.824125415456, 13779.920194073944, 13418.849581909997, 13061.570770637654, 12708.041363820343, 12358.193358432209, 12011.718940291992, 11668.35429065583, 11328.08453068694, 10990.892967344218, 10656.282389846725, 10323.666245846885, 9992.9598216377744, 9664.201332244118, 9337.3653911268993, 9012.4264031402581, 8689.4162094312233, 8368.613449980001, 8050.4524006411266, 7735.0947019042451, 7422.5961680835417, 7113.1000037295735, 6806.6644285353668, 6503.107856676912, 6202.1379607025819, 5903.6184955729223, 5607.7124190540144, 5314.6091843005188, 5024.3518657619488, 4737.0487534217273, 4452.7644388850267, 4171.3480686550502, 3892.6336221371671, 3616.5978083090913, 3343.2676087364471, 3072.5944718766254, 2804.3890955354977, 2538.46734163782, 2274.6417716957349, 2012.7539832113782, 1752.8699864080911, 1495.1635489650282, 1239.7954748266484, 986.82637292768879, 736.25029352497961, 487.89682088402327, 241.65371459693083, 0, -243.21093547707787, -481.72170087446557, 0};

1162 *sfc_lvl = 97;

1163 *trop_lvl = 44;

1164 for( i = 0; i < 101; i++ )

1165  {

1166  *(cmp_t + i) = *(cmp_t_sav + i);

1167  *(cmp_mixr + i) = *(cmp_mixr_sav + i);

1168  *(cmp_h + i) = *(cmp_h_sav + i);

1169  }

1170  return 0;

1171  }

1172 

1173 int mk_cmp_prof( float *merra_p, float *merra_t, float *merra_q,

1174  float *merra_h, float sfcp, float sfct, float sfcrh, float lat,

1175  unsigned char *sfc_lvl, unsigned char *trop_lvl, double *cmp_p,

1176  double *cmp_t, double *cmp_mixr, double *cmp_h )

1177  {

1178  unsigned char get_trop( double *, double *, unsigned char );

1179  int height_profile_( double *, double *, double *, double *, int *,

1180  double * );

1181  int profile_to_101_( double *, double *, double *, int *, double *,

1182  double *, double *, double *, int * );

1183  int32_t nlvl = 42, nlvl_cmp = 101, ioz = 0, ilvl, kstart;

1184  int32_t nlvl_d = nlvl + 1;

1185  float oval, a_factor;

1186  double p_init[nlvl_d], lcl_p[nlvl_d], lcl_t[nlvl_d];

1187  double lcl_w[nlvl_d], lcl_h[nlvl_d], lat_d, sfcp_d;

1188  /*

1189  * Some l1 conventions that need to be addressed:

1190  * the p, q, mixr, h arrays need a extra bottom level at 1100 mb (beyond

1191  * normal seen p) and need to be double.

1192  * the merra_q needs to be converted to mixr via:

1193  * w = q ( 1 - q )^-1 q = spec hum, w = mixing ratio

1194  * and the merra_t needs 273.15 added to it to get deg K

1195  * and the sfc RH needs to be made into a mixing ratio

1196  * lastly, for levels with (-) t, ppropagate higher level to it

1197  */

1198  lat_d = lat;

1199  sfcp_d = sfcp;

1200 

1201  for( ilvl = 0; ilvl < nlvl; ilvl++ )

1202  {

1203  merra_q[ilvl] = 1000. * merra_q[ilvl] / ( 1. - merra_q[ilvl]);

1204  merra_t[ilvl] += 273.15;

1205  }

1206  sfct += 273.15;

1207  /* use cvt_rh_to_q */

1208  met_cvt_rh_to_q( 1, &sfcp, MET_UNITS__P_HPA, &sfct, MET_UNITS__T_K,

1209  &sfcrh, &oval, MET_UNITS__Q_G_KG );

1210 /* WDR test to see that the conversion worked

1211 float ex_val;

1212 met_cvt_q_to_rh(1, &sfcp, MET_UNITS__P_HPA, &sfct, MET_UNITS__T_K,

1213  &oval, MET_UNITS__Q_KG_KG, &ex_val );

1214  END TEST */

1215  sfcrh = oval / ( 1. - oval ); /* now, this is the mixing ratio */

1216 

1217  /* transfer the 42 level profile info to the local double 43 level arrays */

1218  for( ilvl = 0; ilvl < nlvl; ilvl++ )

1219  {

1220  lcl_p[ilvl] = merra_p[ilvl];

1221  lcl_t[ilvl] = merra_t[ilvl];

1222  lcl_w[ilvl] = merra_q[ilvl];

1223  lcl_h[ilvl] = merra_h[ilvl];

1224 

1225  if( lcl_t[ilvl] < 0. )

1226  {

1227  lcl_t[ilvl] = lcl_t[ ilvl - 1 ];

1228  lcl_w[ilvl] = lcl_w[ ilvl - 1 ];

1229  lcl_h[ilvl] = lcl_h[ ilvl - 1 ];

1230  }

1231  }

1232  lcl_p[nlvl] = 1100.;

1233  lcl_t[nlvl] = lcl_t[ nlvl - 1 ];

1234  lcl_w[nlvl] = lcl_w[ nlvl - 1 ];

1235  lcl_h[nlvl] = lcl_h[ nlvl - 1 ];

1236 

1237 /* we will just duplicate the method used in ancillary_module.f90

1238  to set up the incoming profile and placing the surface values

1239  and then interpolation even though I'd probably do it differently

1240 */

1241  memcpy( p_init, lcl_p, nlvl_d * sizeof(double) );

1242 

1243  /* first, if sfc p is > next lowest p level, set the p level as the last */

1244  if( ( sfcp > 0 ) && ( lcl_p[ nlvl_d - 2 ] > sfcp ) )

1245  *sfc_lvl = nlvl_d - 1;

1246  else

1247  *sfc_lvl = 0;

1248 

1249  /* look through the lower half of the levels and find the sfc_lvl for

1250  the rest */

1251  kstart = nlvl_d / 2;

1252 

1253  for( ilvl = kstart; ilvl < nlvl_d; ilvl++ )

1254  {

1255  if( ( sfcp > 0 ) && ( lcl_p[ilvl] > sfcp ) )

1256  {

1257  if( *sfc_lvl == 0 )

1258  {

1259  *sfc_lvl = ilvl;

1260  lcl_t[ilvl] = sfct;

1261  /* now this next part is ver-batum from chimaera code but I don't

1262  understand the reason they did it this way */

1263  if( ( sfcp - lcl_p[ ilvl - 1 ] < 5. ) ||

1264  ( lcl_p[ilvl] - sfcp < 5. ) )

1265  {

1266  lcl_p[ilvl] = ( lcl_p[ilvl] + lcl_p[ ilvl - 1 ] ) / 2.;

1267  }

1268  else

1269  {

1270  lcl_p[ilvl] = sfcp;

1271  }

1272  lcl_w[ilvl] = sfcrh;

1273  }

1274  else

1275  {

1276  lcl_t[ilvl] = sfct;

1277  lcl_w[ilvl] = sfcrh;

1278  lcl_p[ilvl] = p_init[ ilvl - 1 ];

1279  }

1280  }

1281  }

1282  /* more dubvious mods: add the surface level into the coarse profile */

1283  if( *sfc_lvl != nlvl_d )

1284  merra_p[ nlvl_d - 1 ] = p_init[ ilvl - 1 ];

1285  else

1286  lcl_p[ nlvl - 1 ] = sfcp;

1287  lcl_t[ nlvl - 1 ] = sfct;

1288  lcl_w[ nlvl - 1 ] = sfcrh;

1289 

1290  /* on to getting the lowest level of the 101 level profile */

1291  kstart = nlvl_cmp / 2;

1292  for( ilvl = kstart; ilvl < nlvl_cmp; ilvl++ )

1293  {

1294  if( cmp_p[ilvl] >= sfcp )

1295  {

1296  *sfc_lvl = ilvl;

1297  break;

1298  }

1299  }

1300 

1301  profile_to_101_( lcl_p, lcl_t, lcl_w, &nlvl, &lat_d, cmp_p, cmp_t,

1302  cmp_mixr, &ioz );

1303 

1304  /* instead of assigning the sfct, sfcrh (now mixr) to the sfc_lvl, it

1305  interpolates in p coordinates the resulting T, mixr */

1306  a_factor = ( sfcp - cmp_p[ *sfc_lvl - 1 ] ) /

1307  ( cmp_p[*sfc_lvl] - cmp_p[ *sfc_lvl - 1 ] );

1308  cmp_t[*sfc_lvl] = cmp_t[ *sfc_lvl - 1 ] + a_factor *

1309  ( cmp_t[*sfc_lvl] - cmp_t[ *sfc_lvl - 1 ] );

1310  cmp_mixr[*sfc_lvl] = cmp_mixr[ *sfc_lvl - 1 ] + a_factor *

1311  ( cmp_mixr[*sfc_lvl] - cmp_mixr[ *sfc_lvl - 1 ] );

1312 

1313  cmp_p[*sfc_lvl] = sfcp;

1314 

1315  /* and then follow with the height profile */

1316  height_profile_( cmp_p, cmp_t, cmp_mixr, cmp_h, &nlvl_cmp, &sfcp_d );

1317  cmp_h[nlvl_cmp - 1 ] = 0.;

1318 

1319  /* also find the tropopause level - call a chimaera-based algorithm */

1320  *trop_lvl = get_trop( cmp_t, cmp_p, *sfc_lvl );

1321 

1322  return 0;

1323  }

1324 

1325 unsigned char get_trop( double *cmp_t, double *cmp_p, unsigned char sfc_lvl )

1326  {

1327  unsigned char trop_lvl;

1328  int32_t ilev, imin;

1329  float xmin, ptop = 100.;

1330 

1331  xmin = 999999.;

1332  imin = 1;

1333 

1334  for( ilev = 0; ilev < sfc_lvl - 5; ilev++ )

1335  {

1336  if( ( cmp_t[ilev] < xmin ) && ( cmp_p[ilev] > ptop ) )

1337  {

1338  xmin = cmp_t[ilev];

1339  imin = ilev;

1340  }

1341  }

1342 

1343  /* don't allow trop height > 400 mb (level 71) */

1344  trop_lvl = 200;

1345  for( ilev = imin; ilev < 71; ilev++ )

1346  {

1347  if( ( cmp_t[ ilev - 1 ] >= cmp_t[ilev] ) &&

1348  ( cmp_t[ ilev + 1 ] > cmp_t[ilev] ) )

1349  {

1350  trop_lvl = ilev;

1351  break;

1352  }

1353  }

1354 

1355  if( trop_lvl == 200 ) trop_lvl = imin;

1356 

1357  return trop_lvl;

1358  }

1359 

1360 int ncio_grab_ub_ds(int ncid, char *ds_name, unsigned char *data)

1361 /*******************************************************************

1362 

1363  ncio_grab_ub_ds

1364 

1365  purpose: grab dataset and return it in unsigned char format

1366 

1367  Returns type: int - 0 if OK, -1 if problem

1368 

1369  Parameters: (in calling order)

1370  Type Name I/O Description

1371  ---- ---- --- -----------

1372  int ncid I netcdf id of file

1373  char * ds_name I name of dataset to read

1374  unsigned char * data O pointer to a pre-allocated

1375  array to receive the data

1376 

1377  Modification history:

1378  Programmer Date Description of change

1379  ---------- ---- ---------------------

1380  W. Robinson 26 Nov 2018 original development

1381 

1382  *******************************************************************/ {

1383  int status;

1384  int var_id;

1385  /*

1386  * get ID of the dataset

1387  */

1388  if ((status = nc_inq_varid(ncid, ds_name, &var_id)) != NC_NOERR) {

1389  printf("%s, %d: nc_inq_varid returned error %d\n", __FILE__, __LINE__,

1390  status);

1391  return -1;

1392  }

1393  /*

1394  * read the dataset as a unsigned char

1395  */

1396  if ((status = nc_get_var_uchar(ncid, var_id, data)) != NC_NOERR) {

1397  printf("%s, %d: nc_get_var_uchar returned error %d\n",

1398  __FILE__, __LINE__, status);

1399  return -1;

1400  }

1401  return 0;

1402 }

1403 

1404 int ncio_grab_i32_ds(int ncid, char *ds_name, int32_t *data)

1405 /*******************************************************************

1406 

1407  ncio_grab_i32_ds

1408 

1409  purpose: grab dataset and return it in nteger format

1410 

1411  Returns type: int - 0 if OK, -1 if problem

1412 

1413  Parameters: (in calling order)

1414  Type Name I/O Description

1415  ---- ---- --- -----------

1416  int ncid I netcdf id of file

1417  char * ds_name I name of dataset to read

1418  int32_t * data O pointer to a pre-allocated

1419  array to receive the data

1420 

1421  Modification history:

1422  Programmer Date Description of change

1423  ---------- ---- ---------------------

1424  W. Robinson 26 Nov 2018 original development

1425 

1426  *******************************************************************/ {

1427  int status;

1428  int var_id;

1429  /*

1430  * get ID of the dataset

1431  */

1432  if ((status = nc_inq_varid(ncid, ds_name, &var_id)) != NC_NOERR) {

1433  printf("%s, %d: nc_inq_varid returned error %d\n", __FILE__, __LINE__,

1434  status);

1435  return -1;

1436  }

1437  /*

1438  * read the dataset as a unsigned char

1439  */

1440  if ((status = nc_get_var_int(ncid, var_id, (int32_t *)data)) != NC_NOERR) {

1441  printf("%s, %d: nc_get_var_int returned error %d\n",

1442  __FILE__, __LINE__, status);

1443  return -1;

1444  }

1445  return 0;

1446 }