ARTS 2.5.11 (git: 6827797f)
disort.cc
Go to the documentation of this file.
1
11#include "disort.h"
12
13#include <cinttypes>
14#include <cmath>
15#include <cstdlib>
16#include <iostream>
17#include <stdexcept>
18
19#include "agenda_class.h"
20#include "array.h"
21#include "arts_constants.h"
22#include "auto_md.h"
23#include "check_input.h"
24#include "arts_conversions.h"
25#include "debug.h"
26
27#if not ARTS_LGPL
28extern "C" {
29#include "cdisort.h"
30}
31#endif
32
33#include "disort.h"
34#include "interpolation.h"
35#include "logic.h"
36#include "math_funcs.h"
37#include "messages.h"
38#include "rte.h"
39#include "xml_io.h"
40
41#if ARTS_LGPL
42#pragma GCC diagnostic push
43#pragma GCC diagnostic ignored "-Wunused-parameter"
44#endif
45
46inline constexpr Numeric PI=Constant::pi;
47inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
48inline constexpr Numeric PLANCK_CONST=Constant::planck_constant;
51
52void add_normed_phase_functions(Tensor3View pfct1,
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();
59
60
61 ARTS_ASSERT(pfct2.nrows() == nr1);
62
63 ARTS_ASSERT(pfct2.ncols() == nc1);
64
65 for (Index i = 0; i < np1; i++) { // frequncy loop
66 for (Index j = 0; j < nr1 ; j++) { // layer loop
67 for (Index k = 0; k < nc1; k++) // polynomial loop
68
69 pfct1(i, j, k) =
70 (sca1(i, j) * pfct1(i, j, k) + sca2(i, j) * pfct2(j, k)) /
71 (sca1(i, j) + sca2(i, j));
72 }
73 }
74}
75
76void check_disort_input( // Input
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,
84 const ArrayOfIndex& cloudbox_limits,
86 ConstVectorView za_grid,
87 const Index& nstreams) {
88 if (!cloudbox_on) {
89 throw runtime_error(
90 "Cloudbox is off, no scattering calculations to be"
91 "performed.");
92 }
93
94 if (atmfields_checked != 1)
95 throw runtime_error(
96 "The atmospheric fields must be flagged to have "
97 "passed a consistency check (atmfields_checked=1).");
98 if (atmgeom_checked != 1)
99 throw runtime_error(
100 "The atmospheric geometry must be flagged to have "
101 "passed a consistency check (atmgeom_checked=1).");
102 if (cloudbox_checked != 1)
103 throw runtime_error(
104 "The cloudbox must be flagged to have "
105 "passed a consistency check (cloudbox_checked=1).");
106 if (scat_data_checked != 1)
107 throw runtime_error(
108 "The scat_data must be flagged to have "
109 "passed a consistency check (scat_data_checked=1).");
110
111 if (atmosphere_dim != 1)
112 throw runtime_error(
113 "For running DISORT, atmospheric dimensionality "
114 "must be 1.\n");
115
116 if (stokes_dim < 0 || stokes_dim > 1)
117 throw runtime_error(
118 "For running DISORT, the dimension of stokes vector "
119 "must be 1.\n");
120
121 if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
122 throw runtime_error(
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*");
127
128 if (cloudbox_limits[0] != 0) {
129 ostringstream os;
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());
135 }
136
137 if (scat_data.empty())
138 throw runtime_error(
139 "No single scattering data present.\n"
140 "See documentation of WSV *scat_data* for options to "
141 "make single scattering data available.\n");
142
143 // DISORT requires even number of streams:
144 // nstreams is total number of directions, up- and downwelling, and the up-
145 // and downwelling directions need to be symmetrically distributed, i.e. same
146 // number of directions in both hemispheres is required. horizontal direction
147 // (90deg) can not be covered in a plane-parallel atmosphere.
148 if (nstreams / 2 * 2 != nstreams) {
149 ostringstream os;
150 os << "DISORT requires an even number of streams, but yours is " << nstreams
151 << ".\n";
152 throw runtime_error(os.str());
153 }
154
155 // Zenith angle grid.
156 Index nza = za_grid.nelem();
157
158 // za_grid here is only relevant to provide an i_field from which the
159 // sensor los angles can be interpolated by yCalc; it does not the determine
160 // the accuracy of the DISORT output itself at these angles. So we can only
161 // apply a very rough test here, whether the grid is appropriate. However, we
162 // set the threshold fairly high since calculation costs for a higher number
163 // of angles are negligible.
164 if (nza < 20) {
165 ostringstream os;
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());
171 }
172
173 if (za_grid[0] != 0. || za_grid[nza - 1] != 180.)
174 throw runtime_error("The range of *za_grid* must [0 180].");
175
176 if (!is_increasing(za_grid))
177 throw runtime_error("*za_grid* must be increasing.");
178
179 Index i = 1;
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");
183 i++;
184 }
185
186 // DISORT can only handle randomly oriented particles.
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;
191 if (!all_totrand) {
192 ostringstream os;
193 os << "DISORT can only handle scattering elements of type "
194 << PTYPE_TOTAL_RND << " (" << PTypeToString(PTYPE_TOTAL_RND) << "),\n"
195 << "but at least one element of other type (" << PTYPE_AZIMUTH_RND << "="
196 << PTypeToString(PTYPE_AZIMUTH_RND) << " or " << PTYPE_GENERAL << "="
197 << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
198 throw runtime_error(os.str());
199 }
200}
201
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,
208 const ArrayOfArrayOfSingleScatteringData& scat_data,
209 const Index& nstreams) {
210 if (atmfields_checked != 1)
211 throw runtime_error(
212 "The atmospheric fields must be flagged to have "
213 "passed a consistency check (atmfields_checked=1).");
214
215 if (atmgeom_checked != 1)
216 throw runtime_error(
217 "The atmospheric geometry must be flagged to have "
218 "passed a consistency check (atmgeom_checked=1).");
219
220 if (scat_data_checked != 1)
221 throw runtime_error(
222 "The scat_data must be flagged to have "
223 "passed a consistency check (scat_data_checked=1).");
224
225 if (atmosphere_dim != 1)
226 throw runtime_error(
227 "For running DISORT, atmospheric dimensionality "
228 "must be 1.\n");
229
230 if (stokes_dim < 0 || stokes_dim > 1)
231 throw runtime_error(
232 "For running DISORT, the dimension of stokes vector "
233 "must be 1.\n");
234
235 if (scat_data.empty())
236 throw runtime_error(
237 "No single scattering data present.\n"
238 "See documentation of WSV *scat_data* for options to "
239 "make single scattering data available.\n");
240
241 // DISORT requires even number of streams:
242 // nstreams is total number of directions, up- and downwelling, and the up-
243 // and downwelling directions need to be symmetrically distributed, i.e. same
244 // number of directions in both hemispheres is required. horizontal direction
245 // (90deg) can not be covered in a plane-parallel atmosphere.
246 if (nstreams / 2 * 2 != nstreams) {
247 ostringstream os;
248 os << "DISORT requires an even number of streams, but yours is " << nstreams
249 << ".\n";
250 throw runtime_error(os.str());
251 }
252
253 // DISORT can only handle randomly oriented particles.
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;
258 if (!all_totrand) {
259 ostringstream os;
260 os << "DISORT can only handle scattering elements of type "
261 << PTYPE_TOTAL_RND << " (" << PTypeToString(PTYPE_TOTAL_RND) << "),\n"
262 << "but at least one element of other type (" << PTYPE_AZIMUTH_RND << "="
263 << PTypeToString(PTYPE_AZIMUTH_RND) << " or " << PTYPE_GENERAL << "="
264 << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
265 throw runtime_error(os.str());
266 }
267}
268
269void init_ifield( // Output
270 Tensor7& cloudbox_field,
271 // Input
272 const Vector& f_grid,
273 const ArrayOfIndex& cloudbox_limits,
274 const Index& n_za,
275 const Index& n_aa,
276 const Index& stokes_dim) {
277 const Index Nf = f_grid.nelem();
278 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
279 //const Index Nza = za_grid.nelem();
280
281 // Resize and initialize radiation field in the cloudbox
282 cloudbox_field.resize(Nf, Np_cloud, 1, 1, n_za, n_aa, stokes_dim);
283 cloudbox_field = NAN;
284}
285
286void get_disortsurf_props( // Output
287 Vector& albedo,
288 Numeric& btemp,
289 // Input
290 ConstVectorView f_grid,
291 const Numeric& surface_skin_t,
292 ConstVectorView surface_scalar_reflectivity) {
293 // temperature of surface
294 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
295 ostringstream os;
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());
301 }
302 btemp = surface_skin_t;
303
304 // surface albedo
305 if (surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
306 surface_scalar_reflectivity.nelem() != 1) {
307 ostringstream os;
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());
314 }
315
316 if (min(surface_scalar_reflectivity) < 0 ||
317 max(surface_scalar_reflectivity) > 1) {
318 throw runtime_error(
319 "All values in *surface_scalar_reflectivity*"
320 " must be inside [0,1].");
321 }
322
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];
326 else
327 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
328 albedo[f_index] = surface_scalar_reflectivity[0];
329}
330
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();
339
340 ARTS_ASSERT(ext_bulk_gas.nrows() == f_grid.nelem());
341 ARTS_ASSERT(ext_bulk_gas.ncols() == Np);
342
343 // Initialization
344 ext_bulk_gas = 0.;
345
346 // making gas property output containers and input dummies
347 const EnergyLevelMap rtp_nlte_dummy;
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;
353
354 PropagationMatrix propmat_clearsky_local;
355 for (Index ip = 0; ip < Np; ip++) {
357 propmat_clearsky_local,
358 nlte_dummy,
359 partial_dummy,
360 partial_nlte_dummy,
362 {},
363 Vector{f_grid},
364 rtp_mag_dummy,
365 ppath_los_dummy,
366 p_grid[ip],
367 t_profile[ip],
368 rtp_nlte_dummy,
369 Vector{vmr_profiles(joker, ip)},
370 propmat_clearsky_agenda);
371 ext_bulk_gas(joker, ip) += propmat_clearsky_local.Kjj();
372 }
373}
374
376 MatrixView sca_coeff_gas,
377 MatrixView sca_coeff_gas_level,
378 MatrixView pfct_gas,
379 const ConstVectorView& f_grid,
380 const VectorView& p,
381 const VectorView& t,
382 const MatrixView& vmr,
383 const Agenda& gas_scattering_agenda) {
384 const Index Np = p.nelem(); // Number of pressure levels
385 const Index Nl = pfct_gas.ncols(); // Number of legendre polynomials
386 const Index Nf = f_grid.nelem(); // Number of frequencies
387
388 PropagationMatrix K_sca_gas_temp;
389 TransmissionMatrix sca_mat_dummy;
390 Vector dummy;
391 Vector sca_fct_temp;
392 Matrix pmom_gas_level( Np, Nl, 0.);
393 Index N_polys;
394
395 // calculate gas scattering properties on level
396 for (Index ip = 0; ip < Np; ip++) {
398 K_sca_gas_temp,
399 sca_mat_dummy,
400 sca_fct_temp,
401 Vector{f_grid},
402 p[ip],
403 t[ip],
404 Vector{vmr(joker, ip)},
405 dummy,
406 dummy,
407 1,
408 gas_scattering_agenda);
409
410 // gas scattering extinction
411 sca_coeff_gas_level(joker, ip) = K_sca_gas_temp.Kjj(0, 0);
412
413 // gas scattering (phase) function
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];
417 }
418 }
419
420 // layer averages
421 for (Index ip = 0; ip < Np - 1; ip++) {
422 for (Index f_index = 0; f_index < Nf; f_index++) {
423 // extinction
424 sca_coeff_gas(f_index, Np - 2 - ip) =
425 0.5 *
426 (sca_coeff_gas_level(f_index, ip) +
427 sca_coeff_gas_level(f_index, ip + 1));
428 }
429 // phase function
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));
434 }
435 }
436}
437
438void get_paroptprop(MatrixView ext_bulk_par,
439 MatrixView abs_bulk_par,
440 const ArrayOfArrayOfSingleScatteringData& scat_data,
441 ConstMatrixView pnd_profiles,
442 ConstVectorView t_profile,
443 ConstVectorView DEBUG_ONLY(p_grid),
444 const ArrayOfIndex& cloudbox_limits,
445 ConstVectorView f_grid) {
446 const Index Np_cloud = pnd_profiles.ncols();
447 DEBUG_ONLY(const Index Np = p_grid.nelem());
448 const Index nf = f_grid.nelem();
449
450 ARTS_ASSERT(ext_bulk_par.nrows() == nf);
451 ARTS_ASSERT(abs_bulk_par.nrows() == nf);
452 ARTS_ASSERT(ext_bulk_par.ncols() == Np);
453 ARTS_ASSERT(abs_bulk_par.ncols() == Np);
454
455 // Initialization
456 ext_bulk_par = 0.;
457 abs_bulk_par = 0.;
458
459 // preparing input data
460 Vector T_array{t_profile[Range(cloudbox_limits[0], Np_cloud)]};
461 Matrix dir_array(1, 2, 0.); // just a dummy. only tot_random allowed, ie.
462 // optprop are independent of direction.
463
464 // making particle property output containers
465 ArrayOfArrayOfTensor5 ext_mat_Nse;
466 ArrayOfArrayOfTensor4 abs_vec_Nse;
467 ArrayOfArrayOfIndex ptypes_Nse;
468 Matrix t_ok;
469 ArrayOfTensor5 ext_mat_ssbulk;
470 ArrayOfTensor4 abs_vec_ssbulk;
471 ArrayOfIndex ptype_ssbulk;
472 Tensor5 ext_mat_bulk; //nf,nT,ndir,nst,nst
473 Tensor4 abs_vec_bulk;
474 Index ptype_bulk;
475
476 // calculate particle optical properties
477 opt_prop_NScatElems(ext_mat_Nse,
478 abs_vec_Nse,
479 ptypes_Nse,
480 t_ok,
481 scat_data,
482 1,
483 T_array,
484 dir_array,
485 -1);
486 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
487 abs_vec_ssbulk,
488 ptype_ssbulk,
489 ext_mat_Nse,
490 abs_vec_Nse,
491 ptypes_Nse,
492 pnd_profiles,
493 t_ok);
494 opt_prop_Bulk(ext_mat_bulk,
495 abs_vec_bulk,
496 ptype_bulk,
497 ext_mat_ssbulk,
498 abs_vec_ssbulk,
499 ptype_ssbulk);
500
501 Index f_this = 0;
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);
510 }
511}
512
513void get_dtauc_ssalb(MatrixView dtauc,
514 MatrixView ssalb,
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();
521
522 ARTS_ASSERT(dtauc.nrows() == nf);
523 ARTS_ASSERT(ssalb.nrows() == nf);
524 ARTS_ASSERT(dtauc.ncols() == Np - 1);
525 ARTS_ASSERT(ssalb.ncols() == Np - 1);
526
527 // Initialization
528 dtauc = 0.;
529 ssalb = 0.;
530
531 for (Index ip = 0; ip < Np - 1; ip++)
532 // Do layer averaging and derive single scattering albedo & optical depth
533 for (Index f_index = 0; f_index < nf; f_index++) {
534 Numeric ext =
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));
537 if (ext != 0) {
538 Numeric abs =
539 0.5 *
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;
543
544 ARTS_USER_ERROR_IF((ext - abs) / ext > 1,
545 "ssalb > 1 @ \n",
546 "pressure level = ", ip,"\n",
547 "frequency number = ", f_index,"\n");
548
549 }
550
551 dtauc(f_index, Np - 2 - ip) = ext * (z_profile[ip + 1] - z_profile[ip]);
552 }
553}
554
555void get_angs(Vector& pfct_angs,
556 const ArrayOfArrayOfSingleScatteringData& scat_data,
557 const Index& Npfct) {
558 const Index min_nang = 3;
559 Index nang = Npfct;
560
561 if (Npfct < 0) {
562 Index this_ss = 0, this_se = 0;
563 // determine nang and pfct_angs from scat_data with finest za_grid
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--)
566 // considering scat elems within one species mostly sorted from small to
567 // large sizes with large sizes corresponding to large za_grid. that is,
568 // starting searching from back should trigger if-clause (and variable
569 // update) less often.
570 if (nang < scat_data[i_ss][i_se].za_grid.nelem()) {
571 nang = scat_data[i_ss][i_se].za_grid.nelem();
572 this_ss = i_ss;
573 this_se = i_se;
574 }
575 pfct_angs = scat_data[this_ss][this_se].za_grid;
576 } else if (Npfct < min_nang) {
577 ostringstream os;
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());
582 } else {
583 nlinspace(pfct_angs, 0, 180, nang);
584 }
585}
586
587void get_parZ(Tensor3& pha_bulk_par,
588 const ArrayOfArrayOfSingleScatteringData& scat_data,
589 ConstMatrixView pnd_profiles,
590 ConstVectorView t_profile,
591 ConstVectorView pfct_angs,
592 const ArrayOfIndex& cloudbox_limits) {
593 const Index Np_cloud = pnd_profiles.ncols();
594 const Index nang = pfct_angs.nelem();
595
596 // Initialization
597 pha_bulk_par = 0.;
598
599 // preparing input data
600 Vector T_array{t_profile[Range(cloudbox_limits[0], Np_cloud)]};
601 Matrix idir_array(1, 2, 0.); // we want pfct on sca ang grid, hence set
602 // pdir(*,0) to sca ang, all other to 0.
603 Matrix pdir_array(nang, 2, 0.);
604 pdir_array(joker, 0) = pfct_angs;
605
606 // making particle property output containers
607 ArrayOfArrayOfTensor6 pha_mat_Nse;
608 ArrayOfArrayOfIndex ptypes_Nse;
609 Matrix t_ok;
610 ArrayOfTensor6 pha_mat_ssbulk;
611 ArrayOfIndex ptype_ssbulk;
612 Tensor6 pha_mat_bulk; //nf,nT,npdir,nidir,nst,nst
613 Index ptype_bulk;
614
615 // calculate phase matrix
616 // FIXME: might be optimized by instead just executing pha_mat_1ScatElem where
617 // ext_bulk_par or pnd_field are non-zero.
618 pha_mat_NScatElems(pha_mat_Nse,
619 ptypes_Nse,
620 t_ok,
621 scat_data,
622 1,
623 T_array,
624 pdir_array,
625 idir_array,
626 -1);
627 pha_mat_ScatSpecBulk(pha_mat_ssbulk,
628 ptype_ssbulk,
629 pha_mat_Nse,
630 ptypes_Nse,
631 pnd_profiles,
632 t_ok);
633 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
634
635 pha_bulk_par(joker, Range(cloudbox_limits[0], Np_cloud), joker) =
636 pha_mat_bulk(joker, joker, joker, 0, 0, 0);
637}
638
639void get_pfct(Tensor3& pfct_bulk_par,
640 ConstTensor3View pha_bulk_par,
641 ConstMatrixView ext_bulk_par,
642 ConstMatrixView abs_bulk_par,
643 const ArrayOfIndex& cloudbox_limits) {
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();
648
649 ARTS_ASSERT(pfct_bulk_par.npages() == nf);
650 ARTS_ASSERT(pfct_bulk_par.nrows() == Np - 1);
651 ARTS_ASSERT(pfct_bulk_par.ncols() == nang);
652
653 // Initialization
654 pfct_bulk_par = 0.;
655
656 for (Index ip = cloudbox_limits[0]; ip < Np_cloud - 1; ip++)
657 for (Index f_index = 0; f_index < nf; f_index++) {
658 // Calculate layer averaged scattering (omitting factor 0.5 here as
659 // omitting it for layer averaged Z below, too).
660 Numeric sca =
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));
663 if (sca != 0.) {
664 // Calculate layer averaged Z (omitting factor 0.5) and rescale from
665 // Z (Csca) to P (4Pi)
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;
670 }
671 }
672}
673
674void get_pmom(Tensor3View pmom,
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();
681
682 ARTS_ASSERT(nang == pfct_angs.nelem());
683
684 ARTS_ASSERT(pmom.npages() == nf);
685 ARTS_ASSERT(pmom.nrows() == nlyr);
686 ARTS_ASSERT(pmom.ncols() == Nlegendre);
687
688 Numeric pfct_threshold = 0.1;
689
690 // Initialization
691 pmom = 0.;
692
693 // we need the cosine of the pfct angles
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]);
701 ad_angs[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);
711 }
712 }
713
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)};
719
720 // Check if phase function is properly normalized
721 // For highly peaked phasefunctions, integrating over the angle instead
722 // of over the cosine of the angle is numerically more exact though both
723 // ways are analytically equal. Furthermore in *scat_dataCalc* the
724 // integration is also done over the angle.
725 Numeric pint = 0.;
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.));
730 }
731
732 if (abs(pint / 2. - 1.) > pfct_threshold) {
733 ostringstream os;
734 os << "Phase function normalization deviates from expected value by\n"
735 << 1e2 * pint / 2. - 1e2 << "(allowed: " << pfct_threshold * 1e2
736 << "%).\n"
737 << "Occurs at layer #" << il << " and frequency #" << f_index
738 << ".\n"
739 << "Something is wrong with your scattering data. Check!\n";
740 throw runtime_error(os.str());
741 }
742
743 // for the rest, rescale pfct to norm 2
744 pfct *= 2. / pint;
745
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) +=
750 0.25 * ad_angs[ia] *
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.));
754 }
755 }
756 }
757 }
758}
759
760void get_scat_bulk_layer(MatrixView sca_bulk_layer,
761 const MatrixView& ext_bulk,
762 const MatrixView& abs_bulk) {
763 const Index nf = ext_bulk.nrows();
764 const Index Np = ext_bulk.ncols();
765
766 ARTS_ASSERT(sca_bulk_layer.nrows() == nf);
767 ARTS_ASSERT(sca_bulk_layer.ncols() == Np - 1);
768 ARTS_ASSERT(abs_bulk.nrows() == nf);
769 ARTS_ASSERT(abs_bulk.ncols() == Np);
770
771 // Initialization
772 sca_bulk_layer = 0.;
773
774 for (Index ip = 0; ip < Np - 1; ip++)
775 // Do layer averaging and derive single scattering albedo & optical depth
776 for (Index f_index = 0; f_index < nf; f_index++) {
777 Numeric sca =
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));
780
781 sca_bulk_layer(f_index, Np - 2 - ip) = sca;
782 }
783}
784
785// Use a thread_local variable to communicate the Verbosity to the
786// Disort error and warning functions. Ugly workaround, to avoid
787// passing a Verbosity argument throughout the whole cdisort code.
788// We want to avoid changes to the original code to keep it maintainable
789// in respect to upstream updates.
791
792#define MAX_WARNINGS 100
793
795void c_errmsg(const char* messag, int type) {
796 #if not ARTS_LGPL
797 Verbosity verbosity = disort_verbosity;
798 static int warning_limit = FALSE, num_warnings = 0;
799
800 ARTS_USER_ERROR_IF(type == DS_ERROR, "DISORT ERROR >>> ", messag);
801
802 if (warning_limit) return;
803
804 if (++num_warnings <= MAX_WARNINGS) {
806 out1 << "DISORT WARNING >>> " << messag << "\n";
807 } else {
809 out1 << "DISORT TOO MANY WARNING MESSAGES -- They will no longer be "
810 "printed <<<<<<<\n\n";
811 warning_limit = TRUE;
812 }
813
814 #else
815 ARTS_USER_ERROR("Did not compile with -DENABLE_ARTS_LGPL=0")
816 #endif
817
818 return;
819}
820
821#undef MAX_WARNINGS
822
824int c_write_bad_var(int quiet, const char* varnam) {
825 #if not ARTS_LGPL
826 const int maxmsg = 50;
827 static int nummsg = 0;
828
829 nummsg++;
830 if (quiet != QUIET) {
831 Verbosity verbosity = disort_verbosity;
833 out1 << " **** Input variable " << varnam << " in error ****\n";
834 if (nummsg == maxmsg) {
835 c_errmsg("Too many input errors. Aborting...", DS_ERROR);
836 }
837 }
838
839 return TRUE;
840}
841
843int c_write_too_small_dim(int quiet, const char* dimnam, int minval) {
844 if (quiet != QUIET) {
845 Verbosity verbosity = disort_verbosity;
847 out1 << " **** Symbolic dimension " << dimnam
848 << " should be increased to at least " << minval << " ****\n";
849 }
850
851 return TRUE;
852
853 #else
854 ARTS_USER_ERROR("Did not compile with -DENABLE_ARTS_LGPL=0")
855 #endif
856}
857
858void reduced_1datm(Vector& p,
859 Vector& z,
860 Vector& t,
861 Matrix& vmr,
862 Matrix& pnd,
863 ArrayOfIndex& cboxlims,
864 Index& ncboxremoved,
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,
871 const ArrayOfIndex& cloudbox_limits) {
872 // Surface at p_grid[0] and we just need to copy the original data
873 if (abs(z_surface - z_profile[0]) < 1e-3) {
874 p = p_grid;
875 z = z_profile;
876 t = t_profile;
877 vmr = vmr_profiles;
878 pnd = pnd_profiles;
879 cboxlims = cloudbox_limits;
880 ncboxremoved = 0;
881 }
882 // Surface above p_grid[0]
883 else {
884 // Some counters
885 Index np = p_grid.nelem(), ifirst = 0;
886 // Determine where to start with respect to z_profile
887 for (; z_surface >= z_profile[ifirst + 1]; ifirst++) {
888 }
889 np -= ifirst;
890 // Start by copying from ifirst to end
891 Range ind(ifirst, np);
892 p = p_grid[ind];
893 z = z_profile[ind];
894 t = t_profile[ind];
895 vmr = vmr_profiles(joker, ind);
896 // Insert surface altitude
897 z[0] = z_surface;
898 // Prepare interpolation
899 ArrayOfGridPos gp(1);
900 gridpos(gp[0], z_profile, z_surface);
901 Vector itw(2);
902 interpweights(itw, gp[0]);
903 // t and vmr
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]);
907 }
908 // p (we need a matrix version of iwt to use the function *itw2p*)
909 Matrix itw2(1, 2);
910 itw2(0, 0) = itw[0];
911 itw2(0, 1) = itw[1];
912 itw2p(ExhaustiveVectorView{p[0]}, p, gp, itw2);
913 // pnd_field and cloudbox limits need special treatment
914 cboxlims = cloudbox_limits;
915 if (ifirst < cloudbox_limits[0]) { // Surface below cloudbox
916 cboxlims[0] -= ifirst;
917 cboxlims[1] -= ifirst;
918 pnd = pnd_profiles;
919 ncboxremoved = 0;
920 } else { // Surface inside cloudbox
921 ncboxremoved = ifirst - cboxlims[0];
922 cboxlims[0] = 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]);
929 }
930 }
931 }
932}
933
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,
944 const ArrayOfArrayOfSingleScatteringData& scat_data,
945 const ArrayOfSun& suns,
946 const Agenda& propmat_clearsky_agenda,
947 const Agenda& gas_scattering_agenda,
948 const ArrayOfIndex& cloudbox_limits,
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,
956 const ArrayOfString& disort_aux_vars,
957 const Numeric& scale_factor,
958 const Index& nstreams,
959 const Index& Npfct,
960 const Index& only_tro,
961 const Index& quiet,
962 const Index& emission,
963 const Index& intensity_correction,
964 const Verbosity& verbosity) {
965 // Create an atmosphere starting at z_surface
966 Vector p, z, t;
967 Matrix vmr, pnd;
968 ArrayOfIndex cboxlims;
969 Index ncboxremoved;
970 //
972 z,
973 t,
974 vmr,
975 pnd,
976 cboxlims,
977 ncboxremoved,
978 p_grid,
979 z_profile,
980 z_surface,
981 t_profile,
982 vmr_profiles,
983 pnd_profiles,
984 cloudbox_limits);
985
986 #if not ARTS_LGPL
987 disort_state ds;
988 disort_output out;
989
990 if (quiet == 0)
991 disort_verbosity = verbosity;
992 else
993 disort_verbosity = Verbosity(0, 0, 0);
994
995 const Index nf = f_grid.nelem();
996
997 // solar dependent properties if no sun is present
998 // Number of azimuth angles
999 Index nphi = 1;
1000 //local zenith angle of sun
1001 Numeric umu0 = 0.;
1002 //local azimuth angle of sun
1003 Numeric phi0 = 0.;
1004 //Intensity of incident sun beam
1005 Numeric fbeam = 0.;
1006
1007 if (suns_do) {
1008 nphi = aa_grid.nelem();
1009 umu0 = Conversion::cosd(sun_rte_los[0]);
1010 phi0 = sun_rte_los[1];
1011 if (phi0 < 0) {
1012 phi0 = phi0 + 360.;
1013 }
1014 }
1015
1016 ds.accur = 0.005;
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;
1022
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;
1028
1029 ds.nlyr = static_cast<int>(p.nelem() - 1);
1030
1031 ds.flag.brdf_type = BRDF_NONE;
1032
1033 ds.flag.ibcnd = GENERAL_BC;
1034 ds.flag.usrang = TRUE;
1035
1036 if (emission) {
1037 ds.flag.planck = TRUE;
1038 } else {
1039 ds.flag.planck = FALSE;
1040 }
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;
1047 } else {
1048 ds.flag.intensity_correction = FALSE;
1049 ds.flag.old_intensity_correction = FALSE;
1050 }
1051
1052 ds.nstr = static_cast<int>(nstreams);
1053 ds.nphase = ds.nstr;
1054 ds.nmom = ds.nstr;
1055 //ds.ntau = ds.nlyr + 1; // With ds.flag.usrtau = FALSE; set by cdisort
1056 ds.numu = static_cast<int>(za_grid.nelem());
1057 ds.nphi = static_cast<int>(nphi);
1058 Index Nlegendre = nstreams + 1;
1059
1060 /* Allocate memory */
1061 c_disort_state_alloc(&ds);
1062 c_disort_out_alloc(&ds, &out);
1063
1064 // Looking direction of solar beam
1065 ds.bc.umu0 = umu0;
1066 ds.bc.phi0 = phi0;
1067
1068 // Intensity of bottom-boundary isotropic illumination
1069 ds.bc.fluor = 0.;
1070
1071 // fill up azimuth angle and temperature array
1072 for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i];
1073
1074 if (ds.flag.planck==TRUE){
1075 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1076 }
1077
1078 // Transform to mu, starting with negative values
1079 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);
1080
1081 //gas absorption
1082 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1083 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1084
1085 // Get particle bulk properties
1086 Index nang;
1087 Vector pfct_angs;
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;
1091
1092 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1093 nang = Npfct;
1094 nlinspace(pfct_angs, 0, 180, nang);
1095
1096 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1097
1098 ext_bulk_par = 0.0;
1099 abs_bulk_par = 0.0;
1100 pha_bulk_par = 0.0;
1101
1102 Index iflat = 0;
1103
1104 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1105 const Index nse = scat_data[iss].nelem();
1106 ext_abs_pfun_from_tro(ext_bulk_par,
1107 abs_bulk_par,
1108 pha_bulk_par,
1109 scat_data[iss],
1110 iss,
1111 pnd(Range(iflat, nse), joker),
1112 cboxlims,
1113 t,
1114 pfct_angs);
1115 iflat += nse;
1116 }
1117 } else {
1118 get_angs(pfct_angs, scat_data, Npfct);
1119 nang = pfct_angs.nelem();
1120
1121 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1122
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);
1126 }
1127
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);
1130
1131 // Legendre polynomials of phase function
1132 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1133 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1134
1135 if (gas_scattering_do) {
1136 // gas scattering
1137
1138 // layer averaged particle scattering coefficient
1139 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1140 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1141
1142 // call gas_scattering_properties
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.);
1146
1148 sca_coeff_gas_layer,
1149 sca_coeff_gas_level,
1150 pmom_gas,
1151 f_grid,
1152 p,
1153 t,
1154 vmr,
1155 gas_scattering_agenda);
1156
1157 // call add_norm_phase_functions
1159 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1160
1161 // add gas_scat_ext to ext_bulk_par
1162 ext_bulk_par += sca_coeff_gas_level;
1163 }
1164
1165 // Optical depth of layers
1166 Matrix dtauc(nf, ds.nlyr);
1167 // Single scattering albedo of layers
1168 Matrix ssalb(nf, ds.nlyr);
1169 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1170
1171 //upper boundary conditions:
1172 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1173 // emission. Both are applied additively.
1174 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1175 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1176 // are additive we should use either the one or the other.
1177 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1178 // Generally, temis!=1 should be avoided since that technically implies a
1179 // reflective upper boundary (though it seems that this is not exploited in
1180 // DISORT1.2, which we so far use).
1181
1182 // Cosmic background
1183 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1184 ds.bc.fisot = 0;
1185
1186 // Top of the atmosphere temperature and emissivity
1187 ds.bc.ttemp = COSMIC_BG_TEMP;
1188 ds.bc.btemp = surface_skin_t;
1189 ds.bc.temis = 1.;
1190
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);
1193
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);
1200
1201 // Wavenumber in [1/cm]
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;
1205
1206 // set
1207 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1208
1209 // Set irradiance of incident solar beam at top boundary
1210 if (suns_do) {
1211 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1212 (100 * SPEED_OF_LIGHT)*scale_factor;
1213 }
1214 ds.bc.fbeam = fbeam;
1215
1216 std::memcpy(ds.pmom,
1217 pmom(f_index, joker, joker).unsafe_data_handle(),
1218 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1219
1220 enum class Status { FIRST_TRY, RETRY, SUCCESS };
1221 Status tries = Status::FIRST_TRY;
1222 const Numeric eps = 2e-4; //two times the value defined in cdisort.c:3653
1223 do {
1224 try {
1225 c_disort(&ds, &out);
1226 tries = Status::SUCCESS;
1227 } catch (const std::runtime_error& e) {
1228 //catch cases if solar zenith angle=quadrature angle
1229 if (tries == Status::FIRST_TRY) {
1230 // change angle
1231 if (umu0 < 1 - eps) {
1232 umu0 += eps;
1233 } else if (umu0 > 1 - eps) {
1234 umu0 -= eps;
1235 }
1236
1237 const Numeric shift =
1238 abs(Conversion::acosd(umu0) - Conversion::acosd(ds.bc.umu0));
1240 out1
1241 << "Solar zenith angle coincided with one of the quadrature angles\n"
1242 << "We needed to shift the solar sun angle by " << shift
1243 << "deg.\n";
1244
1245 ds.bc.umu0 = umu0;
1246 tries = Status::RETRY;
1247 } else
1248 throw e;
1249 }
1250 } while (tries != Status::SUCCESS);
1251
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)) *
1257 ds.numu] /
1258 (ds.wvnmhi - ds.wvnmlo) / (100 * SPEED_OF_LIGHT);
1259 }
1260 // To avoid potential numerical problems at interpolation of the field,
1261 // we copy the surface field to underground altitudes
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);
1265 }
1266 }
1267 }
1268 }
1269
1270 // Allocate aux data
1271 disort_aux.resize(disort_aux_vars.nelem());
1272 // Allocate and set (if possible here) iy_aux
1273 Index cnt=-1;
1274 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1275
1276
1277 if (disort_aux_vars[i] == "Layer optical thickness"){
1278 cnt += 1;
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);
1284 }
1285 }
1286 disort_aux[cnt] = deltatau;
1287 }
1288 else if (disort_aux_vars[i] == "Single scattering albedo"){
1289 cnt+=1;
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);
1295 }
1296 }
1297 disort_aux[cnt]=snglsctalbedo;
1298 }
1299 else if (disort_aux_vars[i] == "Direct beam") {
1300 cnt += 1;
1301 Matrix directbeam(nf, ds.nlyr + 1, 0);
1302 if (suns_do) {
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;
1306
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);
1311 }
1312 }
1313 }
1314 disort_aux[cnt]=directbeam;
1315 } else {
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");
1322 }
1323 }
1324
1325 /* Free allocated memory */
1326 c_disort_out_free(&ds, &out);
1327 c_disort_state_free(&ds);
1328
1329 #else
1330 ARTS_USER_ERROR("Did not compile with -DENABLE_ARTS_LGPL=0")
1331 #endif
1332}
1333
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,
1344 const ArrayOfArrayOfSingleScatteringData& scat_data,
1345 const ArrayOfSun& suns,
1346 const Agenda& propmat_clearsky_agenda,
1347 const Agenda& gas_scattering_agenda,
1348 const ArrayOfIndex& cloudbox_limits,
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,
1354 const ArrayOfString& disort_aux_vars,
1355 const Numeric& scale_factor,
1356 const Index& nstreams,
1357 const Index& Npfct,
1358 const Index& only_tro,
1359 const Index& quiet,
1360 const Index& emission,
1361 const Index& intensity_correction,
1362 const Verbosity& verbosity) {
1363 // Create an atmosphere starting at z_surface
1364 Vector p, z, t;
1365 Matrix vmr, pnd;
1366 ArrayOfIndex cboxlims;
1367 Index ncboxremoved;
1368 //
1369 reduced_1datm(p,
1370 z,
1371 t,
1372 vmr,
1373 pnd,
1374 cboxlims,
1375 ncboxremoved,
1376 p_grid,
1377 z_profile,
1378 z_surface,
1379 t_profile,
1380 vmr_profiles,
1381 pnd_profiles,
1382 cloudbox_limits);
1383
1384#if not ARTS_LGPL
1385
1386 disort_state ds;
1387 disort_output out;
1388
1389 if (quiet == 0)
1390 disort_verbosity = verbosity;
1391 else
1392 disort_verbosity = Verbosity(0, 0, 0);
1393
1394 const Index nf = f_grid.nelem();
1395
1396 // solar dependent properties if no sun is present
1397 //local zenith angle of sun
1398 Numeric umu0 = 0.;
1399 //local azimuth angle of sun
1400 Numeric phi0 = 0.;
1401 //Intensity of incident sun beam
1402 Numeric fbeam = 0.;
1403
1404 if (suns_do) {
1405 umu0 = Conversion::cosd(sun_rte_los[0]);
1406 phi0 = sun_rte_los[1];
1407 if (phi0 < 0) {
1408 phi0 = phi0 + 360.;
1409 }
1410 }
1411
1412 ds.accur = 0.005;
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;
1418
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;
1424
1425 ds.nlyr = static_cast<int>(p.nelem() - 1);
1426
1427 ds.flag.brdf_type = BRDF_NONE;
1428
1429 ds.flag.ibcnd = GENERAL_BC;
1430 ds.flag.usrang = FALSE;
1431
1432 if (emission) {
1433 ds.flag.planck = TRUE;
1434 } else {
1435 ds.flag.planck = FALSE;
1436 }
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;
1443 } else {
1444 ds.flag.intensity_correction = FALSE;
1445 ds.flag.old_intensity_correction = FALSE;
1446 }
1447
1448 ds.nstr = static_cast<int>(nstreams);
1449 ds.nphase = ds.nstr;
1450 ds.nmom = ds.nstr;
1451
1452 // set grid size of angular variable to dummy value
1453 ds.numu = 1;
1454 ds.nphi = 1;
1455
1456 //Set number of legendre terms
1457 Index Nlegendre = nstreams + 1;
1458
1459 /* Allocate memory */
1460 c_disort_state_alloc(&ds);
1461 c_disort_out_alloc(&ds, &out);
1462
1463 // Looking direction of solar beam
1464 ds.bc.umu0 = umu0;
1465 ds.bc.phi0 = phi0;
1466
1467 // Intensity of bottom-boundary isotropic illumination
1468 ds.bc.fluor = 0.;
1469
1470 // fill up temperature array
1471 if (ds.flag.planck==TRUE){
1472 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1473 }
1474
1475 //gas absorption
1476 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1477 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1478
1479 // Get particle bulk properties
1480 Index nang;
1481 Vector pfct_angs;
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;
1485
1486 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1487 nang = Npfct;
1488 nlinspace(pfct_angs, 0, 180, nang);
1489
1490 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1491
1492 ext_bulk_par = 0.0;
1493 abs_bulk_par = 0.0;
1494 pha_bulk_par = 0.0;
1495
1496 Index iflat = 0;
1497
1498 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1499 const Index nse = scat_data[iss].nelem();
1500 ext_abs_pfun_from_tro(ext_bulk_par,
1501 abs_bulk_par,
1502 pha_bulk_par,
1503 scat_data[iss],
1504 iss,
1505 pnd(Range(iflat, nse), joker),
1506 cboxlims,
1507 t,
1508 pfct_angs);
1509 iflat += nse;
1510 }
1511 } else {
1512 get_angs(pfct_angs, scat_data, Npfct);
1513 nang = pfct_angs.nelem();
1514
1515 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1516
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);
1520 }
1521
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);
1524
1525 // Legendre polynomials of phase function
1526 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1527 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1528
1529 if (gas_scattering_do) {
1530 // gas scattering
1531
1532 // layer averaged particle scattering coefficient
1533 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1534 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1535
1536 // call gas_scattering_properties
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.);
1540
1542 sca_coeff_gas_layer,
1543 sca_coeff_gas_level,
1544 pmom_gas,
1545 f_grid,
1546 p,
1547 t,
1548 vmr,
1549 gas_scattering_agenda);
1550
1551 // call add_norm_phase_functions
1553 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1554
1555 // add gas_scat_ext to ext_bulk_par
1556 ext_bulk_par += sca_coeff_gas_level;
1557 }
1558
1559 // Optical depth of layers
1560 Matrix dtauc(nf, ds.nlyr);
1561 // Single scattering albedo of layers
1562 Matrix ssalb(nf, ds.nlyr);
1563 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1564
1565 //upper boundary conditions:
1566 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1567 // emission. Both are applied additively.
1568 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1569 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1570 // are additive we should use either the one or the other.
1571 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1572 // Generally, temis!=1 should be avoided since that technically implies a
1573 // reflective upper boundary (though it seems that this is not exploited in
1574 // DISORT1.2, which we so far use).
1575
1576 // Cosmic background
1577 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1578 ds.bc.fisot = 0;
1579
1580 // Top of the atmosphere temperature and emissivity
1581 ds.bc.ttemp = COSMIC_BG_TEMP;
1582 ds.bc.btemp = surface_skin_t;
1583 ds.bc.temis = 1.;
1584
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);
1589
1590 if (suns_do){
1591 //Resize direct field
1592 spectral_direct_irradiance_field.resize(nf, ds.nlyr+1);
1593 spectral_direct_irradiance_field = 0;
1594 }
1595
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);
1598
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);
1605
1606 // Wavenumber in [1/cm]
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;
1610
1611 // set
1612 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1613
1614 // Set irradiance of incident solar beam at top boundary
1615 if (suns_do) {
1616 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1617 (100 * SPEED_OF_LIGHT)*scale_factor;
1618 }
1619 ds.bc.fbeam = fbeam;
1620
1621 std::memcpy(ds.pmom,
1622 pmom(f_index, joker, joker).unsafe_data_handle(),
1623 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1624
1625 c_disort(&ds, &out);
1626
1627 //factor for converting it into spectral radiance units
1628 const Numeric conv_fac=(ds.wvnmhi - ds.wvnmlo) * (100 * SPEED_OF_LIGHT);
1629
1630 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1631 if (suns_do){
1632 // downward direct flux
1633 spectral_direct_irradiance_field(f_index, k + ncboxremoved) =
1634 -out.rad[ds.nlyr - k - cboxlims[0]].rfldir/conv_fac;
1635
1636 // downward total flux
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;
1640
1641 } else {
1642 // downward total flux
1643 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1644 -out.rad[ds.nlyr - k - cboxlims[0]].rfldn/conv_fac;
1645 }
1646
1647 // upward flux
1648 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 1) =
1649 out.rad[ds.nlyr - k - cboxlims[0]].flup/conv_fac;
1650
1651 // flux divergence in tau space
1652 dFdtau(f_index, k + ncboxremoved) =
1653 -out.rad[ds.nlyr - k - cboxlims[0]].dfdt;
1654
1655 // k is running over the number of levels but deltatau, ssalb is defined for layers,
1656 // therefore we need to exlude k==0 and remove one from the index.
1657 if (k>0){
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]];
1660 }
1661 }
1662
1663 // To avoid potential numerical problems at interpolation of the field,
1664 // we copy the surface field to underground altitudes
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);
1668
1669 if (suns_do) {
1670 spectral_direct_irradiance_field(f_index, k) =
1671 spectral_direct_irradiance_field(f_index, k + 1);
1672 }
1673 dFdtau(f_index, k) = dFdtau(f_index, k + 1);
1674 }
1675 }
1676
1677 // Allocate aux data
1678 disort_aux.resize(disort_aux_vars.nelem());
1679 // Allocate and set (if possible here) iy_aux
1680 Index cnt=-1;
1681 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1682
1683
1684 if (disort_aux_vars[i] == "Layer optical thickness"){
1685 cnt+=1;
1686 disort_aux[cnt]=deltatau;
1687 }
1688 else if (disort_aux_vars[i] == "Single scattering albedo"){
1689 cnt+=1;
1690 disort_aux[cnt]=snglsctalbedo;
1691 }
1692 else if (disort_aux_vars[i] == "Direct downward spectral irradiance") {
1693 cnt += 1;
1694 disort_aux[cnt] = spectral_direct_irradiance_field;
1695 }
1696 else if (disort_aux_vars[i] == "dFdtau") {
1697 cnt += 1;
1698 disort_aux[cnt] = dFdtau;
1699 }
1700 else {
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"
1706 " \"dFdtau\"\n"
1707 "but you have selected: \"", disort_aux_vars[i], "\"\n");
1708 }
1709 }
1710
1711
1712
1713 /* Free allocated memory */
1714 c_disort_out_free(&ds, &out);
1715 c_disort_state_free(&ds);
1716
1717 #else
1718 ARTS_USER_ERROR("Did not compile with -DENABLE_ARTS_LGPL=0")
1719 #endif
1720}
1721
1723 //Output
1724 VectorView albedo,
1725 Numeric& btemp,
1726 //Input
1727 const Agenda& surface_rtprop_agenda,
1728 ConstVectorView f_grid,
1729 ConstVectorView scat_za_grid,
1730 const Numeric& surf_alt,
1731 const Verbosity& verbosity) {
1732 // Here, we derive an average surface albedo of the setup as given by ARTS'
1733 // surface_rtprop_agenda to use with Disorts's proprietary Lambertian surface.
1734 // In this way, ARTS-Disort can approximately mimick all surface reflection
1735 // types that ARTS itself can handle.
1736 // Surface temperature as derived from surface_rtprop_agenda is also returned.
1737 //
1738 // We derive the reflection matrices over all incident and reflected polar
1739 // angle directions and derive their integrated value (or weighted average).
1740 // The surface_rtprop_agenda handles one reflected direction at a time
1741 // (rtp_los) and for the reflected directions we loop over all (upwelling)
1742 // angles as given by scat_za_grid. The derived reflection matrices already
1743 // represent the reflectivity (or power reflection coefficient) for the given
1744 // reflected direction including proper angle weighting. For integrating/
1745 // averaging over the reflected directions, we use the same approach, i.e.
1746 // weight each angle by its associated range as given by the half distances to
1747 // the neighboring grid angles (using ARTS' scat_za_grid means 0/180deg are
1748 // grid points, 90deg shouldn't be among them (resulting from even number
1749 // requirement for Disort and the (implicitly assumed?) requirement of a
1750 // horizon-symmetric za_grid)).
1751 //
1752 // We do all frequencies here at once (assuming this is the faster variant as
1753 // the agenda anyway (can) provide output for full f_grid at once and as we
1754 // have to apply the same inter/extrapolation to all the frequencies).
1755
1757
1758 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1759
1760 const Index nf = f_grid.nelem();
1761 Index frza = 0;
1762 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
1763 if (frza == scat_za_grid.nelem()) {
1764 ostringstream os;
1765 os << "No upwelling direction found in scat_za_grid.\n";
1766 throw runtime_error(os.str());
1767 }
1768 const Index nrza = scat_za_grid.nelem() - frza;
1769 //cout << nrza << " upwelling directions found, starting from element #"
1770 // << frza << " of scat_za_grid.\n";
1771 Matrix dir_refl_coeff(nrza, nf, 0.);
1772
1773 // Local input of surface_rtprop_agenda.
1774 Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
1775
1776 // first derive the (reflected-)direction dependent power reflection
1777 // coefficient
1778 for (Index rza = 0; rza < nrza; rza++) {
1779 // Local output of surface_rtprop_agenda.
1780 Numeric surface_skin_t;
1781 Matrix surface_los;
1782 Tensor4 surface_rmatrix;
1783 Matrix surface_emission;
1784
1785 Vector rtp_los(1, scat_za_grid[rza + frza]);
1786 out2 << "Doing reflected dir #" << rza << " at " << rtp_los[0] << " degs\n";
1787
1789 surface_skin_t,
1790 surface_emission,
1791 surface_los,
1792 surface_rmatrix,
1793 Vector{f_grid},
1794 rtp_pos,
1795 rtp_los,
1796 surface_rtprop_agenda);
1797 if (surface_los.nrows() > 1) {
1798 ostringstream os;
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());
1802 }
1803
1804 if (rza == 0)
1805 btemp = surface_skin_t;
1806 else if (surface_skin_t != btemp) {
1807 ostringstream os;
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());
1812 }
1813 // Nothing to do if no surface_los (which means blackbody surface)
1814 // as dir_refl_coeff already set to 0
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));
1819 }
1820 out2 << " directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1821 << "\n";
1822 }
1823
1824 if (btemp < 0. || btemp > 1000.) {
1825 ostringstream os;
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());
1829 }
1830
1831 // now integrate/average the (reflected-)direction dependent power reflection
1832 // coefficients
1833 //
1834 // starting with deriving the angles defining the angle ranges
1835 Vector surf_int_grid(nrza + 1);
1836 // the first angle grid point should be around (but above) 90deg and should
1837 // cover the angle range between the 90deg and half-way point towards the next
1838 // angle grid point. we probably also want to check, that we don't
1839 // 'extrapolate' too much.
1840 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1841 ostringstream os;
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());
1845 }
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) {
1850 ostringstream os;
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());
1856 }
1857 if (!is_same_within_epsilon(
1858 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1859 ostringstream os;
1860 os << "Looks like last point in scat_za_grid is not 180deg.\n";
1861 throw runtime_error(os.str());
1862 }
1863
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]);
1869 surf_int_grid *= DEG2RAD;
1870
1871 // now calculating the actual weights and apply them
1872 for (Index rza = 0; rza < nrza; rza++) {
1873 //Numeric coslow = cos(2.*surf_int_grid[rza]);
1874 //Numeric coshigh = cos(2.*surf_int_grid[rza+1]);
1875 //Numeric w = 0.5*(coshigh-coslow);
1876 Numeric w =
1877 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1878 //cout << "at reflLOS[" << rza << "]=" << scat_za_grid[frza+rza] << ":\n";
1879 //cout << " angle weight derives as w = 0.5*(" << coshigh << "-"
1880 // << coslow << ") = " << w << "\n";
1881 //cout << " weighting directional reflection coefficient from "
1882 // << dir_refl_coeff(rza,joker);
1883 dir_refl_coeff(rza, joker) *= w;
1884 //cout << " to " << dir_refl_coeff(rza,joker) << "\n";
1885 }
1886
1887 // eventually sum up the weighted directional power reflection coefficients
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.) {
1893 ostringstream os;
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());
1898 }
1899 }
1900}
1901
1903 //Output
1904 VectorView albedo,
1905 Numeric& btemp,
1906 //Input
1907 const Agenda& surface_rtprop_agenda,
1908 ConstVectorView f_grid,
1909 const Numeric& surf_alt,
1910 const Numeric& inc_angle) {
1911
1912 Vector rtp_pos(1, surf_alt);
1913 Vector rtp_los(1, 180-inc_angle);
1914
1915 Numeric surface_skin_t;
1916 Matrix surface_los;
1917 Tensor4 surface_rmatrix;
1918 Matrix surface_emission;
1919
1921 surface_skin_t,
1922 surface_emission,
1923 surface_los,
1924 surface_rmatrix,
1925 Vector{f_grid},
1926 rtp_pos,
1927 rtp_los,
1928 surface_rtprop_agenda);
1929
1930 if (surface_los.nrows() > 1) {
1931 ostringstream os;
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());
1935 }
1936 if (surface_los.empty()) {
1937 albedo = 0; // Blackbody assumed if no surface_los
1938 } else {
1939 albedo[joker] = surface_rmatrix(0,joker,0,0);
1940 }
1941
1942 btemp = surface_skin_t;
1943 if (btemp < 0. || btemp > 1000.) {
1944 ostringstream os;
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());
1948 }
1949}
1950
1951
1952#if ARTS_LGPL
1953#pragma GCC diagnostic pop
1954#endif
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
Definition: array.h:128
base min(const Array< base > &x)
Min function.
Definition: array.h:144
Constants of physical expressions as constexpr.
Common ARTS conversions.
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)
Definition: auto_md.cc:25049
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)
Definition: auto_md.cc:25794
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)
Definition: auto_md.cc:25978
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:608
The Agenda class.
Definition: agenda_class.h:52
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:75
Workspace class.
Definition: workspace_ng.h:36
Helper macros for debugging.
#define DEBUG_ONLY(...)
Definition: debug.h:53
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
#define ARTS_USER_ERROR(...)
Definition: debug.h:151
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:135
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
Definition: disort.cc:1722
constexpr Numeric PLANCK_CONST
Definition: disort.cc:48
void c_errmsg(const char *messag, int type)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:795
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.
Definition: disort.cc:639
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
Definition: disort.cc:858
constexpr Numeric SPEED_OF_LIGHT
Definition: disort.cc:49
#define MAX_WARNINGS
Definition: disort.cc:792
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.
Definition: disort.cc:934
constexpr Numeric DEG2RAD
Definition: disort.cc:47
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
Definition: disort.cc:375
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.
Definition: disort.cc:587
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.
Definition: disort.cc:1334
void get_disortsurf_props(Vector &albedo, Numeric &btemp, ConstVectorView f_grid, const Numeric &surface_skin_t, ConstVectorView surface_scalar_reflectivity)
get_disortsurf_props.
Definition: disort.cc:286
void get_scat_bulk_layer(MatrixView sca_bulk_layer, const MatrixView &ext_bulk, const MatrixView &abs_bulk)
get_scat_bulk_layer
Definition: disort.cc:760
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
Definition: disort.cc:1902
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.
Definition: disort.cc:202
int c_write_bad_var(int quiet, const char *varnam)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:824
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
Definition: disort.cc:513
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.
Definition: disort.cc:76
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.
Definition: disort.cc:269
thread_local Verbosity disort_verbosity
Definition: disort.cc:790
void get_angs(Vector &pfct_angs, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &Npfct)
get_angs.
Definition: disort.cc:555
void get_pmom(Tensor3View pmom, ConstTensor3View pfct_bulk_par, ConstVectorView pfct_angs, const Index &Nlegendre)
get_pmom
Definition: disort.cc:674
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.
Definition: disort.cc:331
constexpr Numeric COSMIC_BG_TEMP
Definition: disort.cc:50
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)
Definition: disort.cc:438
constexpr Numeric PI
Definition: disort.cc:46
int c_write_too_small_dim(int quiet, const char *dimnam, int minval)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:843
void add_normed_phase_functions(Tensor3View pfct1, const MatrixView &sca1, const MatrixView &pfct2, const MatrixView &sca2)
add_normed_phase_functions
Definition: disort.cc:52
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
Definition: jacobian.h:529
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:208
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:187
#define CREATE_OUT2
Definition: messages.h:188
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.
@ PTYPE_GENERAL
Definition: optproperties.h:36
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:37
@ PTYPE_TOTAL_RND
Definition: optproperties.h:38
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.
#define u
#define w
This file contains basic functions to handle XML data files.