38 #include "ambralsfor.h"
206 #define TWOSTR_NERR 22
242 #define A(i,j) a[i-1+(j-1)*lda]
243 #define AA(j,k) aa[j-1+(k-1)*ia]
244 #define ABD(i,j) abd[i-1+(j-1)*lda]
245 #define ALBMED(iu) out->albmed[iu-1]
246 #define AMB(iq,jq) ab[iq-1+(jq-1)*(ds->nstr/2)].zero
247 #define APB(iq,jq) ab[iq-1+(jq-1)*(ds->nstr/2)].one
249 #define B(iq) b[iq-1]
250 #define BDR(iq,jq) bdr[iq-1+(jq)*(ds->nstr/2)]
251 #define BEM(iq) bem[iq-1]
253 #define CBAND(irow,ncol) cband[irow-1+(ncol-1)*(9*(ds->nstr/2)-2)]
254 #define CC(iq,jq) cc[iq-1+(jq-1)*ds->nstr]
255 #define CH(lc) ch[lc-1]
256 #define CHTAU(ls) chtau[ls]
257 #define CMU(iq) cmu[iq-1]
258 #define CWT(iq) cwt[iq-1]
260 #define DFDT(lu) out->rad[lu-1].dfdt
261 #define DIAG(i) diag[i-1].on
262 #define DTAUC(lc) ds->dtauc[lc-1]
263 #define DTAU_C(lc) dtau_c[lc-1]
264 #define DTAUCPR(lc) dtaucpr[lc-1]
266 #define EMU(iu) emu[iu-1]
267 #define EVAL(j) eval[j-1]
268 #define EVEC(j,k) evec[j-1+(k-1)*ievec]
269 #define EVECC(iq,jq) evecc[iq-1+(jq-1)*ds->nstr]
270 #define EXPBEA(lc) expbea[lc]
272 #define FLDIR(lu) fl[lu-1].zero
273 #define FLDN(lu) fl[lu-1].one
274 #define FLUP(lu) out->rad[lu-1].flup
275 #define FLYR(lc) flyr[lc-1]
277 #define GC(iq,jq,lc) gc[iq-1+(jq-1+(lc-1)*ds->nstr)*ds->nstr]
278 #define GENSRC(maz,lc,iq) ds->gensrc[iq-1+(lc-1+maz*ds->nlyr)*ds->nstr]
279 #define GENSRCU(maz,lc,iu) ds->gensrcu[iu-1+(lc-1+maz*ds->nlyr)*ds->numu]
280 #define GG(lc) gg[lc-1]
281 #define GGPRIM(lc) ggprim[lc-1]
282 #define GL(k,lc) gl[k+(lc-1)*(ds->nstr+1)]
283 #define GMU(k) gmu[k-1]
284 #define GU(iu,iq,lc) gu[iu-1+(iq-1+(lc-1)*ds->nstr)*ds->numu]
285 #define GWT(k) gwt[k-1]
287 #define IERROR(i) ierror[i-1]
288 #define IPVT(k) ipvt[k-1]
290 #define KK(iq,lc) kk[iq-1+(lc-1)*ds->nstr]
292 #define LAYRU(lu) layru[lu-1]
293 #define LL(iq,lc) ll[iq-1+(lc-1)*ds->nstr]
295 #define MU(i) mu[i-1]
297 #define OMEGA(lyr) omega[lyr-1]
298 #define OPRIM(lc) oprim[lc-1]
300 #define PKAG(lc) pkag[lc]
301 #define PKAGC(lc) pkagc[lc-1]
302 #define PHASA(lc) phasa[lc-1]
303 #define PHASE(lc) phase[lc-1]
304 #define PHASM(lc) phasm[lc-1]
305 #define PHAST(lc) phast[lc-1]
306 #define PHI(j) ds->phi[j-1]
307 #define PHIRAD(jp) phirad[jp-1]
308 #define PMOM(k,lc) ds->pmom[k+(lc-1)*(ds->nmom_nstr+1)]
309 #define PRNTU0(i) prntu0[i-1]
310 #define PSI0(iq) psi[iq-1].zero
311 #define PSI1(iq) psi[iq-1].one
313 #define RFLDIR(lu) out->rad[lu-1].rfldir
314 #define RFLDN(lu) out->rad[lu-1].rfldn
315 #define RMU(iu,iq) rmu[iu-1+(iq)*ds->numu]
316 #define RR(lc) rr[lc-1]
318 #define SSALB(lc) ds->ssalb[lc-1]
319 #define SUBD(i) diag[i-1].sub
320 #define SUPERD(i) diag[i-1].super
321 #define SX(i) sx[i-1]
322 #define SY(i) sy[i-1]
324 #define TAU(lc) tau[lc]
325 #define TAUC(lc) tauc[lc]
326 #define TAUCPR(lc) taucpr[lc]
327 #define TEMPER(lc) ds->temper[lc]
328 #define TRNMED(iu) out->trnmed[iu-1]
330 #define U0C(iq,lu) u0c[iq-1+(lu-1)*ds->nstr]
331 #define U0U(iu,lu) out->u0u[iu-1+(lu-1)*ds->numu]
332 #define UAVG(lu) out->rad[lu-1].uavg
333 #define UAVGDN(lu) out->rad[lu-1].uavgdn
334 #define UAVGUP(lu) out->rad[lu-1].uavgup
335 #define UAVGSO(lu) out->rad[lu-1].uavgso
336 #define UMU(iu) ds->umu[iu-1]
337 #define UTAU(lu) ds->utau[lu-1]
338 #define UTAUPR(lu) utaupr[lu-1]
339 #define UUM(iu,lu) uum[iu-1+(lu-1)*ds->numu]
340 #define UU(iu,lu,j) out->uu[iu-1+(lu-1+(j-1)*ds->ntau)*ds->numu]
341 #define OUT_UUM(iu,lu,j) out->uum[iu-1+(lu-1+(j)*ds->ntau)*ds->numu]
343 #define WK(iq) wk[iq-1]
345 #define XBA(lc) xba[lc]
346 #define XB0(iq,lc) xb[iq-1+(lc-1)*ds->nstr].zero
347 #define XB1(iq,lc) xb[iq-1+(lc-1)*ds->nstr].one
348 #define XB_0D(lc) ts[lc-1].xb_0d
349 #define XB_0U(lc) ts[lc-1].xb_0u
350 #define XB_1D(lc) ts[lc-1].xb_1d
351 #define XB_1U(lc) ts[lc-1].xb_1u
352 #define XP_0(lc) ts[lc-1].xp_0
353 #define XP_1(lc) ts[lc-1].xp_1
354 #define XR0(lc) xr[lc-1].zero
355 #define XR1(lc) xr[lc-1].one
357 #define YB_0D(lc) ts[lc-1].yb_0d
358 #define YB_0U(lc) ts[lc-1].yb_0u
359 #define YB_1D(lc) ts[lc-1].yb_1d
360 #define YB_1U(lc) ts[lc-1].yb_1u
361 #define YLM(l,i) ylm[l+(i-1)*(maxmu+1)]
362 #define YLM0(iq) ylm0[iq]
363 #define YLMC(l,iq) ylmc[l+(iq-1)*(ds->nstr+1)]
364 #define YLMU(l,iu) ylmu[l+(iu-1)*(ds->nstr+1)]
365 #define YP_0D(lc) ts[lc-1].yp_0d
366 #define YP_0U(lc) ts[lc-1].yp_0u
367 #define YP_1D(lc) ts[lc-1].yp_1d
368 #define YP_1U(lc) ts[lc-1].yp_1u
371 #define Z0(iu) zee[iu-1].zero
372 #define Z1(iq) zee[iq-1].one
373 #define Z0U(iu,lc) zu[iu-1+(lc-1)*ds->numu].zero
374 #define Z1U(iu,lc) zu[iu-1+(lc-1)*ds->numu].one
375 #define ZB0U(iu,lc) zbu[iu-1+(lc-1)*ds->numu].zero
376 #define ZB1U(iu,lc) zbu[iu-1+(lc-1)*ds->numu].one
377 #define ZBAU(iu,lc) zbu[iu-1+(lc-1)*ds->numu].alpha
378 #define ZB_A(lc) ts[lc-1].zb_a
379 #define ZBEAM(iu,lc) zbeam[iu-1+(lc-1)*ds->numu]
380 #define ZBEAMA(lc) zbeama[lc-1]
381 #define ZBEAM0(iq,lc) zbeamsp[iq-1+(lc-1)*ds->nstr].zero
382 #define ZBEAM1(iq,lc) zbeamsp[iq-1+(lc-1)*ds->nstr].one
383 #define ZBS0(iq) zbs[iq-1].zero
384 #define ZBS1(iq) zbs[iq-1].one
386 #define ZJ(j) zj[j-1]
387 #define ZJG(j) zjg[j-1]
388 #define ZJU(j) zju[j-1]
389 #define ZGU(iu,lc) zgu[iu-1+(lc-1)*ds->numu]
390 #define ZP_A(lc) ts[lc-1].zp_a
391 #define ZPLK0(iq,lc) plk[iq-1+(lc-1)*ds->nstr].zero
392 #define ZPLK1(iq,lc) plk[iq-1+(lc-1)*ds->nstr].one
393 #define ZZ(iq,lc) zz[iq-1+(lc-1)*ds->nstr]
394 #define ZZG(iq,lc) zzg[iq-1+(lc-1)*ds->nstr]
397 #define MUP(it) mu_phase[it-1]
398 #define PHASR(lc) phasr[lc-1]
399 #define PHAS2(it,lc) phas2[it-1+(lc-1)*nphase]
400 #define DSPHASE(it,lc) ds->phase[it-1+(lc-1)*ds->nphase]
401 #define F_PHAS2_ABS(it) f_phas2_abs[it-1]
402 #define MU_EQ(i,lu) mu_eq[i-1+(lu-1)*nf]
403 #define NEG_PHAS(i,lu) neg_phas[i-1+(lu-1)*nf]
404 #define NORM_PHAS(lu) norm_phas[lu-1]
407 #define SDTAUC(i) sdtauc[i-1]
408 #define SUTAU(i) sutau[i-1]
409 #define ZOUT(i) zout[i-1]
410 #define TAUINT(i) tauint[i-1]
411 #define XARR(i) xarr[i-1]
412 #define YARR(i) yarr[i-1]
420 #define FIRST_IPHAS 1
423 #define HENYEY_GREENSTEIN 3
424 #define HAZE_GARCIA_SIEWERT 4
425 #define CLOUD_GARCIA_SIEWERT 5
447 #define BRDF_CAM_NN 4
448 #define BRDF_CAM_SAL 0
449 #define BRDF_CAM_PCL 1
450 #define BRDF_CAM_U10 2
451 #define BRDF_CAM_UPHI 3
465 # define M_E 2.7182818284590452354
466 # define M_LOG2E 1.4426950408889634074
467 # define M_LOG10E 0.43429448190325182765
468 # define M_LN2 0.69314718055994530942
469 # define M_LN10 2.30258509299404568402
470 # define M_PI 3.14159265358979323846
471 # define M_PI_2 1.57079632679489661923
472 # define M_PI_4 0.78539816339744830962
473 # define M_1_PI 0.31830988618379067154
474 # define M_2_PI 0.63661977236758134308
475 # define M_2_SQRTPI 1.12837916709551257390
476 # define M_SQRT2 1.41421356237309504880
477 # define M_SQRT1_2 0.70710678118654752440
480 #define DEG (M_PI/180.)
483 const double _x = (double)(x); \
486 #define MIN(x,y) ({ \
487 const double _x = (double)(x); \
488 const double _y = (double)(y); \
489 _x < _y ? _x : _y; })
491 #define MAX(x,y) ({ \
492 const double _x = (double)(x); \
493 const double _y = (double)(y); \
494 _x > _y ? _x : _y; })
496 #define LIMIT_RANGE(min,x,max) ({ \
497 const double _min = (double)(min); \
498 const double _x = (double)(x); \
499 const double _max = (double)(max); \
500 _x < _min ? _min : ( _x > _max ? _max : _x ); })
502 #define IMIN(i,j) ({ \
503 const int _i = (int)(i); \
504 const int _j = (int)(j); \
505 _i < _j ? _i : _j; })
507 #define IMAX(i,j) ({ \
508 const int _i = (int)(i); \
509 const int _j = (int)(j); \
510 _i > _j ? _i : _j; })
512 #define F77_SIGN(a,b) ((b) >= 0. ? fabs(a) : -fabs(a))
542 double c_dref(
double wvnmlo,
1042 char const *varnam);
void c_albtrans(disort_state *ds, disort_output *out, disort_pair *ab, double *array, double *b, double *bdr, double *cband, double *cc, double *cmu, double *cwt, double *dtaucpr, double *eval, double *evecc, double *gl, double *gc, double *gu, int *ipvt, double *kk, double *ll, int nn, double *taucpr, double *ylmc, double *ylmu, double *z, double *wk)
void c_legendre_poly(int nmu, int m, int maxmu, int twonm1, double *mu, double *ylm)
void c_errmsg(char const *messag, int type)
void c_print_albtrans(disort_state *ds, disort_output *out)
double c_planck_func1(double wnumlo, double wnumhi, double t)
double c_xi_func(double umu1, double umu2, double tau)
int * c_int_vector(int nl, int nh, char const *name)
void prep_double_scat_integr(int nphase, int ntau, int nf, double *mu_phase, double *phas2, double *mu_eq, int *neg_phas, double *norm_phas)
void c_twostr_check_inputs(disort_state *ds, double *gg, int *ierror, double *tauc)
void c_fluxes(disort_state *ds, disort_output *out, double *ch, double *cmu, double *cwt, double *gc, double *kk, int *layru, double *ll, int lyrcut, int ncut, int nn, int prntu0, double *taucpr, double *utaupr, disort_pair *xr, disort_pair *zbeamsp, double *zbeama, double *zz, double *zzg, disort_pair *plk, disort_pair *fl, double *u0c)
void c_disort_out_free(disort_state *ds, disort_output *out)
void c_interp_coefficients_beam_source(disort_state *ds, double *chtau, double delm0, double fbeam, double *gl, int lc, int mazim, int nstr, int numu, double *taucpr, disort_triplet *zbu, double *xba, double *zj, double *ylm0, double *ylmu)
void c_solve0(disort_state *ds, double *b, double *bdr, double *bem, double bplanck, double *cband, double *cmu, double *cwt, double *expbea, int *ipvt, double *ll, int lyrcut, int mazim, int ncol, int ncut, int nn, double tplanck, double *taucpr, double *z, disort_pair *zbeamsp, double *zbeama, double *zz, double *zzg, disort_pair *plk)
double c_ratio(double a, double b)
int c_gaussian_quadrature_test(int nstr, float *sza, double umu0)
double * c_dbl_vector(int nl, int nh, char const *name)
void c_twostr_state_alloc(disort_state *ds)
void c_twostr_fluxes(disort_state *ds, twostr_xyz *ts, double *ch, double cmu, double *kk, int *layru, double *ll, int lyrcut, int ncut, double *oprim, double *rr, double *taucpr, double *utaupr, disort_output *out, double *u0c, disort_pair *fl)
double c_planck_func2(double wnumlo, double wnumhi, double t)
void c_twostr_out_alloc(disort_state *ds, disort_output *out)
void c_twostr_set(disort_state *ds, double *bplanck, double *ch, double *chtau, double *cmu, int deltam, double *dtaucpr, double *expbea, double *flyr, double *gg, double *ggprim, int *layru, int *lyrcut, int *ncut, int *nn, double *oprim, double *pkag, double *pkagc, double radius, double *tauc, double *taucpr, double *tplanck, double *utaupr)
void c_user_intensities(disort_state *ds, double bplanck, double *cmu, double *cwt, double delm0, double *dtaucpr, double *emu, double *expbea, double *gc, double *gu, double *kk, int *layru, double *ll, int lyrcut, int mazim, int ncut, int nn, double *rmu, double *taucpr, double tplanck, double *utaupr, double *wk, disort_triplet *zbu, double *zbeam, disort_pair *zbeamsp, double *zbeama, double *zgu, disort_pair *zu, double *zz, double *zzg, disort_pair *plk, double *uum)
void c_upbeam_pseudo_spherical(disort_state *ds, int lc, double *array, double *cc, double *cmu, int *ipvt, int nn, double *wk, disort_pair *xb, double *xba, disort_pair *zbs, double *zbsa, disort_pair *zbeamsp, double *zbeama)
void c_sgesl(double *a, int lda, int n, int *ipvt, double *b, int job)
void c_sgbfa(double *abd, int lda, int n, int ml, int mu, int *ipvt, int *info)
void c_print_avg_intensities(disort_state *ds, disort_output *out)
void c_twostr_out_free(disort_state *ds, disort_output *out)
void c_new_intensity_correction(disort_state *ds, disort_output *out, double dither, double *flyr, int *layru, int lyrcut, int ncut, double *oprim, double *phasa, double *phast, double *phasm, double *phirad, double *tauc, double *taucpr, double *utaupr)
double c_secondary_scat(disort_state *ds, int iu, int lu, double ctheta, double *flyr, int layru, double *tauc)
void c_albtrans_intensity(disort_state *ds, disort_output *out, double *gu, double *kk, double *ll, int nn, double *taucpr, double *wk)
void c_check_inputs(disort_state *ds, int scat_yes, int deltam, int corint, double *tauc, int callnum)
void c_disort_out_alloc(disort_state *ds, disort_output *out)
void c_getmom(int iphas, double gg, int nmom, double *pmom)
void c_disort(disort_state *ds, disort_output *out)
void c_saxpy(int n, double sa, double *sx, double *sy)
void c_twostr(disort_state *ds, disort_output *out, int deltam, double *gg, int *ierror, double radius)
void c_set_coefficients_beam_source(disort_state *ds, double *ch, double *chtau, double *cmu, double delm0, double fbeam, double *gl, int lc, int mazim, int nstr, double *taucpr, double *xba, disort_pair *xb, double *ylm0, double *ylmc, double *zj)
int c_write_bad_var(int quiet, char const *varnam)
void c_sgbsl(double *abd, int lda, int n, int ml, int mu, int *ipvt, double *b, int job)
double c_chapman_simpler(int lc, double taup, int nlyr, double *zd, double *dtau_c, double zenang, double r)
int c_write_too_small_dim(int quiet, char const *dimnam, int minval)
void c_disort_state_alloc(disort_state *ds)
void c_solve_eigen(disort_state *ds, int lc, disort_pair *ab, double *array, double *cmu, double *cwt, double *gl, int mazim, int nn, double *ylmc, double *cc, double *evecc, double *eval, double *kk, double *gc, double *wk)
double c_bidir_reflectivity_rpv(rpv_brdf_spec *brdf, double mu1, double mu2, double phi, double badmu)
void c_upbeam(disort_state *ds, int lc, double *array, double *cc, double *cmu, double delm0, double *gl, int *ipvt, int mazim, int nn, double *wk, double *ylm0, double *ylmc, double *zj, double *zz)
void c_sgefa(double *a, int lda, int n, int *ipvt, int *info)
void c_twostr_solve_bc(disort_state *ds, twostr_xyz *ts, double bplanck, double *cband, double cmu, double *expbea, int lyrcut, int nn, int ncut, double tplanck, double *taucpr, double *kk, double *rr, int *ipvt, double *b, double *ll, twostr_diag *diag)
void c_twostr_print_inputs(disort_state *ds, int deltam, double *flyr, double *gg, int lyrcut, double *oprim, double *tauc, double *taucpr)
void c_free_dbl_vector(double *m, int nl, int nh)
double c_inter(int npoints, int itype, double arg, float *xarr, double *yarr, double *hh)
void c_self_test(int compare, int *prntu0, disort_state *ds, disort_output *out)
void c_intensity_correction(disort_state *ds, disort_output *out, double dither, double *flyr, int *layru, int lyrcut, int ncut, double *oprim, double *phasa, double *phast, double *phasm, double *phirad, double *tauc, double *taucpr, double *utaupr)
void c_asymmetric_matrix(double *aa, double *evec, double *eval, int m, int ia, int ievec, int *ier, double *wk)
void c_upisot(disort_state *ds, int lc, double *array, double *cc, double *cmu, int *ipvt, int nn, double *oprim, double *wk, disort_pair *xr, disort_pair *zee, disort_pair *plk)
double c_dref(double wvnmlo, double wvnmhi, double mu, int brdf_type, disort_brdf *brdf, int callnum)
void c_twostr_solns(disort_state *ds, double *ch, double *chtau, double cmu, int ncut, double *oprim, double *pkag, double *pkagc, double *taucpr, double *ggprim, double *kk, double *rr, twostr_xyz *ts)
void c_disort_state_free(disort_state *ds)
int c_setout(float *sdtauc, int nlyr, int ntau, float *sutau, float *z, float *zout)
int c_fcmp(double x1, double x2)
void c_gaussian_quadrature(int m, double *gmu, double *gwt)
double calc_phase_squared(int nphase, int lu, double ctheta, int nf, double *mu_phase, double *phas2, double *mu_eq, int *neg_phas, double norm_phas)
void c_set_matrix(disort_state *ds, double *bdr, double *cband, double *cmu, double *cwt, double delm0, double *dtaucpr, double *gc, double *kk, int lyrcut, int *ncol, int ncut, double *taucpr, double *wk)
void c_print_intensities(disort_state *ds, disort_output *out)
void c_intensity_components(disort_state *ds, double *gc, double *kk, int *layru, double *ll, int lyrcut, int mazim, int ncut, int nn, double *taucpr, double *utaupr, double *zz, disort_pair *plk, double *uum)
double c_sdot(int n, double *sx, double *sy)
double c_new_secondary_scat(disort_state *ds, int iu, int lu, int it, double ctheta, double *flyr, int layru, double *tauc, int nf, double *phas2, double *mu_eq, int *neg_phas, double norm_phas)
void c_print_inputs(disort_state *ds, double *dtaucpr, int scat_yes, int deltam, int corint, double *flyr, int lyrcut, double *oprim, double *tauc, double *taucpr)
void c_sscal(int n, double sa, double *sx)
void c_upbeam_general_source(disort_state *ds, int lc, int maz, double *array, double *cc, int *ipvt, int nn, double *wk, double *zjg, double *zzg)
double c_sasum(int n, double *sx)
struct disort_state disort_state
struct disort_output disort_output
void c_albtrans_spherical(disort_state *ds, double *cmu, double *cwt, double *gc, double *kk, double *ll, int nn, double *taucpr, double *sflup, double *sfldn)
double c_chapman(int lc, double taup, double *tauc, int nlyr, double *zd, double *dtau_c, double zenang, double r)
void c_interp_source(disort_state *ds, int lc, double *cwt, double delm0, double *gl, int mazim, double *oprim, double *ylm0, double *ylmc, double *ylmu, disort_pair *psi, disort_pair *xr, disort_pair *zee, double *zj, double *zjg, double *zbeam, disort_triplet *zbu, disort_pair *zbs, double zbsa, double *zgu, disort_pair *zu)
void c_interp_eigenvec(disort_state *ds, int lc, double *cwt, double *evecc, double *gl, double *gu, int mazim, int nn, double *wk, double *ylmc, double *ylmu)
int c_isamax(int n, double *sx)
double c_bidir_reflectivity(double wvnmlo, double wvnmhi, double mu, double mup, double dphi, int brdf_type, disort_brdf *brdf, int callnum)
double c_single_scat(double dither, int layru, int nlyr, double *phase, double *omega, double *tau, double umu, double umu0, double utau, double fbeam)
void c_twostr_state_free(disort_state *ds)
double c_bidir_reflectivity_hapke(double wvnmlo, double wvnmhi, double mu, double mup, double dphi)
void c_solve1(disort_state *ds, double *cband, int ihom, int *ipvt, int ncol, int ncut, int nn, double *b, double *ll)
void c_disort_set(disort_state *ds, double *ch, double *chtau, double *cmu, double *cwt, int deltam, double *dtaucpr, double *expbea, double *flyr, double *gl, int *layru, int *lyrcut, int *ncut, int *nn, int *corint, double *oprim, double *tauc, double *taucpr, double *utaupr)
void c_sgbco(double *abd, int lda, int n, int ml, int mu, int *ipvt, double *rcond, double *z)
void print_test(disort_state *ds_calc, disort_output *calc, disort_state *ds_good, disort_output *good)
void c_surface_bidir(disort_state *ds, double delm0, double *cmu, int mazim, int nn, double *bdr, double *emu, double *bem, double *rmu, int callnum)
void c_sgeco(double *a, int lda, int n, int *ipvt, double *rcond, double *z)
int old_intensity_correction