42#pragma GCC diagnostic push
43#pragma GCC diagnostic ignored "-Wunused-parameter"
53 const MatrixView& sca1,
54 const MatrixView& pfct2,
55 const MatrixView& sca2) {
56 const Index np1 = pfct1.npages();
57 const Index nr1 = pfct1.nrows();
58 const Index nc1 = pfct1.ncols();
65 for (Index i = 0; i < np1; i++) {
66 for (Index j = 0; j < nr1 ; j++) {
67 for (Index k = 0; k < nc1; k++)
70 (sca1(i, j) * pfct1(i, j, k) + sca2(i, j) * pfct2(j, k)) /
71 (sca1(i, j) + sca2(i, j));
77 const Index& cloudbox_on,
78 const Index& atmfields_checked,
79 const Index& atmgeom_checked,
80 const Index& cloudbox_checked,
81 const Index& scat_data_checked,
82 const Index& atmosphere_dim,
83 const Index& stokes_dim,
86 ConstVectorView za_grid,
87 const Index& nstreams) {
90 "Cloudbox is off, no scattering calculations to be"
94 if (atmfields_checked != 1)
96 "The atmospheric fields must be flagged to have "
97 "passed a consistency check (atmfields_checked=1).");
98 if (atmgeom_checked != 1)
100 "The atmospheric geometry must be flagged to have "
101 "passed a consistency check (atmgeom_checked=1).");
102 if (cloudbox_checked != 1)
104 "The cloudbox must be flagged to have "
105 "passed a consistency check (cloudbox_checked=1).");
106 if (scat_data_checked != 1)
108 "The scat_data must be flagged to have "
109 "passed a consistency check (scat_data_checked=1).");
111 if (atmosphere_dim != 1)
113 "For running DISORT, atmospheric dimensionality "
116 if (stokes_dim < 0 || stokes_dim > 1)
118 "For running DISORT, the dimension of stokes vector "
121 if (cloudbox_limits.
nelem() != 2 * atmosphere_dim)
123 "*cloudbox_limits* is a vector which contains the"
124 "upper and lower limit of the cloud for all "
125 "atmospheric dimensions. So its dimension must"
126 "be 2 x *atmosphere_dim*");
128 if (cloudbox_limits[0] != 0) {
130 os <<
"DISORT calculations currently only possible with "
131 <<
"lower cloudbox limit\n"
132 <<
"at 0th atmospheric level "
133 <<
"(assumes surface there, ignoring z_surface).\n";
134 throw runtime_error(os.str());
137 if (scat_data.empty())
139 "No single scattering data present.\n"
140 "See documentation of WSV *scat_data* for options to "
141 "make single scattering data available.\n");
148 if (nstreams / 2 * 2 != nstreams) {
150 os <<
"DISORT requires an even number of streams, but yours is " << nstreams
152 throw runtime_error(os.str());
156 Index nza = za_grid.nelem();
166 os <<
"We require size of za_grid to be >= 20, to ensure a\n"
167 <<
"reasonable interpolation of the calculated cloudbox field.\n"
168 <<
"Note that for DISORT additional computation costs for\n"
169 <<
"larger numbers of angles are negligible.";
170 throw runtime_error(os.str());
173 if (za_grid[0] != 0. || za_grid[nza - 1] != 180.)
174 throw runtime_error(
"The range of *za_grid* must [0 180].");
176 if (!is_increasing(za_grid))
177 throw runtime_error(
"*za_grid* must be increasing.");
180 while (za_grid[i] <= 90) {
181 if (za_grid[i] == 90)
182 throw runtime_error(
"*za_grid* is not allowed to contain the value 90");
187 bool all_totrand =
true;
188 for (Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
189 for (Index i_se = 0; i_se < scat_data[i_ss].
nelem(); i_se++)
190 if (scat_data[i_ss][i_se].ptype !=
PTYPE_TOTAL_RND) all_totrand =
false;
193 os <<
"DISORT can only handle scattering elements of type "
198 throw runtime_error(os.str());
203 const Index& atmfields_checked,
204 const Index& atmgeom_checked,
205 const Index& scat_data_checked,
206 const Index& atmosphere_dim,
207 const Index& stokes_dim,
209 const Index& nstreams) {
210 if (atmfields_checked != 1)
212 "The atmospheric fields must be flagged to have "
213 "passed a consistency check (atmfields_checked=1).");
215 if (atmgeom_checked != 1)
217 "The atmospheric geometry must be flagged to have "
218 "passed a consistency check (atmgeom_checked=1).");
220 if (scat_data_checked != 1)
222 "The scat_data must be flagged to have "
223 "passed a consistency check (scat_data_checked=1).");
225 if (atmosphere_dim != 1)
227 "For running DISORT, atmospheric dimensionality "
230 if (stokes_dim < 0 || stokes_dim > 1)
232 "For running DISORT, the dimension of stokes vector "
235 if (scat_data.empty())
237 "No single scattering data present.\n"
238 "See documentation of WSV *scat_data* for options to "
239 "make single scattering data available.\n");
246 if (nstreams / 2 * 2 != nstreams) {
248 os <<
"DISORT requires an even number of streams, but yours is " << nstreams
250 throw runtime_error(os.str());
254 bool all_totrand =
true;
255 for (Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
256 for (Index i_se = 0; i_se < scat_data[i_ss].
nelem(); i_se++)
257 if (scat_data[i_ss][i_se].ptype !=
PTYPE_TOTAL_RND) all_totrand =
false;
260 os <<
"DISORT can only handle scattering elements of type "
265 throw runtime_error(os.str());
270 Tensor7& cloudbox_field,
272 const Vector& f_grid,
276 const Index& stokes_dim) {
277 const Index Nf = f_grid.nelem();
278 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
282 cloudbox_field.resize(Nf, Np_cloud, 1, 1, n_za, n_aa, stokes_dim);
283 cloudbox_field = NAN;
290 ConstVectorView f_grid,
291 const Numeric& surface_skin_t,
292 ConstVectorView surface_scalar_reflectivity) {
294 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
296 os <<
"Surface temperature has been set or derived as " << btemp <<
" K,\n"
297 <<
"which is not considered a meaningful value.\n"
298 <<
"For surface method 'L', *surface_skin_t* needs to\n"
299 <<
"be set and passed explicitly. Maybe you didn't do this?";
300 throw runtime_error(os.str());
302 btemp = surface_skin_t;
305 if (surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
306 surface_scalar_reflectivity.nelem() != 1) {
308 os <<
"The number of elements in *surface_scalar_reflectivity*\n"
309 <<
"should match length of *f_grid* or be 1."
310 <<
"\n length of *f_grid* : " << f_grid.nelem()
311 <<
"\n length of *surface_scalar_reflectivity* : "
312 << surface_scalar_reflectivity.nelem() <<
"\n";
313 throw runtime_error(os.str());
316 if (
min(surface_scalar_reflectivity) < 0 ||
317 max(surface_scalar_reflectivity) > 1) {
319 "All values in *surface_scalar_reflectivity*"
320 " must be inside [0,1].");
323 if (surface_scalar_reflectivity.nelem() > 1)
324 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
325 albedo[f_index] = surface_scalar_reflectivity[f_index];
327 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
328 albedo[f_index] = surface_scalar_reflectivity[0];
332 MatrixView ext_bulk_gas,
333 const Agenda& propmat_clearsky_agenda,
334 ConstVectorView t_profile,
335 ConstMatrixView vmr_profiles,
336 ConstVectorView p_grid,
337 ConstVectorView f_grid) {
338 const Index Np = p_grid.nelem();
340 ARTS_ASSERT(ext_bulk_gas.nrows() == f_grid.nelem());
348 const Vector rtp_mag_dummy(3, 0);
349 const Vector ppath_los_dummy;
350 StokesVector nlte_dummy;
351 ArrayOfStokesVector partial_nlte_dummy;
352 ArrayOfPropagationMatrix partial_dummy;
354 PropagationMatrix propmat_clearsky_local;
355 for (Index ip = 0; ip < Np; ip++) {
357 propmat_clearsky_local,
369 Vector{vmr_profiles(joker, ip)},
370 propmat_clearsky_agenda);
371 ext_bulk_gas(joker, ip) += propmat_clearsky_local.Kjj();
376 MatrixView sca_coeff_gas,
377 MatrixView sca_coeff_gas_level,
379 const ConstVectorView& f_grid,
382 const MatrixView& vmr,
383 const Agenda& gas_scattering_agenda) {
384 const Index Np = p.nelem();
385 const Index Nl = pfct_gas.ncols();
386 const Index Nf = f_grid.nelem();
388 PropagationMatrix K_sca_gas_temp;
392 Matrix pmom_gas_level( Np, Nl, 0.);
396 for (Index ip = 0; ip < Np; ip++) {
404 Vector{vmr(joker, ip)},
408 gas_scattering_agenda);
411 sca_coeff_gas_level(joker, ip) = K_sca_gas_temp.Kjj(0, 0);
414 N_polys =
min(Nl, sca_fct_temp.nelem());
415 for (Index k = 0; k < N_polys; k++) {
416 pmom_gas_level( ip, k) = sca_fct_temp[k];
421 for (Index ip = 0; ip < Np - 1; ip++) {
422 for (Index f_index = 0; f_index < Nf; f_index++) {
424 sca_coeff_gas(f_index, Np - 2 - ip) =
426 (sca_coeff_gas_level(f_index, ip) +
427 sca_coeff_gas_level(f_index, ip + 1));
430 for (Index l_index = 0; l_index < Nl; l_index++) {
431 pfct_gas( Np - 2 - ip, l_index) =
432 0.5 * (pmom_gas_level( ip, l_index) +
433 pmom_gas_level( ip + 1, l_index));
439 MatrixView abs_bulk_par,
441 ConstMatrixView pnd_profiles,
442 ConstVectorView t_profile,
445 ConstVectorView f_grid) {
446 const Index Np_cloud = pnd_profiles.ncols();
448 const Index nf = f_grid.nelem();
460 Vector T_array{t_profile[Range(cloudbox_limits[0], Np_cloud)]};
461 Matrix dir_array(1, 2, 0.);
465 ArrayOfArrayOfTensor5 ext_mat_Nse;
466 ArrayOfArrayOfTensor4 abs_vec_Nse;
469 ArrayOfTensor5 ext_mat_ssbulk;
470 ArrayOfTensor4 abs_vec_ssbulk;
472 Tensor5 ext_mat_bulk;
473 Tensor4 abs_vec_bulk;
502 bool pf = (abs_vec_bulk.nbooks() != 1);
503 for (Index ip = 0; ip < Np_cloud; ip++)
504 for (Index f_index = 0; f_index < nf; f_index++) {
505 if (pf) f_this = f_index;
506 ext_bulk_par(f_index, ip + cloudbox_limits[0]) =
507 ext_mat_bulk(f_this, ip, 0, 0, 0);
508 abs_bulk_par(f_index, ip + cloudbox_limits[0]) =
509 abs_vec_bulk(f_this, ip, 0, 0);
515 ConstMatrixView ext_bulk_gas,
516 ConstMatrixView ext_bulk_par,
517 ConstMatrixView abs_bulk_par,
518 ConstVectorView z_profile) {
519 const Index nf = ext_bulk_gas.nrows();
520 const Index Np = ext_bulk_gas.ncols();
531 for (Index ip = 0; ip < Np - 1; ip++)
533 for (Index f_index = 0; f_index < nf; f_index++) {
535 0.5 * (ext_bulk_gas(f_index, ip) + ext_bulk_par(f_index, ip) +
536 ext_bulk_gas(f_index, ip + 1) + ext_bulk_par(f_index, ip + 1));
540 (ext_bulk_gas(f_index, ip) + abs_bulk_par(f_index, ip) +
541 ext_bulk_gas(f_index, ip + 1) + abs_bulk_par(f_index, ip + 1));
542 ssalb(f_index, Np - 2 - ip) = (ext - abs) / ext;
546 "pressure level = ", ip,
"\n",
547 "frequency number = ", f_index,
"\n");
551 dtauc(f_index, Np - 2 - ip) = ext * (z_profile[ip + 1] - z_profile[ip]);
557 const Index& Npfct) {
558 const Index min_nang = 3;
562 Index this_ss = 0, this_se = 0;
564 for (Index i_ss = 0; i_ss < scat_data.
nelem(); i_ss++)
565 for (Index i_se = scat_data[i_ss].nelem() - 1; i_se >= 0; i_se--)
570 if (nang < scat_data[i_ss][i_se].za_grid.
nelem()) {
571 nang = scat_data[i_ss][i_se].za_grid.
nelem();
575 pfct_angs = scat_data[this_ss][this_se].za_grid;
576 }
else if (Npfct < min_nang) {
578 os <<
"Number of requested angular grid points (Npfct=" << Npfct
579 <<
") is insufficient.\n"
580 <<
"At least " << min_nang <<
" points required.\n";
581 throw runtime_error(os.str());
589 ConstMatrixView pnd_profiles,
590 ConstVectorView t_profile,
591 ConstVectorView pfct_angs,
593 const Index Np_cloud = pnd_profiles.ncols();
594 const Index nang = pfct_angs.nelem();
600 Vector T_array{t_profile[Range(cloudbox_limits[0], Np_cloud)]};
601 Matrix idir_array(1, 2, 0.);
603 Matrix pdir_array(nang, 2, 0.);
604 pdir_array(joker, 0) = pfct_angs;
607 ArrayOfArrayOfTensor6 pha_mat_Nse;
610 ArrayOfTensor6 pha_mat_ssbulk;
612 Tensor6 pha_mat_bulk;
633 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
635 pha_bulk_par(joker, Range(cloudbox_limits[0], Np_cloud), joker) =
636 pha_mat_bulk(joker, joker, joker, 0, 0, 0);
640 ConstTensor3View pha_bulk_par,
641 ConstMatrixView ext_bulk_par,
642 ConstMatrixView abs_bulk_par,
644 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
645 const Index Np = pha_bulk_par.nrows();
646 const Index nf = pha_bulk_par.npages();
647 Index nang = pha_bulk_par.ncols();
656 for (Index ip = cloudbox_limits[0]; ip < Np_cloud - 1; ip++)
657 for (Index f_index = 0; f_index < nf; f_index++) {
661 (ext_bulk_par(f_index, ip) + ext_bulk_par(f_index, ip + 1)) -
662 (abs_bulk_par(f_index, ip) + abs_bulk_par(f_index, ip + 1));
666 for (Index ia = 0; ia < nang; ia++)
667 pfct_bulk_par(f_index, Np - 2 - ip, ia) +=
668 pha_bulk_par(f_index, ip, ia) + pha_bulk_par(f_index, ip + 1, ia);
669 pfct_bulk_par(f_index, Np - 2 - ip, joker) *= 4 *
PI / sca;
675 ConstTensor3View pfct_bulk_par,
676 ConstVectorView pfct_angs,
677 const Index& Nlegendre) {
678 const Index nf = pfct_bulk_par.npages();
679 const Index nlyr = pfct_bulk_par.nrows();
680 const Index nang = pfct_bulk_par.ncols();
688 Numeric pfct_threshold = 0.1;
694 Vector
u(nang), adu(nang - 1), ad_angs(nang - 1);
695 Tensor3 px(nang - 1, Nlegendre, 2, 0.);
696 u[0] = cos(pfct_angs[0] *
PI / 180.);
697 px(joker, 0, joker) = 1.;
698 for (Index ia = 1; ia < nang; ia++) {
699 u[ia] = cos(pfct_angs[ia] *
PI / 180.);
700 adu[ia - 1] = abs(
u[ia] -
u[ia - 1]);
702 abs(pfct_angs[ia] *
PI / 180. - pfct_angs[ia - 1] *
PI / 180.);
703 px(ia - 1, 1, 0) =
u[ia - 1];
704 px(ia - 1, 1, 1) =
u[ia];
705 for (Index l = 2; l < Nlegendre; l++) {
706 Numeric dl = (double)l;
707 px(ia - 1, l, 0) = (2 * dl - 1) / dl *
u[ia - 1] * px(ia - 1, l - 1, 0) -
708 (dl - 1) / dl * px(ia - 1, l - 2, 0);
709 px(ia - 1, l, 1) = (2 * dl - 1) / dl *
u[ia] * px(ia - 1, l - 1, 1) -
710 (dl - 1) / dl * px(ia - 1, l - 2, 1);
714 for (Index il = 0; il < nlyr; il++)
715 if (sum(pfct_bulk_par(joker, il, 0)) != 0.)
716 for (Index f_index = 0; f_index < nf; f_index++) {
717 if (pfct_bulk_par(f_index, il, 0) != 0) {
718 Vector pfct{pfct_bulk_par(f_index, il, joker)};
726 for (Index ia = 0; ia < nang - 1; ia++) {
727 pint += 0.5 * ad_angs[ia] *
728 (pfct[ia] * sin(pfct_angs[ia] *
PI / 180.) +
729 pfct[ia + 1] * sin(pfct_angs[ia + 1] *
PI / 180.));
732 if (abs(pint / 2. - 1.) > pfct_threshold) {
734 os <<
"Phase function normalization deviates from expected value by\n"
735 << 1e2 * pint / 2. - 1e2 <<
"(allowed: " << pfct_threshold * 1e2
737 <<
"Occurs at layer #" << il <<
" and frequency #" << f_index
739 <<
"Something is wrong with your scattering data. Check!\n";
740 throw runtime_error(os.str());
746 pmom(f_index, il, 0) = 1.;
747 for (Index ia = 0; ia < nang - 1; ia++) {
748 for (Index l = 1; l < Nlegendre; l++) {
749 pmom(f_index, il, l) +=
751 (px(ia, l, 0) * pfct[ia] * sin(pfct_angs[ia] *
PI / 180.) +
752 px(ia, l, 1) * pfct[ia + 1] *
753 sin(pfct_angs[ia + 1] *
PI / 180.));
761 const MatrixView& ext_bulk,
762 const MatrixView& abs_bulk) {
763 const Index nf = ext_bulk.nrows();
764 const Index Np = ext_bulk.ncols();
774 for (Index ip = 0; ip < Np - 1; ip++)
776 for (Index f_index = 0; f_index < nf; f_index++) {
778 0.5 * (ext_bulk(f_index, ip) - abs_bulk(f_index, ip) +
779 ext_bulk(f_index, ip + 1) - abs_bulk(f_index, ip + 1));
781 sca_bulk_layer(f_index, Np - 2 - ip) = sca;
792#define MAX_WARNINGS 100
798 static int warning_limit = FALSE, num_warnings = 0;
802 if (warning_limit)
return;
806 out1 <<
"DISORT WARNING >>> " << messag <<
"\n";
809 out1 <<
"DISORT TOO MANY WARNING MESSAGES -- They will no longer be "
810 "printed <<<<<<<\n\n";
811 warning_limit = TRUE;
826 const int maxmsg = 50;
827 static int nummsg = 0;
830 if (quiet != QUIET) {
833 out1 <<
" **** Input variable " << varnam <<
" in error ****\n";
834 if (nummsg == maxmsg) {
835 c_errmsg(
"Too many input errors. Aborting...", DS_ERROR);
844 if (quiet != QUIET) {
847 out1 <<
" **** Symbolic dimension " << dimnam
848 <<
" should be increased to at least " << minval <<
" ****\n";
865 ConstVectorView p_grid,
866 ConstVectorView z_profile,
867 const Numeric& z_surface,
868 ConstVectorView t_profile,
869 ConstMatrixView vmr_profiles,
870 ConstMatrixView pnd_profiles,
873 if (abs(z_surface - z_profile[0]) < 1e-3) {
879 cboxlims = cloudbox_limits;
885 Index np = p_grid.
nelem(), ifirst = 0;
887 for (; z_surface >= z_profile[ifirst + 1]; ifirst++) {
891 Range ind(ifirst, np);
895 vmr = vmr_profiles(joker, ind);
900 gridpos(gp[0], z_profile, z_surface);
904 t[0] =
interp(itw, t, gp[0]);
905 for (
int i = 0; i < vmr.nrows(); i++) {
906 vmr(i, 0) =
interp(itw, vmr(i, joker), gp[0]);
912 itw2p(ExhaustiveVectorView{p[0]}, p, gp, itw2);
914 cboxlims = cloudbox_limits;
915 if (ifirst < cloudbox_limits[0]) {
916 cboxlims[0] -= ifirst;
917 cboxlims[1] -= ifirst;
921 ncboxremoved = ifirst - cboxlims[0];
923 cboxlims[1] = cloudbox_limits[1] - cloudbox_limits[0] - ncboxremoved;
924 ind = Range(ncboxremoved, cboxlims[1] + 1);
925 pnd = pnd_profiles(joker, ind);
926 gp[0].idx -= cloudbox_limits[0] + ncboxremoved;
927 for (
int i = 0; i < pnd.nrows(); i++) {
928 pnd(i, 0) =
interp(itw, pnd(i, joker), gp[0]);
935 Tensor7& cloudbox_field,
936 ArrayOfMatrix& disort_aux,
937 ConstVectorView f_grid,
938 ConstVectorView p_grid,
939 ConstVectorView z_profile,
940 const Numeric& z_surface,
941 ConstVectorView t_profile,
942 ConstMatrixView vmr_profiles,
943 ConstMatrixView pnd_profiles,
946 const Agenda& propmat_clearsky_agenda,
947 const Agenda& gas_scattering_agenda,
949 const Numeric& surface_skin_t,
950 const Vector& surface_scalar_reflectivity,
951 ConstVectorView za_grid,
952 ConstVectorView aa_grid,
953 ConstVectorView sun_rte_los,
954 const Index& gas_scattering_do,
955 const Index& suns_do,
957 const Numeric& scale_factor,
958 const Index& nstreams,
960 const Index& only_tro,
962 const Index& emission,
963 const Index& intensity_correction,
995 const Index nf = f_grid.nelem();
1008 nphi = aa_grid.nelem();
1010 phi0 = sun_rte_los[1];
1017 ds.flag.prnt[0] = FALSE;
1018 ds.flag.prnt[1] = FALSE;
1019 ds.flag.prnt[2] = FALSE;
1020 ds.flag.prnt[3] = FALSE;
1021 ds.flag.prnt[4] = TRUE;
1023 ds.flag.usrtau = FALSE;
1024 ds.flag.usrang = TRUE;
1025 ds.flag.spher = FALSE;
1026 ds.flag.general_source = FALSE;
1027 ds.flag.output_uum = FALSE;
1029 ds.nlyr =
static_cast<int>(p.nelem() - 1);
1031 ds.flag.brdf_type = BRDF_NONE;
1033 ds.flag.ibcnd = GENERAL_BC;
1034 ds.flag.usrang = TRUE;
1037 ds.flag.planck = TRUE;
1039 ds.flag.planck = FALSE;
1041 ds.flag.onlyfl = FALSE;
1042 ds.flag.lamber = TRUE;
1043 ds.flag.quiet = FALSE;
1044 if (intensity_correction) {
1045 ds.flag.intensity_correction = TRUE;
1046 ds.flag.old_intensity_correction = FALSE;
1048 ds.flag.intensity_correction = FALSE;
1049 ds.flag.old_intensity_correction = FALSE;
1052 ds.nstr =
static_cast<int>(nstreams);
1053 ds.nphase = ds.nstr;
1056 ds.numu =
static_cast<int>(za_grid.nelem());
1057 ds.nphi =
static_cast<int>(nphi);
1058 Index Nlegendre = nstreams + 1;
1061 c_disort_state_alloc(&ds);
1062 c_disort_out_alloc(&ds, &out);
1072 for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i];
1074 if (ds.flag.planck==TRUE){
1075 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1079 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] *
PI / 180);
1082 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1083 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1088 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1089 Index nf_ssd = scat_data[0][0].f_grid.
nelem();
1090 Tensor3 pha_bulk_par;
1092 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1096 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1104 for (Index iss = 0; iss < scat_data.
nelem(); iss++) {
1105 const Index nse = scat_data[iss].
nelem();
1111 pnd(Range(iflat, nse), joker),
1118 get_angs(pfct_angs, scat_data, Npfct);
1119 nang = pfct_angs.nelem();
1121 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1124 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1125 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1128 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1129 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1132 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1133 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1135 if (gas_scattering_do) {
1139 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1143 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1144 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1145 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1148 sca_coeff_gas_layer,
1149 sca_coeff_gas_level,
1155 gas_scattering_agenda);
1159 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1162 ext_bulk_par += sca_coeff_gas_level;
1166 Matrix dtauc(nf, ds.nlyr);
1168 Matrix ssalb(nf, ds.nlyr);
1169 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1188 ds.bc.btemp = surface_skin_t;
1191 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1192 snprintf(ds.header, 128,
"ARTS Calc f_index = %" PRId64, f_index);
1194 std::memcpy(ds.dtauc,
1195 dtauc(f_index, joker).unsafe_data_handle(),
1196 sizeof(Numeric) * ds.nlyr);
1197 std::memcpy(ds.ssalb,
1198 ssalb(f_index, joker).unsafe_data_handle(),
1199 sizeof(Numeric) * ds.nlyr);
1202 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. *
SPEED_OF_LIGHT);
1203 ds.wvnmhi += ds.wvnmhi * 1e-7;
1204 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1207 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1211 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1214 ds.bc.fbeam = fbeam;
1216 std::memcpy(ds.pmom,
1217 pmom(f_index, joker, joker).unsafe_data_handle(),
1218 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1220 enum class Status { FIRST_TRY, RETRY, SUCCESS };
1221 Status tries = Status::FIRST_TRY;
1222 const Numeric eps = 2e-4;
1225 c_disort(&ds, &out);
1226 tries = Status::SUCCESS;
1227 }
catch (
const std::runtime_error& e) {
1229 if (tries == Status::FIRST_TRY) {
1231 if (umu0 < 1 - eps) {
1233 }
else if (umu0 > 1 - eps) {
1237 const Numeric shift =
1241 <<
"Solar zenith angle coincided with one of the quadrature angles\n"
1242 <<
"We needed to shift the solar sun angle by " << shift
1246 tries = Status::RETRY;
1250 }
while (tries != Status::SUCCESS);
1252 for (Index i = 0; i < ds.nphi; i++) {
1253 for (Index j = 0; j < ds.numu; j++) {
1254 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1255 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, i, 0) =
1256 out.uu[j + ((ds.nlyr - k - cboxlims[0]) + i * (ds.nlyr + 1)) *
1262 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1263 cloudbox_field(f_index, k, 0, 0, j, i, 0) =
1264 cloudbox_field(f_index, k + 1, 0, 0, j, i, 0);
1271 disort_aux.resize(disort_aux_vars.
nelem());
1274 for (Index i = 0; i < disort_aux_vars.
nelem(); i++) {
1277 if (disort_aux_vars[i] ==
"Layer optical thickness"){
1279 Matrix deltatau(nf, ds.nlyr, 0);
1280 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1281 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1282 deltatau(f_index, k - 1 + ncboxremoved) =
1283 dtauc(f_index, ds.nlyr - k + ncboxremoved);
1286 disort_aux[cnt] = deltatau;
1288 else if (disort_aux_vars[i] ==
"Single scattering albedo"){
1290 Matrix snglsctalbedo(nf, ds.nlyr,0);
1291 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1292 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1293 snglsctalbedo(f_index, k - 1 + ncboxremoved) =
1294 ssalb(f_index, ds.nlyr - k + ncboxremoved);
1297 disort_aux[cnt]=snglsctalbedo;
1299 else if (disort_aux_vars[i] ==
"Direct beam") {
1301 Matrix directbeam(nf, ds.nlyr + 1, 0);
1303 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1304 directbeam(f_index, cboxlims[1] - cboxlims[0] + ncboxremoved) =
1305 suns[0].spectrum(f_index, 0)/
PI;
1307 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1308 directbeam(f_index, k - 1 + ncboxremoved) =
1309 directbeam(f_index, k + ncboxremoved) *
1310 exp(-dtauc(f_index, ds.nlyr - k + ncboxremoved)/umu0);
1314 disort_aux[cnt]=directbeam;
1317 "The only allowed strings in *disort_aux_vars* are:\n"
1318 " \"Layer optical thickness\"\n"
1319 " \"Single scattering albedo\"\n"
1320 " \"Direct beam\"\n"
1321 "but you have selected: \"", disort_aux_vars[i],
"\"\n");
1326 c_disort_out_free(&ds, &out);
1327 c_disort_state_free(&ds);
1335 Tensor5& spectral_irradiance_field,
1336 ArrayOfMatrix& disort_aux,
1337 ConstVectorView f_grid,
1338 ConstVectorView p_grid,
1339 ConstVectorView z_profile,
1340 const Numeric& z_surface,
1341 ConstVectorView t_profile,
1342 ConstMatrixView vmr_profiles,
1343 ConstMatrixView pnd_profiles,
1346 const Agenda& propmat_clearsky_agenda,
1347 const Agenda& gas_scattering_agenda,
1349 const Numeric& surface_skin_t,
1350 const Vector& surface_scalar_reflectivity,
1351 ConstVectorView sun_rte_los,
1352 const Index& gas_scattering_do,
1353 const Index& suns_do,
1355 const Numeric& scale_factor,
1356 const Index& nstreams,
1358 const Index& only_tro,
1360 const Index& emission,
1361 const Index& intensity_correction,
1394 const Index nf = f_grid.nelem();
1406 phi0 = sun_rte_los[1];
1413 ds.flag.prnt[0] = FALSE;
1414 ds.flag.prnt[1] = FALSE;
1415 ds.flag.prnt[2] = FALSE;
1416 ds.flag.prnt[3] = FALSE;
1417 ds.flag.prnt[4] = TRUE;
1419 ds.flag.usrtau = FALSE;
1420 ds.flag.usrang = TRUE;
1421 ds.flag.spher = FALSE;
1422 ds.flag.general_source = FALSE;
1423 ds.flag.output_uum = FALSE;
1425 ds.nlyr =
static_cast<int>(p.nelem() - 1);
1427 ds.flag.brdf_type = BRDF_NONE;
1429 ds.flag.ibcnd = GENERAL_BC;
1430 ds.flag.usrang = FALSE;
1433 ds.flag.planck = TRUE;
1435 ds.flag.planck = FALSE;
1437 ds.flag.onlyfl = TRUE;
1438 ds.flag.lamber = TRUE;
1439 ds.flag.quiet = FALSE;
1440 if (intensity_correction) {
1441 ds.flag.intensity_correction = TRUE;
1442 ds.flag.old_intensity_correction = FALSE;
1444 ds.flag.intensity_correction = FALSE;
1445 ds.flag.old_intensity_correction = FALSE;
1448 ds.nstr =
static_cast<int>(nstreams);
1449 ds.nphase = ds.nstr;
1457 Index Nlegendre = nstreams + 1;
1460 c_disort_state_alloc(&ds);
1461 c_disort_out_alloc(&ds, &out);
1471 if (ds.flag.planck==TRUE){
1472 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1476 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1477 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1482 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1483 Index nf_ssd = scat_data[0][0].f_grid.
nelem();
1484 Tensor3 pha_bulk_par;
1486 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1490 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1498 for (Index iss = 0; iss < scat_data.
nelem(); iss++) {
1499 const Index nse = scat_data[iss].
nelem();
1505 pnd(Range(iflat, nse), joker),
1512 get_angs(pfct_angs, scat_data, Npfct);
1513 nang = pfct_angs.nelem();
1515 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1518 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1519 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1522 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1523 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1526 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1527 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1529 if (gas_scattering_do) {
1533 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1537 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1538 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1539 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1542 sca_coeff_gas_layer,
1543 sca_coeff_gas_level,
1549 gas_scattering_agenda);
1553 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1556 ext_bulk_par += sca_coeff_gas_level;
1560 Matrix dtauc(nf, ds.nlyr);
1562 Matrix ssalb(nf, ds.nlyr);
1563 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1582 ds.bc.btemp = surface_skin_t;
1585 Matrix spectral_direct_irradiance_field;
1586 Matrix dFdtau(nf, ds.nlyr+1,0);
1587 Matrix deltatau(nf, ds.nlyr,0);
1588 Matrix snglsctalbedo(nf, ds.nlyr,0);
1592 spectral_direct_irradiance_field.resize(nf, ds.nlyr+1);
1593 spectral_direct_irradiance_field = 0;
1596 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1597 snprintf(ds.header, 128,
"ARTS Calc f_index = %" PRId64, f_index);
1599 std::memcpy(ds.dtauc,
1600 dtauc(f_index, joker).unsafe_data_handle(),
1601 sizeof(Numeric) * ds.nlyr);
1602 std::memcpy(ds.ssalb,
1603 ssalb(f_index, joker).unsafe_data_handle(),
1604 sizeof(Numeric) * ds.nlyr);
1607 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. *
SPEED_OF_LIGHT);
1608 ds.wvnmhi += ds.wvnmhi * 1e-7;
1609 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1612 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1616 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1619 ds.bc.fbeam = fbeam;
1621 std::memcpy(ds.pmom,
1622 pmom(f_index, joker, joker).unsafe_data_handle(),
1623 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1625 c_disort(&ds, &out);
1628 const Numeric conv_fac=(ds.wvnmhi - ds.wvnmlo) * (100 *
SPEED_OF_LIGHT);
1630 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1633 spectral_direct_irradiance_field(f_index, k + ncboxremoved) =
1634 -out.rad[ds.nlyr - k - cboxlims[0]].rfldir/conv_fac;
1637 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1638 -(out.rad[ds.nlyr - k - cboxlims[0]].rfldir +
1639 out.rad[ds.nlyr - k - cboxlims[0]].rfldn)/conv_fac;
1643 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1644 -out.rad[ds.nlyr - k - cboxlims[0]].rfldn/conv_fac;
1648 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 1) =
1649 out.rad[ds.nlyr - k - cboxlims[0]].flup/conv_fac;
1652 dFdtau(f_index, k + ncboxremoved) =
1653 -out.rad[ds.nlyr - k - cboxlims[0]].dfdt;
1658 deltatau(f_index, k - 1 + ncboxremoved) = ds.dtauc[ds.nlyr - k - 1 - cboxlims[0]];
1659 snglsctalbedo(f_index, k - 1 + ncboxremoved) = ds.ssalb[ds.nlyr - k - 1 - cboxlims[0]];
1665 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1666 spectral_irradiance_field(f_index, k, 0, 0, joker) =
1667 spectral_irradiance_field(f_index, k + 1, 0, 0, joker);
1670 spectral_direct_irradiance_field(f_index, k) =
1671 spectral_direct_irradiance_field(f_index, k + 1);
1673 dFdtau(f_index, k) = dFdtau(f_index, k + 1);
1678 disort_aux.resize(disort_aux_vars.
nelem());
1681 for (Index i = 0; i < disort_aux_vars.
nelem(); i++) {
1684 if (disort_aux_vars[i] ==
"Layer optical thickness"){
1686 disort_aux[cnt]=deltatau;
1688 else if (disort_aux_vars[i] ==
"Single scattering albedo"){
1690 disort_aux[cnt]=snglsctalbedo;
1692 else if (disort_aux_vars[i] ==
"Direct downward spectral irradiance") {
1694 disort_aux[cnt] = spectral_direct_irradiance_field;
1696 else if (disort_aux_vars[i] ==
"dFdtau") {
1698 disort_aux[cnt] = dFdtau;
1702 "The only allowed strings in *disort_aux_vars* are:\n"
1703 " \"Layer optical thickness\"\n"
1704 " \"Single scattering albedo\"\n"
1705 " \"Direct downward spectral irradiance\"\n"
1707 "but you have selected: \"", disort_aux_vars[i],
"\"\n");
1714 c_disort_out_free(&ds, &out);
1715 c_disort_state_free(&ds);
1727 const Agenda& surface_rtprop_agenda,
1728 ConstVectorView f_grid,
1729 ConstVectorView scat_za_grid,
1730 const Numeric& surf_alt,
1758 chk_not_empty(
"surface_rtprop_agenda", surface_rtprop_agenda);
1760 const Index nf = f_grid.nelem();
1762 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
1763 if (frza == scat_za_grid.nelem()) {
1765 os <<
"No upwelling direction found in scat_za_grid.\n";
1766 throw runtime_error(os.str());
1768 const Index nrza = scat_za_grid.nelem() - frza;
1771 Matrix dir_refl_coeff(nrza, nf, 0.);
1774 Vector rtp_pos(1, surf_alt);
1778 for (Index rza = 0; rza < nrza; rza++) {
1780 Numeric surface_skin_t;
1782 Tensor4 surface_rmatrix;
1783 Matrix surface_emission;
1785 Vector rtp_los(1, scat_za_grid[rza + frza]);
1786 out2 <<
"Doing reflected dir #" << rza <<
" at " << rtp_los[0] <<
" degs\n";
1796 surface_rtprop_agenda);
1797 if (surface_los.nrows() > 1) {
1799 os <<
"For this method, *surface_rtprop_agenda* must be set up to\n"
1800 <<
"return zero or one direction in *surface_los*.";
1801 throw runtime_error(os.str());
1805 btemp = surface_skin_t;
1806 else if (surface_skin_t != btemp) {
1808 os <<
"Something went wrong.\n"
1809 <<
" *surface_rtprop_agenda* returned different surface_skin_t\n"
1810 <<
" for different LOS.\n";
1811 throw runtime_error(os.str());
1815 if (!surface_los.empty()) {
1816 for (Index f_index = 0; f_index < nf; f_index++)
1817 dir_refl_coeff(rza, f_index) =
1818 sum(surface_rmatrix(joker, f_index, 0, 0));
1820 out2 <<
" directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1824 if (btemp < 0. || btemp > 1000.) {
1826 os <<
"Surface temperature has been derived as " << btemp <<
" K,\n"
1827 <<
"which is not considered a meaningful value.\n";
1828 throw runtime_error(os.str());
1835 Vector surf_int_grid(nrza + 1);
1840 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1842 os <<
"Looks like scat_za_grid contains the 90deg direction,\n"
1843 <<
"which it shouldn't for running Disort.\n";
1844 throw runtime_error(os.str());
1846 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1847 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1848 const Numeric ok_extrapol = 0.5;
1849 if ((za_extrapol - ok_extrapol) > 1e-6) {
1851 os <<
"Extrapolation range from shallowest scat_za_grid point\n"
1852 <<
"to horizon is too big.\n"
1853 <<
" Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1854 <<
", which is " << za_extrapol - 0.5 <<
" too big.\n";
1855 throw runtime_error(os.str());
1857 if (!is_same_within_epsilon(
1858 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1860 os <<
"Looks like last point in scat_za_grid is not 180deg.\n";
1861 throw runtime_error(os.str());
1864 surf_int_grid[0] = 90.;
1865 surf_int_grid[nrza] = 180.;
1866 for (Index rza = 1; rza < nrza; rza++)
1867 surf_int_grid[rza] =
1868 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1872 for (Index rza = 0; rza < nrza; rza++) {
1877 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1883 dir_refl_coeff(rza, joker) *=
w;
1888 for (Index f_index = 0; f_index < nf; f_index++) {
1889 albedo[f_index] = sum(dir_refl_coeff(joker, f_index));
1890 out2 <<
"at f=" << f_grid[f_index] * 1e-9
1891 <<
" GHz, ending up with albedo=" << albedo[f_index] <<
"\n";
1892 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1894 os <<
"Something went wrong: Albedo must be inside [0,1],\n"
1895 <<
" but is not at freq #" << f_index <<
" , where it is "
1896 << albedo[f_index] <<
".\n";
1897 throw runtime_error(os.str());
1907 const Agenda& surface_rtprop_agenda,
1908 ConstVectorView f_grid,
1909 const Numeric& surf_alt,
1910 const Numeric& inc_angle) {
1912 Vector rtp_pos(1, surf_alt);
1913 Vector rtp_los(1, 180-inc_angle);
1915 Numeric surface_skin_t;
1917 Tensor4 surface_rmatrix;
1918 Matrix surface_emission;
1928 surface_rtprop_agenda);
1930 if (surface_los.nrows() > 1) {
1932 os <<
"For this method, *surface_rtprop_agenda* must be set up to\n"
1933 <<
"return zero or one direction in *surface_los*.";
1934 throw runtime_error(os.str());
1936 if (surface_los.empty()) {
1939 albedo[joker] = surface_rmatrix(0,joker,0,0);
1942 btemp = surface_skin_t;
1943 if (btemp < 0. || btemp > 1000.) {
1945 os <<
"Surface temperature has been derived as " << btemp <<
" K,\n"
1946 <<
"which is not considered a meaningful value.\n";
1947 throw runtime_error(os.str());
1953#pragma GCC diagnostic pop
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
Constants of physical expressions as constexpr.
void gas_scattering_agendaExecute(Workspace &ws, PropagationMatrix &gas_scattering_coef, TransmissionMatrix &gas_scattering_mat, Vector &gas_scattering_fct_legendre, const Vector &f_grid, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Vector &gas_scattering_los_in, const Vector &gas_scattering_los_out, const Index gas_scattering_output_type, const Agenda &input_agenda)
void propmat_clearsky_agendaExecute(Workspace &ws, PropagationMatrix &propmat_clearsky, StokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_source_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfSpeciesTag &select_abs_species, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &input_agenda)
void surface_rtprop_agendaExecute(Workspace &ws, Numeric &surface_skin_t, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Index nelem() const ARTS_NOEXCEPT
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
void surf_albedoCalc(Workspace &ws, VectorView albedo, Numeric &btemp, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, ConstVectorView scat_za_grid, const Numeric &surf_alt, const Verbosity &verbosity)
surf_albedoCalc
constexpr Numeric PLANCK_CONST
void c_errmsg(const char *messag, int type)
Verbosity enabled replacement for the original cdisort function.
void get_pfct(Tensor3 &pfct_bulk_par, ConstTensor3View pha_bulk_par, ConstMatrixView ext_bulk_par, ConstMatrixView abs_bulk_par, const ArrayOfIndex &cloudbox_limits)
get_pfct.
void reduced_1datm(Vector &p, Vector &z, Vector &t, Matrix &vmr, Matrix &pnd, ArrayOfIndex &cboxlims, Index &ncboxremoved, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfIndex &cloudbox_limits)
reduced_1datm
constexpr Numeric SPEED_OF_LIGHT
void run_cdisort(Workspace &ws, Tensor7 &cloudbox_field, ArrayOfMatrix &disort_aux, ConstVectorView f_grid, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &gas_scattering_agenda, const ArrayOfIndex &cloudbox_limits, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView sun_rte_los, const Index &gas_scattering_do, const Index &suns_do, const ArrayOfString &disort_aux_vars, const Numeric &scale_factor, const Index &nstreams, const Index &Npfct, const Index &only_tro, const Index &quiet, const Index &emission, const Index &intensity_correction, const Verbosity &verbosity)
Calculate doit_i_field with Disort including a sun source.
constexpr Numeric DEG2RAD
void get_gas_scattering_properties(Workspace &ws, MatrixView sca_coeff_gas, MatrixView sca_coeff_gas_level, MatrixView pfct_gas, const ConstVectorView &f_grid, const VectorView &p, const VectorView &t, const MatrixView &vmr, const Agenda &gas_scattering_agenda)
get_gas_scattering_properties
void get_parZ(Tensor3 &pha_bulk_par, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstMatrixView pnd_profiles, ConstVectorView t_profile, ConstVectorView pfct_angs, const ArrayOfIndex &cloudbox_limits)
get_parZ.
void run_cdisort_flux(Workspace &ws, Tensor5 &spectral_irradiance_field, ArrayOfMatrix &disort_aux, ConstVectorView f_grid, ConstVectorView p_grid, ConstVectorView z_profile, const Numeric &z_surface, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstMatrixView pnd_profiles, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfSun &suns, const Agenda &propmat_clearsky_agenda, const Agenda &gas_scattering_agenda, const ArrayOfIndex &cloudbox_limits, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, ConstVectorView sun_rte_los, const Index &gas_scattering_do, const Index &suns_do, const ArrayOfString &disort_aux_vars, const Numeric &scale_factor, const Index &nstreams, const Index &Npfct, const Index &only_tro, const Index &quiet, const Index &emission, const Index &intensity_correction, const Verbosity &verbosity)
Calculate spectral_irradiance_field with Disort including a sun source.
void get_disortsurf_props(Vector &albedo, Numeric &btemp, ConstVectorView f_grid, const Numeric &surface_skin_t, ConstVectorView surface_scalar_reflectivity)
get_disortsurf_props.
void get_scat_bulk_layer(MatrixView sca_bulk_layer, const MatrixView &ext_bulk, const MatrixView &abs_bulk)
get_scat_bulk_layer
void surf_albedoCalcSingleAngle(Workspace &ws, VectorView albedo, Numeric &btemp, const Agenda &surface_rtprop_agenda, ConstVectorView f_grid, const Numeric &surf_alt, const Numeric &inc_angle)
surf_albedoCalcSingleAngle
void check_disort_irradiance_input(const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &nstreams)
check_disort_input.
int c_write_bad_var(int quiet, const char *varnam)
Verbosity enabled replacement for the original cdisort function.
void get_dtauc_ssalb(MatrixView dtauc, MatrixView ssalb, ConstMatrixView ext_bulk_gas, ConstMatrixView ext_bulk_par, ConstMatrixView abs_bulk_par, ConstVectorView z_profile)
get_dtauc_ssalb
void check_disort_input(const Index &cloudbox_on, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &atmosphere_dim, const Index &stokes_dim, const ArrayOfIndex &cloudbox_limits, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstVectorView za_grid, const Index &nstreams)
check_disort_input.
void init_ifield(Tensor7 &cloudbox_field, const Vector &f_grid, const ArrayOfIndex &cloudbox_limits, const Index &n_za, const Index &n_aa, const Index &stokes_dim)
init_ifield.
thread_local Verbosity disort_verbosity
void get_angs(Vector &pfct_angs, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &Npfct)
get_angs.
void get_pmom(Tensor3View pmom, ConstTensor3View pfct_bulk_par, ConstVectorView pfct_angs, const Index &Nlegendre)
get_pmom
void get_gasoptprop(Workspace &ws, MatrixView ext_bulk_gas, const Agenda &propmat_clearsky_agenda, ConstVectorView t_profile, ConstMatrixView vmr_profiles, ConstVectorView p_grid, ConstVectorView f_grid)
get_gasoptprop.
constexpr Numeric COSMIC_BG_TEMP
void get_paroptprop(MatrixView ext_bulk_par, MatrixView abs_bulk_par, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstMatrixView pnd_profiles, ConstVectorView t_profile, ConstVectorView DEBUG_ONLY(p_grid), const ArrayOfIndex &cloudbox_limits, ConstVectorView f_grid)
int c_write_too_small_dim(int quiet, const char *dimnam, int minval)
Verbosity enabled replacement for the original cdisort function.
void add_normed_phase_functions(Tensor3View pfct1, const MatrixView &sca1, const MatrixView &pfct2, const MatrixView &sca2)
add_normed_phase_functions
Functions for disort interface.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Declarations having to do with the four output streams.
constexpr Numeric cosmic_microwave_background_temperature
Global constant, Planck temperature for cosmic background radiation [K].
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric speed_of_light
Speed of light [m/s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 201...
constexpr Numeric planck_constant
Planck constant [J s] From: https://en.wikipedia.org/wiki/2019_redefinition_of_SI_base_units Date: 20...
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto acosd(auto x) noexcept
Returns rad2deg of the arc-cosine of the input.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void pha_mat_ScatSpecBulk(ArrayOfTensor6 &pha_mat, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor6 &pha_mat_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk phase matrices.
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
void pha_mat_NScatElems(ArrayOfArrayOfTensor6 &pha_mat, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_index, const Index &t_interp_order)
Phase matrices from all scattering elements.
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
void ext_abs_pfun_from_tro(MatrixView ext_data, MatrixView abs_data, Tensor3View pfun_data, const ArrayOfSingleScatteringData &scat_data, const Index &iss, ConstMatrixView pnd_data, ArrayOfIndex &cloudbox_limits, ConstVectorView T_grid, ConstVectorView sa_grid)
Extinction, absorption and phase function for one scattering species, 1D and TRO.
Declaration of functions in rte.cc.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
This file contains basic functions to handle XML data files.