ARTS 2.5.4 (git: 31ce4f0e)
disort.cc
Go to the documentation of this file.
1/* Copyright (C) 2006-2012 Claudia Emde <claudia.emde@dlr.de>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA.
17*/
18
29#include "disort.h"
30
31#include <cinttypes>
32#include <cmath>
33#include <cstdlib>
34#include <iostream>
35#include <stdexcept>
36
37#include "agenda_class.h"
38#include "array.h"
39#include "arts_constants.h"
40#include "auto_md.h"
41#include "check_input.h"
42#include "arts_conversions.h"
43
44extern "C" {
45#include "cdisort.h"
46}
47
48#include "disort.h"
49#include "interpolation.h"
50#include "logic.h"
51#include "math_funcs.h"
52#include "messages.h"
53#include "rte.h"
54#include "xml_io.h"
55
56inline constexpr Numeric PI=Constant::pi;
61
63 const MatrixView& sca1,
64 const MatrixView& pfct2,
65 const MatrixView& sca2) {
66 const Index np1 = pfct1.npages();
67 const Index nr1 = pfct1.nrows();
68 const Index nc1 = pfct1.ncols();
69
70
71 ARTS_ASSERT(pfct2.nrows() == nr1);
72
73 ARTS_ASSERT(pfct2.ncols() == nc1);
74
75 for (Index i = 0; i < np1; i++) { // frequncy loop
76 for (Index j = 0; j < nr1 ; j++) { // layer loop
77 for (Index k = 0; k < nc1; k++) // polynomial loop
78
79 pfct1(i, j, k) =
80 (sca1(i, j) * pfct1(i, j, k) + sca2(i, j) * pfct2(j, k)) /
81 (sca1(i, j) + sca2(i, j));
82 }
83 }
84}
85
86void check_disort_input( // Input
87 const Index& cloudbox_on,
88 const Index& atmfields_checked,
89 const Index& atmgeom_checked,
90 const Index& cloudbox_checked,
91 const Index& scat_data_checked,
92 const Index& atmosphere_dim,
93 const Index& stokes_dim,
94 const ArrayOfIndex& cloudbox_limits,
96 ConstVectorView za_grid,
97 const Index& nstreams) {
98 if (!cloudbox_on) {
99 throw runtime_error(
100 "Cloudbox is off, no scattering calculations to be"
101 "performed.");
102 }
103
104 if (atmfields_checked != 1)
105 throw runtime_error(
106 "The atmospheric fields must be flagged to have "
107 "passed a consistency check (atmfields_checked=1).");
108 if (atmgeom_checked != 1)
109 throw runtime_error(
110 "The atmospheric geometry must be flagged to have "
111 "passed a consistency check (atmgeom_checked=1).");
112 if (cloudbox_checked != 1)
113 throw runtime_error(
114 "The cloudbox must be flagged to have "
115 "passed a consistency check (cloudbox_checked=1).");
116 if (scat_data_checked != 1)
117 throw runtime_error(
118 "The scat_data must be flagged to have "
119 "passed a consistency check (scat_data_checked=1).");
120
121 if (atmosphere_dim != 1)
122 throw runtime_error(
123 "For running DISORT, atmospheric dimensionality "
124 "must be 1.\n");
125
126 if (stokes_dim < 0 || stokes_dim > 1)
127 throw runtime_error(
128 "For running DISORT, the dimension of stokes vector "
129 "must be 1.\n");
130
131 if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
132 throw runtime_error(
133 "*cloudbox_limits* is a vector which contains the"
134 "upper and lower limit of the cloud for all "
135 "atmospheric dimensions. So its dimension must"
136 "be 2 x *atmosphere_dim*");
137
138 if (cloudbox_limits[0] != 0) {
139 ostringstream os;
140 os << "DISORT calculations currently only possible with "
141 << "lower cloudbox limit\n"
142 << "at 0th atmospheric level "
143 << "(assumes surface there, ignoring z_surface).\n";
144 throw runtime_error(os.str());
145 }
146
147 if (scat_data.empty())
148 throw runtime_error(
149 "No single scattering data present.\n"
150 "See documentation of WSV *scat_data* for options to "
151 "make single scattering data available.\n");
152
153 // DISORT requires even number of streams:
154 // nstreams is total number of directions, up- and downwelling, and the up-
155 // and downwelling directions need to be symmetrically distributed, i.e. same
156 // number of directions in both hemispheres is required. horizontal direction
157 // (90deg) can not be covered in a plane-parallel atmosphere.
158 if (nstreams / 2 * 2 != nstreams) {
159 ostringstream os;
160 os << "DISORT requires an even number of streams, but yours is " << nstreams
161 << ".\n";
162 throw runtime_error(os.str());
163 }
164
165 // Zenith angle grid.
166 Index nza = za_grid.nelem();
167
168 // za_grid here is only relevant to provide an i_field from which the
169 // sensor los angles can be interpolated by yCalc; it does not the determine
170 // the accuracy of the DISORT output itself at these angles. So we can only
171 // apply a very rough test here, whether the grid is appropriate. However, we
172 // set the threshold fairly high since calculation costs for a higher number
173 // of angles are negligible.
174 if (nza < 20) {
175 ostringstream os;
176 os << "We require size of za_grid to be >= 20, to ensure a\n"
177 << "reasonable interpolation of the calculated cloudbox field.\n"
178 << "Note that for DISORT additional computation costs for\n"
179 << "larger numbers of angles are negligible.";
180 throw runtime_error(os.str());
181 }
182
183 if (za_grid[0] != 0. || za_grid[nza - 1] != 180.)
184 throw runtime_error("The range of *za_grid* must [0 180].");
185
186 if (!is_increasing(za_grid))
187 throw runtime_error("*za_grid* must be increasing.");
188
189 Index i = 1;
190 while (za_grid[i] <= 90) {
191 if (za_grid[i] == 90)
192 throw runtime_error("*za_grid* is not allowed to contain the value 90");
193 i++;
194 }
195
196 // DISORT can only handle randomly oriented particles.
197 bool all_totrand = true;
198 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
199 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++)
200 if (scat_data[i_ss][i_se].ptype != PTYPE_TOTAL_RND) all_totrand = false;
201 if (!all_totrand) {
202 ostringstream os;
203 os << "DISORT can only handle scattering elements of type "
204 << PTYPE_TOTAL_RND << " (" << PTypeToString(PTYPE_TOTAL_RND) << "),\n"
205 << "but at least one element of other type (" << PTYPE_AZIMUTH_RND << "="
206 << PTypeToString(PTYPE_AZIMUTH_RND) << " or " << PTYPE_GENERAL << "="
207 << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
208 throw runtime_error(os.str());
209 }
210}
211
213 const Index& atmfields_checked,
214 const Index& atmgeom_checked,
215 const Index& scat_data_checked,
216 const Index& atmosphere_dim,
217 const Index& stokes_dim,
218 const ArrayOfArrayOfSingleScatteringData& scat_data,
219 const Index& nstreams) {
220 if (atmfields_checked != 1)
221 throw runtime_error(
222 "The atmospheric fields must be flagged to have "
223 "passed a consistency check (atmfields_checked=1).");
224
225 if (atmgeom_checked != 1)
226 throw runtime_error(
227 "The atmospheric geometry must be flagged to have "
228 "passed a consistency check (atmgeom_checked=1).");
229
230 if (scat_data_checked != 1)
231 throw runtime_error(
232 "The scat_data must be flagged to have "
233 "passed a consistency check (scat_data_checked=1).");
234
235 if (atmosphere_dim != 1)
236 throw runtime_error(
237 "For running DISORT, atmospheric dimensionality "
238 "must be 1.\n");
239
240 if (stokes_dim < 0 || stokes_dim > 1)
241 throw runtime_error(
242 "For running DISORT, the dimension of stokes vector "
243 "must be 1.\n");
244
245 if (scat_data.empty())
246 throw runtime_error(
247 "No single scattering data present.\n"
248 "See documentation of WSV *scat_data* for options to "
249 "make single scattering data available.\n");
250
251 // DISORT requires even number of streams:
252 // nstreams is total number of directions, up- and downwelling, and the up-
253 // and downwelling directions need to be symmetrically distributed, i.e. same
254 // number of directions in both hemispheres is required. horizontal direction
255 // (90deg) can not be covered in a plane-parallel atmosphere.
256 if (nstreams / 2 * 2 != nstreams) {
257 ostringstream os;
258 os << "DISORT requires an even number of streams, but yours is " << nstreams
259 << ".\n";
260 throw runtime_error(os.str());
261 }
262
263 // DISORT can only handle randomly oriented particles.
264 bool all_totrand = true;
265 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
266 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++)
267 if (scat_data[i_ss][i_se].ptype != PTYPE_TOTAL_RND) all_totrand = false;
268 if (!all_totrand) {
269 ostringstream os;
270 os << "DISORT can only handle scattering elements of type "
271 << PTYPE_TOTAL_RND << " (" << PTypeToString(PTYPE_TOTAL_RND) << "),\n"
272 << "but at least one element of other type (" << PTYPE_AZIMUTH_RND << "="
273 << PTypeToString(PTYPE_AZIMUTH_RND) << " or " << PTYPE_GENERAL << "="
274 << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
275 throw runtime_error(os.str());
276 }
277}
278
279void init_ifield( // Output
280 Tensor7& cloudbox_field,
281 // Input
282 const Vector& f_grid,
283 const ArrayOfIndex& cloudbox_limits,
284 const Index& n_za,
285 const Index& n_aa,
286 const Index& stokes_dim) {
287 const Index Nf = f_grid.nelem();
288 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
289 //const Index Nza = za_grid.nelem();
290
291 // Resize and initialize radiation field in the cloudbox
292 cloudbox_field.resize(Nf, Np_cloud, 1, 1, n_za, n_aa, stokes_dim);
293 cloudbox_field = NAN;
294}
295
296void get_disortsurf_props( // Output
297 Vector& albedo,
298 Numeric& btemp,
299 // Input
300 ConstVectorView f_grid,
301 const Numeric& surface_skin_t,
302 ConstVectorView surface_scalar_reflectivity) {
303 // temperature of surface
304 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
305 ostringstream os;
306 os << "Surface temperature has been set or derived as " << btemp << " K,\n"
307 << "which is not considered a meaningful value.\n"
308 << "For surface method 'L', *surface_skin_t* needs to\n"
309 << "be set and passed explicitly. Maybe you didn't do this?";
310 throw runtime_error(os.str());
311 }
312 btemp = surface_skin_t;
313
314 // surface albedo
315 if (surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
316 surface_scalar_reflectivity.nelem() != 1) {
317 ostringstream os;
318 os << "The number of elements in *surface_scalar_reflectivity*\n"
319 << "should match length of *f_grid* or be 1."
320 << "\n length of *f_grid* : " << f_grid.nelem()
321 << "\n length of *surface_scalar_reflectivity* : "
322 << surface_scalar_reflectivity.nelem() << "\n";
323 throw runtime_error(os.str());
324 }
325
326 if (min(surface_scalar_reflectivity) < 0 ||
327 max(surface_scalar_reflectivity) > 1) {
328 throw runtime_error(
329 "All values in *surface_scalar_reflectivity*"
330 " must be inside [0,1].");
331 }
332
333 if (surface_scalar_reflectivity.nelem() > 1)
334 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
335 albedo[f_index] = surface_scalar_reflectivity[f_index];
336 else
337 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
338 albedo[f_index] = surface_scalar_reflectivity[0];
339}
340
342 MatrixView ext_bulk_gas,
343 const Agenda& propmat_clearsky_agenda,
344 ConstVectorView t_profile,
345 ConstMatrixView vmr_profiles,
346 ConstVectorView p_grid,
347 ConstVectorView f_grid) {
348 const Index Np = p_grid.nelem();
349
350 ARTS_ASSERT(ext_bulk_gas.nrows() == f_grid.nelem());
351 ARTS_ASSERT(ext_bulk_gas.ncols() == Np);
352
353 // Initialization
354 ext_bulk_gas = 0.;
355
356 // making gas property output containers and input dummies
357 const EnergyLevelMap rtp_nlte_dummy;
358 const Vector rtp_mag_dummy(3, 0);
359 const Vector ppath_los_dummy;
360 StokesVector nlte_dummy;
361 ArrayOfStokesVector partial_nlte_dummy;
362 ArrayOfPropagationMatrix partial_dummy;
363
364 PropagationMatrix propmat_clearsky_local;
365 for (Index ip = 0; ip < Np; ip++) {
367 propmat_clearsky_local,
368 nlte_dummy,
369 partial_dummy,
370 partial_nlte_dummy,
372 {},
373 f_grid,
374 rtp_mag_dummy,
375 ppath_los_dummy,
376 p_grid[ip],
377 t_profile[ip],
378 rtp_nlte_dummy,
379 vmr_profiles(joker, ip),
380 propmat_clearsky_agenda);
381 ext_bulk_gas(joker, ip) += propmat_clearsky_local.Kjj();
382 }
383}
384
386 MatrixView& sca_coeff_gas,
387 MatrixView& sca_coeff_gas_level,
388 MatrixView& pfct_gas,
389 const ConstVectorView& f_grid,
390 const VectorView& p,
391 const VectorView& t,
392 const MatrixView& vmr,
393 const Agenda& gas_scattering_agenda) {
394 const Index Np = p.nelem(); // Number of pressure levels
395 const Index Nl = pfct_gas.ncols(); // Number of legendre polynomials
396 const Index Nf = f_grid.nelem(); // Number of frequencies
397
398 PropagationMatrix K_sca_gas_temp;
399 TransmissionMatrix sca_mat_dummy;
400 Vector dummy;
401 Vector sca_fct_temp;
402 Matrix pmom_gas_level( Np, Nl, 0.);
403 Index N_polys;
404
405 // calculate gas scattering properties on level
406 for (Index ip = 0; ip < Np; ip++) {
408 K_sca_gas_temp,
409 sca_mat_dummy,
410 sca_fct_temp,
411 f_grid,
412 p[ip],
413 t[ip],
414 vmr(joker, ip),
415 dummy,
416 dummy,
417 1,
418 gas_scattering_agenda);
419
420 // gas scattering extinction
421 sca_coeff_gas_level(joker, ip) = K_sca_gas_temp.Kjj(0, 0);
422
423 // gas scattering (phase) function
424 N_polys = min(Nl, sca_fct_temp.nelem());
425 for (Index k = 0; k < N_polys; k++) {
426 pmom_gas_level( ip, k) = sca_fct_temp[k];
427 }
428 }
429
430 // layer averages
431 for (Index ip = 0; ip < Np - 1; ip++) {
432 for (Index f_index = 0; f_index < Nf; f_index++) {
433 // extinction
434 sca_coeff_gas(f_index, Np - 2 - ip) =
435 0.5 *
436 (sca_coeff_gas_level(f_index, ip) +
437 sca_coeff_gas_level(f_index, ip + 1));
438 }
439 // phase function
440 for (Index l_index = 0; l_index < Nl; l_index++) {
441 pfct_gas( Np - 2 - ip, l_index) =
442 0.5 * (pmom_gas_level( ip, l_index) +
443 pmom_gas_level( ip + 1, l_index));
444 }
445 }
446}
447
448void get_paroptprop(MatrixView ext_bulk_par,
449 MatrixView abs_bulk_par,
450 const ArrayOfArrayOfSingleScatteringData& scat_data,
451 ConstMatrixView pnd_profiles,
452 ConstVectorView t_profile,
454 const ArrayOfIndex& cloudbox_limits,
455 ConstVectorView f_grid) {
456 const Index Np_cloud = pnd_profiles.ncols();
457 DEBUG_ONLY(const Index Np = p_grid.nelem());
458 const Index nf = f_grid.nelem();
459
460 ARTS_ASSERT(ext_bulk_par.nrows() == nf);
461 ARTS_ASSERT(abs_bulk_par.nrows() == nf);
462 ARTS_ASSERT(ext_bulk_par.ncols() == Np);
463 ARTS_ASSERT(abs_bulk_par.ncols() == Np);
464
465 // Initialization
466 ext_bulk_par = 0.;
467 abs_bulk_par = 0.;
468
469 // preparing input data
470 Vector T_array = t_profile[Range(cloudbox_limits[0], Np_cloud)];
471 Matrix dir_array(1, 2, 0.); // just a dummy. only tot_random allowed, ie.
472 // optprop are independent of direction.
473
474 // making particle property output containers
475 ArrayOfArrayOfTensor5 ext_mat_Nse;
476 ArrayOfArrayOfTensor4 abs_vec_Nse;
477 ArrayOfArrayOfIndex ptypes_Nse;
478 Matrix t_ok;
479 ArrayOfTensor5 ext_mat_ssbulk;
480 ArrayOfTensor4 abs_vec_ssbulk;
481 ArrayOfIndex ptype_ssbulk;
482 Tensor5 ext_mat_bulk; //nf,nT,ndir,nst,nst
483 Tensor4 abs_vec_bulk;
484 Index ptype_bulk;
485
486 // calculate particle optical properties
487 opt_prop_NScatElems(ext_mat_Nse,
488 abs_vec_Nse,
489 ptypes_Nse,
490 t_ok,
491 scat_data,
492 1,
493 T_array,
494 dir_array,
495 -1);
496 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
497 abs_vec_ssbulk,
498 ptype_ssbulk,
499 ext_mat_Nse,
500 abs_vec_Nse,
501 ptypes_Nse,
502 pnd_profiles,
503 t_ok);
504 opt_prop_Bulk(ext_mat_bulk,
505 abs_vec_bulk,
506 ptype_bulk,
507 ext_mat_ssbulk,
508 abs_vec_ssbulk,
509 ptype_ssbulk);
510
511 Index f_this = 0;
512 bool pf = (abs_vec_bulk.nbooks() != 1);
513 for (Index ip = 0; ip < Np_cloud; ip++)
514 for (Index f_index = 0; f_index < nf; f_index++) {
515 if (pf) f_this = f_index;
516 ext_bulk_par(f_index, ip + cloudbox_limits[0]) =
517 ext_mat_bulk(f_this, ip, 0, 0, 0);
518 abs_bulk_par(f_index, ip + cloudbox_limits[0]) =
519 abs_vec_bulk(f_this, ip, 0, 0);
520 }
521}
522
524 MatrixView ssalb,
525 ConstMatrixView ext_bulk_gas,
526 ConstMatrixView ext_bulk_par,
527 ConstMatrixView abs_bulk_par,
528 ConstVectorView z_profile) {
529 const Index nf = ext_bulk_gas.nrows();
530 const Index Np = ext_bulk_gas.ncols();
531
532 ARTS_ASSERT(dtauc.nrows() == nf);
533 ARTS_ASSERT(ssalb.nrows() == nf);
534 ARTS_ASSERT(dtauc.ncols() == Np - 1);
535 ARTS_ASSERT(ssalb.ncols() == Np - 1);
536
537 // Initialization
538 dtauc = 0.;
539 ssalb = 0.;
540
541 for (Index ip = 0; ip < Np - 1; ip++)
542 // Do layer averaging and derive single scattering albedo & optical depth
543 for (Index f_index = 0; f_index < nf; f_index++) {
544 Numeric ext =
545 0.5 * (ext_bulk_gas(f_index, ip) + ext_bulk_par(f_index, ip) +
546 ext_bulk_gas(f_index, ip + 1) + ext_bulk_par(f_index, ip + 1));
547 if (ext != 0) {
548 Numeric abs =
549 0.5 *
550 (ext_bulk_gas(f_index, ip) + abs_bulk_par(f_index, ip) +
551 ext_bulk_gas(f_index, ip + 1) + abs_bulk_par(f_index, ip + 1));
552 ssalb(f_index, Np - 2 - ip) = (ext - abs) / ext;
553
554 ARTS_USER_ERROR_IF((ext - abs) / ext > 1,
555 "ssalb > 1 @ \n",
556 "pressure level = ", ip,"\n",
557 "frequency number = ", f_index,"\n");
558
559 }
560
561 dtauc(f_index, Np - 2 - ip) = ext * (z_profile[ip + 1] - z_profile[ip]);
562 }
563}
564
565void get_angs(Vector& pfct_angs,
566 const ArrayOfArrayOfSingleScatteringData& scat_data,
567 const Index& Npfct) {
568 const Index min_nang = 3;
569 Index nang = Npfct;
570
571 if (Npfct < 0) {
572 Index this_ss = 0, this_se = 0;
573 // determine nang and pfct_angs from scat_data with finest za_grid
574 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
575 for (Index i_se = scat_data[i_ss].nelem() - 1; i_se >= 0; i_se--)
576 // considering scat elems within one species mostly sorted from small to
577 // large sizes with large sizes corresponding to large za_grid. that is,
578 // starting searching from back should trigger if-clause (and variable
579 // update) less often.
580 if (nang < scat_data[i_ss][i_se].za_grid.nelem()) {
581 nang = scat_data[i_ss][i_se].za_grid.nelem();
582 this_ss = i_ss;
583 this_se = i_se;
584 }
585 pfct_angs = scat_data[this_ss][this_se].za_grid;
586 } else if (Npfct < min_nang) {
587 ostringstream os;
588 os << "Number of requested angular grid points (Npfct=" << Npfct
589 << ") is insufficient.\n"
590 << "At least " << min_nang << " points required.\n";
591 throw runtime_error(os.str());
592 } else {
593 nlinspace(pfct_angs, 0, 180, nang);
594 }
595}
596
597void get_parZ(Tensor3& pha_bulk_par,
598 const ArrayOfArrayOfSingleScatteringData& scat_data,
599 ConstMatrixView pnd_profiles,
600 ConstVectorView t_profile,
601 ConstVectorView pfct_angs,
602 const ArrayOfIndex& cloudbox_limits) {
603 const Index Np_cloud = pnd_profiles.ncols();
604 const Index nang = pfct_angs.nelem();
605
606 // Initialization
607 pha_bulk_par = 0.;
608
609 // preparing input data
610 Vector T_array = t_profile[Range(cloudbox_limits[0], Np_cloud)];
611 Matrix idir_array(1, 2, 0.); // we want pfct on sca ang grid, hence set
612 // pdir(*,0) to sca ang, all other to 0.
613 Matrix pdir_array(nang, 2, 0.);
614 pdir_array(joker, 0) = pfct_angs;
615
616 // making particle property output containers
617 ArrayOfArrayOfTensor6 pha_mat_Nse;
618 ArrayOfArrayOfIndex ptypes_Nse;
619 Matrix t_ok;
620 ArrayOfTensor6 pha_mat_ssbulk;
621 ArrayOfIndex ptype_ssbulk;
622 Tensor6 pha_mat_bulk; //nf,nT,npdir,nidir,nst,nst
623 Index ptype_bulk;
624
625 // calculate phase matrix
626 // FIXME: might be optimized by instead just executing pha_mat_1ScatElem where
627 // ext_bulk_par or pnd_field are non-zero.
628 pha_mat_NScatElems(pha_mat_Nse,
629 ptypes_Nse,
630 t_ok,
631 scat_data,
632 1,
633 T_array,
634 pdir_array,
635 idir_array,
636 -1);
637 pha_mat_ScatSpecBulk(pha_mat_ssbulk,
638 ptype_ssbulk,
639 pha_mat_Nse,
640 ptypes_Nse,
641 pnd_profiles,
642 t_ok);
643 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
644
645 pha_bulk_par(joker, Range(cloudbox_limits[0], Np_cloud), joker) =
646 pha_mat_bulk(joker, joker, joker, 0, 0, 0);
647}
648
649void get_pfct(Tensor3& pfct_bulk_par,
650 ConstTensor3View& pha_bulk_par,
651 ConstMatrixView ext_bulk_par,
652 ConstMatrixView abs_bulk_par,
653 const ArrayOfIndex& cloudbox_limits) {
654 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
655 const Index Np = pha_bulk_par.nrows();
656 const Index nf = pha_bulk_par.npages();
657 Index nang = pha_bulk_par.ncols();
658
659 ARTS_ASSERT(pfct_bulk_par.npages() == nf);
660 ARTS_ASSERT(pfct_bulk_par.nrows() == Np - 1);
661 ARTS_ASSERT(pfct_bulk_par.ncols() == nang);
662
663 // Initialization
664 pfct_bulk_par = 0.;
665
666 for (Index ip = cloudbox_limits[0]; ip < Np_cloud - 1; ip++)
667 for (Index f_index = 0; f_index < nf; f_index++) {
668 // Calculate layer averaged scattering (omitting factor 0.5 here as
669 // omitting it for layer averaged Z below, too).
670 Numeric sca =
671 (ext_bulk_par(f_index, ip) + ext_bulk_par(f_index, ip + 1)) -
672 (abs_bulk_par(f_index, ip) + abs_bulk_par(f_index, ip + 1));
673 if (sca != 0.) {
674 // Calculate layer averaged Z (omitting factor 0.5) and rescale from
675 // Z (Csca) to P (4Pi)
676 for (Index ia = 0; ia < nang; ia++)
677 pfct_bulk_par(f_index, Np - 2 - ip, ia) +=
678 pha_bulk_par(f_index, ip, ia) + pha_bulk_par(f_index, ip + 1, ia);
679 pfct_bulk_par(f_index, Np - 2 - ip, joker) *= 4 * PI / sca;
680 }
681 }
682}
683
685 ConstTensor3View pfct_bulk_par,
686 ConstVectorView pfct_angs,
687 const Index& Nlegendre) {
688 const Index nf = pfct_bulk_par.npages();
689 const Index nlyr = pfct_bulk_par.nrows();
690 const Index nang = pfct_bulk_par.ncols();
691
692 ARTS_ASSERT(nang == pfct_angs.nelem());
693
694 ARTS_ASSERT(pmom.npages() == nf);
695 ARTS_ASSERT(pmom.nrows() == nlyr);
696 ARTS_ASSERT(pmom.ncols() == Nlegendre);
697
698 Numeric pfct_threshold = 0.1;
699
700 // Initialization
701 pmom = 0.;
702
703 // we need the cosine of the pfct angles
704 Vector u(nang), adu(nang - 1), ad_angs(nang - 1);
705 Tensor3 px(nang - 1, Nlegendre, 2, 0.);
706 u[0] = cos(pfct_angs[0] * PI / 180.);
707 px(joker, 0, joker) = 1.;
708 for (Index ia = 1; ia < nang; ia++) {
709 u[ia] = cos(pfct_angs[ia] * PI / 180.);
710 adu[ia - 1] = abs(u[ia] - u[ia - 1]);
711 ad_angs[ia - 1] =
712 abs(pfct_angs[ia] * PI / 180. - pfct_angs[ia - 1] * PI / 180.);
713 px(ia - 1, 1, 0) = u[ia - 1];
714 px(ia - 1, 1, 1) = u[ia];
715 for (Index l = 2; l < Nlegendre; l++) {
716 Numeric dl = (double)l;
717 px(ia - 1, l, 0) = (2 * dl - 1) / dl * u[ia - 1] * px(ia - 1, l - 1, 0) -
718 (dl - 1) / dl * px(ia - 1, l - 2, 0);
719 px(ia - 1, l, 1) = (2 * dl - 1) / dl * u[ia] * px(ia - 1, l - 1, 1) -
720 (dl - 1) / dl * px(ia - 1, l - 2, 1);
721 }
722 }
723
724 for (Index il = 0; il < nlyr; il++)
725 if (pfct_bulk_par(joker, il, 0).sum() != 0.)
726 for (Index f_index = 0; f_index < nf; f_index++) {
727 if (pfct_bulk_par(f_index, il, 0) != 0) {
728 Vector pfct = pfct_bulk_par(f_index, il, joker);
729
730 // Check if phase function is properly normalized
731 // For highly peaked phasefunctions, integrating over the angle instead
732 // of over the cosine of the angle is numerically more exact though both
733 // ways are analytically equal. Furthermore in *scat_dataCalc* the
734 // integration is also done over the angle.
735 Numeric pint = 0.;
736 for (Index ia = 0; ia < nang - 1; ia++) {
737 pint += 0.5 * ad_angs[ia] *
738 (pfct[ia] * sin(pfct_angs[ia] * PI / 180.) +
739 pfct[ia + 1] * sin(pfct_angs[ia + 1] * PI / 180.));
740 }
741
742 if (abs(pint / 2. - 1.) > pfct_threshold) {
743 ostringstream os;
744 os << "Phase function normalization deviates from expected value by\n"
745 << 1e2 * pint / 2. - 1e2 << "(allowed: " << pfct_threshold * 1e2
746 << "%).\n"
747 << "Occurs at layer #" << il << " and frequency #" << f_index
748 << ".\n"
749 << "Something is wrong with your scattering data. Check!\n";
750 throw runtime_error(os.str());
751 }
752
753 // for the rest, rescale pfct to norm 2
754 pfct *= 2. / pint;
755
756 pmom(f_index, il, 0) = 1.;
757 for (Index ia = 0; ia < nang - 1; ia++) {
758 for (Index l = 1; l < Nlegendre; l++) {
759 pmom(f_index, il, l) +=
760 0.25 * ad_angs[ia] *
761 (px(ia, l, 0) * pfct[ia] * sin(pfct_angs[ia] * PI / 180.) +
762 px(ia, l, 1) * pfct[ia + 1] *
763 sin(pfct_angs[ia + 1] * PI / 180.));
764 }
765 }
766 }
767 }
768}
769
770void get_scat_bulk_layer(MatrixView& sca_bulk_layer,
771 const MatrixView& ext_bulk,
772 const MatrixView& abs_bulk) {
773 const Index nf = ext_bulk.nrows();
774 const Index Np = ext_bulk.ncols();
775
776 ARTS_ASSERT(sca_bulk_layer.nrows() == nf);
777 ARTS_ASSERT(sca_bulk_layer.ncols() == Np - 1);
778 ARTS_ASSERT(abs_bulk.nrows() == nf);
779 ARTS_ASSERT(abs_bulk.ncols() == Np);
780
781 // Initialization
782 sca_bulk_layer = 0.;
783
784 for (Index ip = 0; ip < Np - 1; ip++)
785 // Do layer averaging and derive single scattering albedo & optical depth
786 for (Index f_index = 0; f_index < nf; f_index++) {
787 Numeric sca =
788 0.5 * (ext_bulk(f_index, ip) - abs_bulk(f_index, ip) +
789 ext_bulk(f_index, ip + 1) - abs_bulk(f_index, ip + 1));
790
791 sca_bulk_layer(f_index, Np - 2 - ip) = sca;
792 }
793}
794
795// Use a thread_local variable to communicate the Verbosity to the
796// Disort error and warning functions. Ugly workaround, to avoid
797// passing a Verbosity argument throughout the whole cdisort code.
798// We want to avoid changes to the original code to keep it maintainable
799// in respect to upstream updates.
801
802#define MAX_WARNINGS 100
803
805void c_errmsg(const char* messag, int type) {
806 Verbosity verbosity = disort_verbosity;
807 static int warning_limit = FALSE, num_warnings = 0;
808
809 if (type == DS_ERROR) {
811 out0 << " ******* ERROR >>>>>> " << messag << "\n";
812 arts_exit(1);
813 }
814
815 if (warning_limit) return;
816
817 if (++num_warnings <= MAX_WARNINGS) {
819 out1 << " ******* WARNING >>>>>> " << messag << "\n";
820 } else {
822 out1 << " >>>>>> TOO MANY WARNING MESSAGES -- They will no longer be "
823 "printed <<<<<<<\n\n";
824 warning_limit = TRUE;
825 }
826
827 return;
828}
829
830#undef MAX_WARNINGS
831
833int c_write_bad_var(int quiet, const char* varnam) {
834 const int maxmsg = 50;
835 static int nummsg = 0;
836
837 nummsg++;
838 if (quiet != QUIET) {
839 Verbosity verbosity = disort_verbosity;
841 out1 << " **** Input variable " << varnam << " in error ****\n";
842 if (nummsg == maxmsg) {
843 c_errmsg("Too many input errors. Aborting...", DS_ERROR);
844 }
845 }
846
847 return TRUE;
848}
849
851int c_write_too_small_dim(int quiet, const char* dimnam, int minval) {
852 if (quiet != QUIET) {
853 Verbosity verbosity = disort_verbosity;
855 out1 << " **** Symbolic dimension " << dimnam
856 << " should be increased to at least " << minval << " ****\n";
857 }
858
859 return TRUE;
860}
861
863 Vector& z,
864 Vector& t,
865 Matrix& vmr,
866 Matrix& pnd,
867 ArrayOfIndex& cboxlims,
868 Index& ncboxremoved,
869 ConstVectorView p_grid,
870 ConstVectorView z_profile,
871 const Numeric& z_surface,
872 ConstVectorView t_profile,
873 ConstMatrixView vmr_profiles,
874 ConstMatrixView pnd_profiles,
875 const ArrayOfIndex& cloudbox_limits) {
876 // Surface at p_grid[0] and we just need to copy the original data
877 if (abs(z_surface - z_profile[0]) < 1e-3) {
878 p = p_grid;
879 z = z_profile;
880 t = t_profile;
881 vmr = vmr_profiles;
882 pnd = pnd_profiles;
883 cboxlims = cloudbox_limits;
884 ncboxremoved = 0;
885 }
886 // Surface above p_grid[0]
887 else {
888 // Some counters
889 Index np = p_grid.nelem(), ifirst = 0;
890 // Determine where to start with respect to z_profile
891 for (; z_surface >= z_profile[ifirst + 1]; ifirst++) {
892 }
893 np -= ifirst;
894 // Start by copying from ifirst to end
895 Range ind(ifirst, np);
896 p = p_grid[ind];
897 z = z_profile[ind];
898 t = t_profile[ind];
899 vmr = vmr_profiles(joker, ind);
900 // Insert surface altitude
901 z[0] = z_surface;
902 // Prepare interpolation
903 ArrayOfGridPos gp(1);
904 gridpos(gp[0], z_profile, z_surface);
905 Vector itw(2);
906 interpweights(itw, gp[0]);
907 // t and vmr
908 t[0] = interp(itw, t, gp[0]);
909 for (int i = 0; i < vmr.nrows(); i++) {
910 vmr(i, 0) = interp(itw, vmr(i, joker), gp[0]);
911 }
912 // p (we need a matrix version of iwt to use the function *itw2p*)
913 Matrix itw2(1, 2);
914 itw2(0, 0) = itw[0];
915 itw2(0, 1) = itw[1];
916 itw2p(p[0], p, gp, itw2);
917 // pnd_field and cloudbox limits need special treatment
918 cboxlims = cloudbox_limits;
919 if (ifirst < cloudbox_limits[0]) { // Surface below cloudbox
920 cboxlims[0] -= ifirst;
921 cboxlims[1] -= ifirst;
922 pnd = pnd_profiles;
923 ncboxremoved = 0;
924 } else { // Surface inside cloudbox
925 ncboxremoved = ifirst - cboxlims[0];
926 cboxlims[0] = 0;
927 cboxlims[1] = cloudbox_limits[1] - cloudbox_limits[0] - ncboxremoved;
928 ind = Range(ncboxremoved, cboxlims[1] + 1);
929 pnd = pnd_profiles(joker, ind);
930 gp[0].idx -= cloudbox_limits[0] + ncboxremoved;
931 for (int i = 0; i < pnd.nrows(); i++) {
932 pnd(i, 0) = interp(itw, pnd(i, joker), gp[0]);
933 }
934 }
935 }
936}
937
939 Tensor7& cloudbox_field,
940 ArrayOfMatrix& disort_aux,
941 ConstVectorView f_grid,
942 ConstVectorView p_grid,
943 ConstVectorView z_profile,
944 const Numeric& z_surface,
945 ConstVectorView t_profile,
946 ConstMatrixView vmr_profiles,
947 ConstMatrixView pnd_profiles,
948 const ArrayOfArrayOfSingleScatteringData& scat_data,
949 const ArrayOfStar& stars,
950 const Agenda& propmat_clearsky_agenda,
951 const Agenda& gas_scattering_agenda,
952 const ArrayOfIndex& cloudbox_limits,
953 const Numeric& surface_skin_t,
954 const Vector& surface_scalar_reflectivity,
955 ConstVectorView za_grid,
956 ConstVectorView aa_grid,
957 ConstVectorView star_rte_los,
958 const Index& gas_scattering_do,
959 const Index& stars_do,
960 const ArrayOfString& disort_aux_vars,
961 const Numeric& scale_factor,
962 const Index& nstreams,
963 const Index& Npfct,
964 const Index& only_tro,
965 const Index& quiet,
966 const Index& emission,
967 const Index& intensity_correction,
968 const Verbosity& verbosity) {
969 // Create an atmosphere starting at z_surface
970 Vector p, z, t;
971 Matrix vmr, pnd;
972 ArrayOfIndex cboxlims;
973 Index ncboxremoved;
974 //
976 z,
977 t,
978 vmr,
979 pnd,
980 cboxlims,
981 ncboxremoved,
982 p_grid,
983 z_profile,
984 z_surface,
985 t_profile,
986 vmr_profiles,
987 pnd_profiles,
988 cloudbox_limits);
989
990 disort_state ds;
991 disort_output out;
992
993 if (quiet == 0)
994 disort_verbosity = verbosity;
995 else
996 disort_verbosity = Verbosity(0, 0, 0);
997
998 const Index nf = f_grid.nelem();
999
1000 // solar dependent properties if no sun is present
1001 // Number of azimuth angles
1002 Index nphi = 1;
1003 //local zenith angle of sun
1004 Numeric umu0 = 0.;
1005 //local azimuth angle of sun
1006 Numeric phi0 = 0.;
1007 //Intensity of incident sun beam
1008 Numeric fbeam = 0.;
1009
1010 if (stars_do) {
1011 nphi = aa_grid.nelem();
1012 umu0 = Conversion::cosd(star_rte_los[0]);
1013 phi0 = star_rte_los[1];
1014 if (phi0 < 0) {
1015 phi0 = phi0 + 360.;
1016 }
1017 }
1018
1019 ds.accur = 0.005;
1020 ds.flag.prnt[0] = FALSE;
1021 ds.flag.prnt[1] = FALSE;
1022 ds.flag.prnt[2] = FALSE;
1023 ds.flag.prnt[3] = FALSE;
1024 ds.flag.prnt[4] = TRUE;
1025
1026 ds.flag.usrtau = FALSE;
1027 ds.flag.usrang = TRUE;
1028 ds.flag.spher = FALSE;
1029 ds.flag.general_source = FALSE;
1030 ds.flag.output_uum = FALSE;
1031
1032 ds.nlyr = static_cast<int>(p.nelem() - 1);
1033
1034 ds.flag.brdf_type = BRDF_NONE;
1035
1036 ds.flag.ibcnd = GENERAL_BC;
1037 ds.flag.usrang = TRUE;
1038
1039 if (emission) {
1040 ds.flag.planck = TRUE;
1041 } else {
1042 ds.flag.planck = FALSE;
1043 }
1044 ds.flag.onlyfl = FALSE;
1045 ds.flag.lamber = TRUE;
1046 ds.flag.quiet = FALSE;
1047 if (intensity_correction) {
1048 ds.flag.intensity_correction = TRUE;
1049 ds.flag.old_intensity_correction = FALSE;
1050 } else {
1051 ds.flag.intensity_correction = FALSE;
1052 ds.flag.old_intensity_correction = FALSE;
1053 }
1054
1055 ds.nstr = static_cast<int>(nstreams);
1056 ds.nphase = ds.nstr;
1057 ds.nmom = ds.nstr;
1058 //ds.ntau = ds.nlyr + 1; // With ds.flag.usrtau = FALSE; set by cdisort
1059 ds.numu = static_cast<int>(za_grid.nelem());
1060 ds.nphi = static_cast<int>(nphi);
1061 Index Nlegendre = nstreams + 1;
1062
1063 /* Allocate memory */
1064 c_disort_state_alloc(&ds);
1065 c_disort_out_alloc(&ds, &out);
1066
1067 // Looking direction of solar beam
1068 ds.bc.umu0 = umu0;
1069 ds.bc.phi0 = phi0;
1070
1071 // Intensity of bottom-boundary isotropic illumination
1072 ds.bc.fluor = 0.;
1073
1074 // fill up azimuth angle and temperature array
1075 for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i];
1076
1077 if (ds.flag.planck==TRUE){
1078 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1079 }
1080
1081 // Transform to mu, starting with negative values
1082 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);
1083
1084 //gas absorption
1085 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1086 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1087
1088 // Get particle bulk properties
1089 Index nang;
1090 Vector pfct_angs;
1091 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1092 Index nf_ssd = scat_data[0][0].f_grid.nelem();
1093 Tensor3 pha_bulk_par;
1094
1095 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1096 nang = Npfct;
1097 nlinspace(pfct_angs, 0, 180, nang);
1098
1099 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1100
1101 ext_bulk_par = 0.0;
1102 abs_bulk_par = 0.0;
1103 pha_bulk_par = 0.0;
1104
1105 Index iflat = 0;
1106
1107 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1108 const Index nse = scat_data[iss].nelem();
1109 ext_abs_pfun_from_tro(ext_bulk_par,
1110 abs_bulk_par,
1111 pha_bulk_par,
1112 scat_data[iss],
1113 iss,
1114 pnd(Range(iflat, nse), joker),
1115 cboxlims,
1116 t,
1117 pfct_angs);
1118 iflat += nse;
1119 }
1120 } else {
1121 get_angs(pfct_angs, scat_data, Npfct);
1122 nang = pfct_angs.nelem();
1123
1124 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1125
1127 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1128 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1129 }
1130
1131 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1132 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1133
1134 // Legendre polynomials of phase function
1135 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1136 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1137
1138 if (gas_scattering_do) {
1139 // gas scattering
1140
1141 // layer averaged particle scattering coefficient
1142 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1143 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1144
1145 // call gas_scattering_properties
1146 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1147 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1148 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1149
1151 sca_coeff_gas_layer,
1152 sca_coeff_gas_level,
1153 pmom_gas,
1154 f_grid,
1155 p,
1156 t,
1157 vmr,
1158 gas_scattering_agenda);
1159
1160 // call add_norm_phase_functions
1162 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1163
1164 // add gas_scat_ext to ext_bulk_par
1165 ext_bulk_par += sca_coeff_gas_level;
1166 }
1167
1168 // Optical depth of layers
1169 Matrix dtauc(nf, ds.nlyr);
1170 // Single scattering albedo of layers
1171 Matrix ssalb(nf, ds.nlyr);
1172 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1173
1174 //upper boundary conditions:
1175 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1176 // emission. Both are applied additively.
1177 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1178 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1179 // are additive we should use either the one or the other.
1180 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1181 // Generally, temis!=1 should be avoided since that technically implies a
1182 // reflective upper boundary (though it seems that this is not exploited in
1183 // DISORT1.2, which we so far use).
1184
1185 // Cosmic background
1186 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1187 ds.bc.fisot = 0;
1188
1189 // Top of the atmosphere temperature and emissivity
1190 ds.bc.ttemp = COSMIC_BG_TEMP;
1191 ds.bc.btemp = surface_skin_t;
1192 ds.bc.temis = 1.;
1193
1194 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1195 sprintf(ds.header, "ARTS Calc f_index = %" PRId64, f_index);
1196
1197 std::memcpy(ds.dtauc,
1198 dtauc(f_index, joker).get_c_array(),
1199 sizeof(Numeric) * ds.nlyr);
1200 std::memcpy(ds.ssalb,
1201 ssalb(f_index, joker).get_c_array(),
1202 sizeof(Numeric) * ds.nlyr);
1203
1204 // Wavenumber in [1/cm]
1205 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
1206 ds.wvnmhi += ds.wvnmhi * 1e-7;
1207 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1208
1209 // set
1210 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1211
1212 // Set irradiance of incident solar beam at top boundary
1213 if (stars_do) {
1214 fbeam = stars[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1215 (100 * SPEED_OF_LIGHT)*scale_factor;
1216 }
1217 ds.bc.fbeam = fbeam;
1218
1219 std::memcpy(ds.pmom,
1220 pmom(f_index, joker, joker).get_c_array(),
1221 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1222
1223 c_disort(&ds, &out);
1224
1225 for (Index i = 0; i < ds.nphi; i++) {
1226 for (Index j = 0; j < ds.numu; j++) {
1227 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1228 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, i, 0) =
1229 out.uu[j + ((ds.nlyr - k - cboxlims[0]) + i * (ds.nlyr + 1)) *
1230 ds.numu] /
1231 (ds.wvnmhi - ds.wvnmlo) / (100 * SPEED_OF_LIGHT);
1232 }
1233 // To avoid potential numerical problems at interpolation of the field,
1234 // we copy the surface field to underground altitudes
1235 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1236 cloudbox_field(f_index, k, 0, 0, j, i, 0) =
1237 cloudbox_field(f_index, k + 1, 0, 0, j, i, 0);
1238 }
1239 }
1240 }
1241 }
1242
1243 // Allocate aux data
1244 disort_aux.resize(disort_aux_vars.nelem());
1245 // Allocate and set (if possible here) iy_aux
1246 Index cnt=-1;
1247 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1248
1249
1250 if (disort_aux_vars[i] == "Layer optical thickness"){
1251 cnt += 1;
1252 Matrix deltatau(nf, ds.nlyr, 0);
1253 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1254 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1255 deltatau(f_index, k - 1 + ncboxremoved) =
1256 dtauc(f_index, ds.nlyr - k + ncboxremoved);
1257 }
1258 }
1259 disort_aux[cnt] = deltatau;
1260 }
1261 else if (disort_aux_vars[i] == "Single scattering albedo"){
1262 cnt+=1;
1263 Matrix snglsctalbedo(nf, ds.nlyr,0);
1264 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1265 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1266 snglsctalbedo(f_index, k - 1 + ncboxremoved) =
1267 ssalb(f_index, ds.nlyr - k + ncboxremoved);
1268 }
1269 }
1270 disort_aux[cnt]=snglsctalbedo;
1271 }
1272 else if (disort_aux_vars[i] == "Direct beam") {
1273 cnt += 1;
1274 Matrix directbeam(nf, ds.nlyr + 1, 0);
1275 if (stars_do) {
1276 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1277 directbeam(f_index, cboxlims[1] - cboxlims[0] + ncboxremoved) =
1278 stars[0].spectrum(f_index, 0)/PI;
1279
1280 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1281 directbeam(f_index, k - 1 + ncboxremoved) =
1282 directbeam(f_index, k + ncboxremoved) *
1283 exp(-dtauc(f_index, ds.nlyr - k + ncboxremoved)/umu0);
1284 }
1285 }
1286 }
1287 disort_aux[cnt]=directbeam;
1288 } else {
1290 "The only allowed strings in *disort_aux_vars* are:\n"
1291 " \"Layer optical thickness\"\n"
1292 " \"Single scattering albedo\"\n"
1293 " \"Direct beam\"\n"
1294 "but you have selected: \"", disort_aux_vars[i], "\"\n");
1295 }
1296 }
1297
1298 /* Free allocated memory */
1299 c_disort_out_free(&ds, &out);
1300 c_disort_state_free(&ds);
1301}
1302
1304 Tensor5& spectral_irradiance_field,
1305 ArrayOfMatrix& disort_aux,
1306 ConstVectorView f_grid,
1307 ConstVectorView p_grid,
1308 ConstVectorView z_profile,
1309 const Numeric& z_surface,
1310 ConstVectorView t_profile,
1311 ConstMatrixView vmr_profiles,
1312 ConstMatrixView pnd_profiles,
1313 const ArrayOfArrayOfSingleScatteringData& scat_data,
1314 const ArrayOfStar& stars,
1315 const Agenda& propmat_clearsky_agenda,
1316 const Agenda& gas_scattering_agenda,
1317 const ArrayOfIndex& cloudbox_limits,
1318 const Numeric& surface_skin_t,
1319 const Vector& surface_scalar_reflectivity,
1320 ConstVectorView star_rte_los,
1321 const Index& gas_scattering_do,
1322 const Index& stars_do,
1323 const ArrayOfString& disort_aux_vars,
1324 const Numeric& scale_factor,
1325 const Index& nstreams,
1326 const Index& Npfct,
1327 const Index& only_tro,
1328 const Index& quiet,
1329 const Index& emission,
1330 const Index& intensity_correction,
1331 const Verbosity& verbosity) {
1332 // Create an atmosphere starting at z_surface
1333 Vector p, z, t;
1334 Matrix vmr, pnd;
1335 ArrayOfIndex cboxlims;
1336 Index ncboxremoved;
1337 //
1338 reduced_1datm(p,
1339 z,
1340 t,
1341 vmr,
1342 pnd,
1343 cboxlims,
1344 ncboxremoved,
1345 p_grid,
1346 z_profile,
1347 z_surface,
1348 t_profile,
1349 vmr_profiles,
1350 pnd_profiles,
1351 cloudbox_limits);
1352
1353 disort_state ds;
1354 disort_output out;
1355
1356 if (quiet == 0)
1357 disort_verbosity = verbosity;
1358 else
1359 disort_verbosity = Verbosity(0, 0, 0);
1360
1361 const Index nf = f_grid.nelem();
1362
1363 // solar dependent properties if no sun is present
1364 //local zenith angle of sun
1365 Numeric umu0 = 0.;
1366 //local azimuth angle of sun
1367 Numeric phi0 = 0.;
1368 //Intensity of incident sun beam
1369 Numeric fbeam = 0.;
1370
1371 if (stars_do) {
1372 umu0 = Conversion::cosd(star_rte_los[0]);
1373 phi0 = star_rte_los[1];
1374 if (phi0 < 0) {
1375 phi0 = phi0 + 360.;
1376 }
1377 }
1378
1379 ds.accur = 0.005;
1380 ds.flag.prnt[0] = FALSE;
1381 ds.flag.prnt[1] = FALSE;
1382 ds.flag.prnt[2] = FALSE;
1383 ds.flag.prnt[3] = FALSE;
1384 ds.flag.prnt[4] = TRUE;
1385
1386 ds.flag.usrtau = FALSE;
1387 ds.flag.usrang = TRUE;
1388 ds.flag.spher = FALSE;
1389 ds.flag.general_source = FALSE;
1390 ds.flag.output_uum = FALSE;
1391
1392 ds.nlyr = static_cast<int>(p.nelem() - 1);
1393
1394 ds.flag.brdf_type = BRDF_NONE;
1395
1396 ds.flag.ibcnd = GENERAL_BC;
1397 ds.flag.usrang = FALSE;
1398
1399 if (emission) {
1400 ds.flag.planck = TRUE;
1401 } else {
1402 ds.flag.planck = FALSE;
1403 }
1404 ds.flag.onlyfl = TRUE;
1405 ds.flag.lamber = TRUE;
1406 ds.flag.quiet = FALSE;
1407 if (intensity_correction) {
1408 ds.flag.intensity_correction = TRUE;
1409 ds.flag.old_intensity_correction = FALSE;
1410 } else {
1411 ds.flag.intensity_correction = FALSE;
1412 ds.flag.old_intensity_correction = FALSE;
1413 }
1414
1415 ds.nstr = static_cast<int>(nstreams);
1416 ds.nphase = ds.nstr;
1417 ds.nmom = ds.nstr;
1418
1419 // set grid size of angular variable to dummy value
1420 ds.numu = 1;
1421 ds.nphi = 1;
1422
1423 //Set number of legendre terms
1424 Index Nlegendre = nstreams + 1;
1425
1426 /* Allocate memory */
1427 c_disort_state_alloc(&ds);
1428 c_disort_out_alloc(&ds, &out);
1429
1430 // Looking direction of solar beam
1431 ds.bc.umu0 = umu0;
1432 ds.bc.phi0 = phi0;
1433
1434 // Intensity of bottom-boundary isotropic illumination
1435 ds.bc.fluor = 0.;
1436
1437 // fill up temperature array
1438 if (ds.flag.planck==TRUE){
1439 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1440 }
1441
1442 //gas absorption
1443 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1444 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1445
1446 // Get particle bulk properties
1447 Index nang;
1448 Vector pfct_angs;
1449 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1450 Index nf_ssd = scat_data[0][0].f_grid.nelem();
1451 Tensor3 pha_bulk_par;
1452
1453 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1454 nang = Npfct;
1455 nlinspace(pfct_angs, 0, 180, nang);
1456
1457 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1458
1459 ext_bulk_par = 0.0;
1460 abs_bulk_par = 0.0;
1461 pha_bulk_par = 0.0;
1462
1463 Index iflat = 0;
1464
1465 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1466 const Index nse = scat_data[iss].nelem();
1467 ext_abs_pfun_from_tro(ext_bulk_par,
1468 abs_bulk_par,
1469 pha_bulk_par,
1470 scat_data[iss],
1471 iss,
1472 pnd(Range(iflat, nse), joker),
1473 cboxlims,
1474 t,
1475 pfct_angs);
1476 iflat += nse;
1477 }
1478 } else {
1479 get_angs(pfct_angs, scat_data, Npfct);
1480 nang = pfct_angs.nelem();
1481
1482 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1483
1485 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1486 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1487 }
1488
1489 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1490 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1491
1492 // Legendre polynomials of phase function
1493 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1494 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1495
1496 if (gas_scattering_do) {
1497 // gas scattering
1498
1499 // layer averaged particle scattering coefficient
1500 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1501 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1502
1503 // call gas_scattering_properties
1504 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1505 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1506 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1507
1509 sca_coeff_gas_layer,
1510 sca_coeff_gas_level,
1511 pmom_gas,
1512 f_grid,
1513 p,
1514 t,
1515 vmr,
1516 gas_scattering_agenda);
1517
1518 // call add_norm_phase_functions
1520 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1521
1522 // add gas_scat_ext to ext_bulk_par
1523 ext_bulk_par += sca_coeff_gas_level;
1524 }
1525
1526 // Optical depth of layers
1527 Matrix dtauc(nf, ds.nlyr);
1528 // Single scattering albedo of layers
1529 Matrix ssalb(nf, ds.nlyr);
1530 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1531
1532 //upper boundary conditions:
1533 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1534 // emission. Both are applied additively.
1535 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1536 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1537 // are additive we should use either the one or the other.
1538 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1539 // Generally, temis!=1 should be avoided since that technically implies a
1540 // reflective upper boundary (though it seems that this is not exploited in
1541 // DISORT1.2, which we so far use).
1542
1543 // Cosmic background
1544 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1545 ds.bc.fisot = 0;
1546
1547 // Top of the atmosphere temperature and emissivity
1548 ds.bc.ttemp = COSMIC_BG_TEMP;
1549 ds.bc.btemp = surface_skin_t;
1550 ds.bc.temis = 1.;
1551
1552 Matrix spectral_direct_irradiance_field;
1553 Matrix dFdtau(nf, ds.nlyr+1,0);
1554 Matrix deltatau(nf, ds.nlyr,0);
1555 Matrix snglsctalbedo(nf, ds.nlyr,0);
1556
1557 if (stars_do){
1558 //Resize direct field
1559 spectral_direct_irradiance_field.resize(nf, ds.nlyr+1);
1560 spectral_direct_irradiance_field = 0;
1561 }
1562
1563 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1564 sprintf(ds.header, "ARTS Calc f_index = %" PRId64, f_index);
1565
1566 std::memcpy(ds.dtauc,
1567 dtauc(f_index, joker).get_c_array(),
1568 sizeof(Numeric) * ds.nlyr);
1569 std::memcpy(ds.ssalb,
1570 ssalb(f_index, joker).get_c_array(),
1571 sizeof(Numeric) * ds.nlyr);
1572
1573 // Wavenumber in [1/cm]
1574 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
1575 ds.wvnmhi += ds.wvnmhi * 1e-7;
1576 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1577
1578 // set
1579 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1580
1581 // Set irradiance of incident solar beam at top boundary
1582 if (stars_do) {
1583 fbeam = stars[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1584 (100 * SPEED_OF_LIGHT)*scale_factor;
1585 }
1586 ds.bc.fbeam = fbeam;
1587
1588 std::memcpy(ds.pmom,
1589 pmom(f_index, joker, joker).get_c_array(),
1590 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1591
1592 c_disort(&ds, &out);
1593
1594 //factor for converting it into spectral radiance units
1595 const Numeric conv_fac=(ds.wvnmhi - ds.wvnmlo) * (100 * SPEED_OF_LIGHT);
1596
1597 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1598 if (stars_do){
1599 // downward direct flux
1600 spectral_direct_irradiance_field(f_index, k + ncboxremoved) =
1601 -out.rad[ds.nlyr - k - cboxlims[0]].rfldir/conv_fac;
1602
1603 // downward total flux
1604 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1605 -(out.rad[ds.nlyr - k - cboxlims[0]].rfldir +
1606 out.rad[ds.nlyr - k - cboxlims[0]].rfldn)/conv_fac;
1607
1608 } else {
1609 // downward total flux
1610 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1611 -out.rad[ds.nlyr - k - cboxlims[0]].rfldn/conv_fac;
1612 }
1613
1614 // upward flux
1615 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 1) =
1616 out.rad[ds.nlyr - k - cboxlims[0]].flup/conv_fac;
1617
1618 // flux divergence in tau space
1619 dFdtau(f_index, k + ncboxremoved) =
1620 -out.rad[ds.nlyr - k - cboxlims[0]].dfdt;
1621
1622 // k is running over the number of levels but deltatau, ssalb is defined for layers,
1623 // therefore we need to exlude k==0 and remove one from the index.
1624 if (k>0){
1625 deltatau(f_index, k - 1 + ncboxremoved) = ds.dtauc[ds.nlyr - k - 1 - cboxlims[0]];
1626 snglsctalbedo(f_index, k - 1 + ncboxremoved) = ds.ssalb[ds.nlyr - k - 1 - cboxlims[0]];
1627 }
1628 }
1629
1630 // To avoid potential numerical problems at interpolation of the field,
1631 // we copy the surface field to underground altitudes
1632 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1633 spectral_irradiance_field(f_index, k, 0, 0, joker) =
1634 spectral_irradiance_field(f_index, k + 1, 0, 0, joker);
1635
1636 if (stars_do) {
1637 spectral_direct_irradiance_field(f_index, k) =
1638 spectral_direct_irradiance_field(f_index, k + 1);
1639 }
1640 dFdtau(f_index, k) = dFdtau(f_index, k + 1);
1641 }
1642 }
1643
1644 // Allocate aux data
1645 disort_aux.resize(disort_aux_vars.nelem());
1646 // Allocate and set (if possible here) iy_aux
1647 Index cnt=-1;
1648 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1649
1650
1651 if (disort_aux_vars[i] == "Layer optical thickness"){
1652 cnt+=1;
1653 disort_aux[cnt]=deltatau;
1654 }
1655 else if (disort_aux_vars[i] == "Single scattering albedo"){
1656 cnt+=1;
1657 disort_aux[cnt]=snglsctalbedo;
1658 }
1659 else if (disort_aux_vars[i] == "Direct downward spectral irradiance") {
1660 cnt += 1;
1661 disort_aux[cnt] = spectral_direct_irradiance_field;
1662 }
1663 else if (disort_aux_vars[i] == "dFdtau") {
1664 cnt += 1;
1665 disort_aux[cnt] = dFdtau;
1666 }
1667 else {
1669 "The only allowed strings in *disort_aux_vars* are:\n"
1670 " \"Layer optical thickness\"\n"
1671 " \"Single scattering albedo\"\n"
1672 " \"Direct downward spectral irradiance\"\n"
1673 " \"dFdtau\"\n"
1674 "but you have selected: \"", disort_aux_vars[i], "\"\n");
1675 }
1676 }
1677
1678
1679
1680 /* Free allocated memory */
1681 c_disort_out_free(&ds, &out);
1682 c_disort_state_free(&ds);
1683}
1684
1686 //Output
1687 VectorView albedo,
1688 Numeric& btemp,
1689 //Input
1690 const Agenda& surface_rtprop_agenda,
1691 ConstVectorView f_grid,
1692 ConstVectorView scat_za_grid,
1693 const Numeric& surf_alt,
1694 const Verbosity& verbosity) {
1695 // Here, we derive an average surface albedo of the setup as given by ARTS'
1696 // surface_rtprop_agenda to use with Disorts's proprietary Lambertian surface.
1697 // In this way, ARTS-Disort can approximately mimick all surface reflection
1698 // types that ARTS itself can handle.
1699 // Surface temperature as derived from surface_rtprop_agenda is also returned.
1700 //
1701 // We derive the reflection matrices over all incident and reflected polar
1702 // angle directions and derive their integrated value (or weighted average).
1703 // The surface_rtprop_agenda handles one reflected direction at a time
1704 // (rtp_los) and for the reflected directions we loop over all (upwelling)
1705 // angles as given by scat_za_grid. The derived reflection matrices already
1706 // represent the reflectivity (or power reflection coefficient) for the given
1707 // reflected direction including proper angle weighting. For integrating/
1708 // averaging over the reflected directions, we use the same approach, i.e.
1709 // weight each angle by its associated range as given by the half distances to
1710 // the neighboring grid angles (using ARTS' scat_za_grid means 0/180deg are
1711 // grid points, 90deg shouldn't be among them (resulting from even number
1712 // requirement for Disort and the (implicitly assumed?) requirement of a
1713 // horizon-symmetric za_grid)).
1714 //
1715 // We do all frequencies here at once (assuming this is the faster variant as
1716 // the agenda anyway (can) provide output for full f_grid at once and as we
1717 // have to apply the same inter/extrapolation to all the frequencies).
1718
1720
1721 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1722
1723 const Index nf = f_grid.nelem();
1724 Index frza = 0;
1725 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
1726 if (frza == scat_za_grid.nelem()) {
1727 ostringstream os;
1728 os << "No upwelling direction found in scat_za_grid.\n";
1729 throw runtime_error(os.str());
1730 }
1731 const Index nrza = scat_za_grid.nelem() - frza;
1732 //cout << nrza << " upwelling directions found, starting from element #"
1733 // << frza << " of scat_za_grid.\n";
1734 Matrix dir_refl_coeff(nrza, nf, 0.);
1735
1736 // Local input of surface_rtprop_agenda.
1737 Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
1738
1739 // first derive the (reflected-)direction dependent power reflection
1740 // coefficient
1741 for (Index rza = 0; rza < nrza; rza++) {
1742 // Local output of surface_rtprop_agenda.
1743 Numeric surface_skin_t;
1744 Matrix surface_los;
1745 Tensor4 surface_rmatrix;
1746 Matrix surface_emission;
1747
1748 Vector rtp_los(1, scat_za_grid[rza + frza]);
1749 out2 << "Doing reflected dir #" << rza << " at " << rtp_los[0] << " degs\n";
1750
1752 surface_skin_t,
1753 surface_emission,
1754 surface_los,
1755 surface_rmatrix,
1756 f_grid,
1757 rtp_pos,
1758 rtp_los,
1759 surface_rtprop_agenda);
1760 if (surface_los.nrows() > 1) {
1761 ostringstream os;
1762 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1763 << "return zero or one direction in *surface_los*.";
1764 throw runtime_error(os.str());
1765 }
1766
1767 if (rza == 0)
1768 btemp = surface_skin_t;
1769 else if (surface_skin_t != btemp) {
1770 ostringstream os;
1771 os << "Something went wrong.\n"
1772 << " *surface_rtprop_agenda* returned different surface_skin_t\n"
1773 << " for different LOS.\n";
1774 throw runtime_error(os.str());
1775 }
1776 // Nothing to do if no surface_los (which means blackbody surface)
1777 // as dir_refl_coeff already set to 0
1778 if (!surface_los.empty()) {
1779 for (Index f_index = 0; f_index < nf; f_index++)
1780 dir_refl_coeff(rza, f_index) =
1781 surface_rmatrix(joker, f_index, 0, 0).sum();
1782 }
1783 out2 << " directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1784 << "\n";
1785 }
1786
1787 if (btemp < 0. || btemp > 1000.) {
1788 ostringstream os;
1789 os << "Surface temperature has been derived as " << btemp << " K,\n"
1790 << "which is not considered a meaningful value.\n";
1791 throw runtime_error(os.str());
1792 }
1793
1794 // now integrate/average the (reflected-)direction dependent power reflection
1795 // coefficients
1796 //
1797 // starting with deriving the angles defining the angle ranges
1798 Vector surf_int_grid(nrza + 1);
1799 // the first angle grid point should be around (but above) 90deg and should
1800 // cover the angle range between the 90deg and half-way point towards the next
1801 // angle grid point. we probably also want to check, that we don't
1802 // 'extrapolate' too much.
1803 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1804 ostringstream os;
1805 os << "Looks like scat_za_grid contains the 90deg direction,\n"
1806 << "which it shouldn't for running Disort.\n";
1807 throw runtime_error(os.str());
1808 }
1809 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1810 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1811 const Numeric ok_extrapol = 0.5;
1812 if ((za_extrapol - ok_extrapol) > 1e-6) {
1813 ostringstream os;
1814 os << "Extrapolation range from shallowest scat_za_grid point\n"
1815 << "to horizon is too big.\n"
1816 << " Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1817 << ", which is " << za_extrapol - 0.5 << " too big.\n";
1818 throw runtime_error(os.str());
1819 }
1821 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1822 ostringstream os;
1823 os << "Looks like last point in scat_za_grid is not 180deg.\n";
1824 throw runtime_error(os.str());
1825 }
1826
1827 surf_int_grid[0] = 90.;
1828 surf_int_grid[nrza] = 180.;
1829 for (Index rza = 1; rza < nrza; rza++)
1830 surf_int_grid[rza] =
1831 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1832 surf_int_grid *= DEG2RAD;
1833
1834 // now calculating the actual weights and apply them
1835 for (Index rza = 0; rza < nrza; rza++) {
1836 //Numeric coslow = cos(2.*surf_int_grid[rza]);
1837 //Numeric coshigh = cos(2.*surf_int_grid[rza+1]);
1838 //Numeric w = 0.5*(coshigh-coslow);
1839 Numeric w =
1840 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1841 //cout << "at reflLOS[" << rza << "]=" << scat_za_grid[frza+rza] << ":\n";
1842 //cout << " angle weight derives as w = 0.5*(" << coshigh << "-"
1843 // << coslow << ") = " << w << "\n";
1844 //cout << " weighting directional reflection coefficient from "
1845 // << dir_refl_coeff(rza,joker);
1846 dir_refl_coeff(rza, joker) *= w;
1847 //cout << " to " << dir_refl_coeff(rza,joker) << "\n";
1848 }
1849
1850 // eventually sum up the weighted directional power reflection coefficients
1851 for (Index f_index = 0; f_index < nf; f_index++) {
1852 albedo[f_index] = dir_refl_coeff(joker, f_index).sum();
1853 out2 << "at f=" << f_grid[f_index] * 1e-9
1854 << " GHz, ending up with albedo=" << albedo[f_index] << "\n";
1855 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1856 ostringstream os;
1857 os << "Something went wrong: Albedo must be inside [0,1],\n"
1858 << " but is not at freq #" << f_index << " , where it is "
1859 << albedo[f_index] << ".\n";
1860 throw runtime_error(os.str());
1861 }
1862 }
1863}
1864
1866 //Output
1867 VectorView albedo,
1868 Numeric& btemp,
1869 //Input
1870 const Agenda& surface_rtprop_agenda,
1871 ConstVectorView f_grid,
1872 const Numeric& surf_alt,
1873 const Numeric& inc_angle) {
1874
1875 Vector rtp_pos(1, surf_alt);
1876 Vector rtp_los(1, 180-inc_angle);
1877
1878 Numeric surface_skin_t;
1879 Matrix surface_los;
1880 Tensor4 surface_rmatrix;
1881 Matrix surface_emission;
1882
1884 surface_skin_t,
1885 surface_emission,
1886 surface_los,
1887 surface_rmatrix,
1888 f_grid,
1889 rtp_pos,
1890 rtp_los,
1891 surface_rtprop_agenda);
1892
1893 if (surface_los.nrows() > 1) {
1894 ostringstream os;
1895 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1896 << "return zero or one direction in *surface_los*.";
1897 throw runtime_error(os.str());
1898 }
1899 if (surface_los.empty()) {
1900 albedo = 0; // Blackbody assumed if no surface_los
1901 } else {
1902 albedo[joker] = surface_rmatrix(0,joker,0,0);
1903 }
1904
1905 btemp = surface_skin_t;
1906 if (btemp < 0. || btemp > 1000.) {
1907 ostringstream os;
1908 os << "Surface temperature has been derived as " << btemp << " K,\n"
1909 << "which is not considered a meaningful value.\n";
1910 throw runtime_error(os.str());
1911 }
1912}
Declarations for agendas.
This file contains the definition of Array.
void arts_exit(int status)
This is the exit function of ARTS.
Definition: arts.cc:42
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:24844
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:25655
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:25839
void chk_not_empty(const String &x_name, const Agenda &x)
chk_not_empty
Definition: check_input.cc:627
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1043
bool empty() const noexcept
Definition: matpackI.h:1062
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Tensor3.
Definition: matpackIII.h:130
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index nbooks() const noexcept
Definition: matpackIV.h:139
A constant view of a Vector.
Definition: matpackI.h:512
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:159
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3View class.
Definition: matpackIII.h:246
The Tensor3 class.
Definition: matpackIII.h:346
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:429
The Tensor5 class.
Definition: matpackV.h:516
The Tensor6 class.
Definition: matpackVI.h:1099
The Tensor7 class.
Definition: matpackVII.h:2399
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5488
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
Workspace class.
Definition: workspace_ng.h:53
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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:1685
constexpr Numeric PLANCK_CONST
Definition: disort.cc:58
void c_errmsg(const char *messag, int type)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:805
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:385
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:649
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:862
constexpr Numeric SPEED_OF_LIGHT
Definition: disort.cc:59
#define MAX_WARNINGS
Definition: disort.cc:802
constexpr Numeric DEG2RAD
Definition: disort.cc:57
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:597
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:296
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:1865
void get_scat_bulk_layer(MatrixView &sca_bulk_layer, const MatrixView &ext_bulk, const MatrixView &abs_bulk)
get_scat_bulk_layer
Definition: disort.cc:770
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:212
int c_write_bad_var(int quiet, const char *varnam)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:833
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:523
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:86
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:279
thread_local Verbosity disort_verbosity
Definition: disort.cc:800
void get_angs(Vector &pfct_angs, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &Npfct)
get_angs.
Definition: disort.cc:565
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 ArrayOfStar &stars, 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 star_rte_los, const Index &gas_scattering_do, const Index &stars_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 star source.
Definition: disort.cc:1303
void get_pmom(Tensor3View pmom, ConstTensor3View pfct_bulk_par, ConstVectorView pfct_angs, const Index &Nlegendre)
get_pmom
Definition: disort.cc:684
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:341
void add_normed_phase_functions(Tensor3View &pfct1, const MatrixView &sca1, const MatrixView &pfct2, const MatrixView &sca2)
add_normed_phase_functions
Definition: disort.cc:62
constexpr Numeric COSMIC_BG_TEMP
Definition: disort.cc:60
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:448
constexpr Numeric PI
Definition: disort.cc:56
int c_write_too_small_dim(int quiet, const char *dimnam, int minval)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:851
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 ArrayOfStar &stars, 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 star_rte_los, const Index &gas_scattering_do, const Index &stars_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 star source.
Definition: disort.cc:938
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:539
#define abs(x)
#define min(a, b)
#define max(a, b)
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:354
Header file for logic.cc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:204
#define CREATE_OUT2
Definition: messages.h:205
#define CREATE_OUT0
Definition: messages.h:203
Index nelem(const Lines &l)
Number of lines.
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 k
Boltzmann constant convenience name [J/K].
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 Numeric e
Elementary charge convenience name [C].
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr Numeric dl(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
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:53
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
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.