ARTS 2.5.10 (git: 2f1c442c)
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 ARTS_USER_ERROR_IF(type == DS_ERROR, "DISORT ERROR >>> ", messag);
810
811 if (warning_limit) return;
812
813 if (++num_warnings <= MAX_WARNINGS) {
815 out1 << "DISORT WARNING >>> " << messag << "\n";
816 } else {
818 out1 << "DISORT TOO MANY WARNING MESSAGES -- They will no longer be "
819 "printed <<<<<<<\n\n";
820 warning_limit = TRUE;
821 }
822
823 return;
824}
825
826#undef MAX_WARNINGS
827
829int c_write_bad_var(int quiet, const char* varnam) {
830 const int maxmsg = 50;
831 static int nummsg = 0;
832
833 nummsg++;
834 if (quiet != QUIET) {
835 Verbosity verbosity = disort_verbosity;
837 out1 << " **** Input variable " << varnam << " in error ****\n";
838 if (nummsg == maxmsg) {
839 c_errmsg("Too many input errors. Aborting...", DS_ERROR);
840 }
841 }
842
843 return TRUE;
844}
845
847int c_write_too_small_dim(int quiet, const char* dimnam, int minval) {
848 if (quiet != QUIET) {
849 Verbosity verbosity = disort_verbosity;
851 out1 << " **** Symbolic dimension " << dimnam
852 << " should be increased to at least " << minval << " ****\n";
853 }
854
855 return TRUE;
856}
857
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(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 disort_state ds;
987 disort_output out;
988
989 if (quiet == 0)
990 disort_verbosity = verbosity;
991 else
992 disort_verbosity = Verbosity(0, 0, 0);
993
994 const Index nf = f_grid.nelem();
995
996 // solar dependent properties if no sun is present
997 // Number of azimuth angles
998 Index nphi = 1;
999 //local zenith angle of sun
1000 Numeric umu0 = 0.;
1001 //local azimuth angle of sun
1002 Numeric phi0 = 0.;
1003 //Intensity of incident sun beam
1004 Numeric fbeam = 0.;
1005
1006 if (suns_do) {
1007 nphi = aa_grid.nelem();
1008 umu0 = Conversion::cosd(sun_rte_los[0]);
1009 phi0 = sun_rte_los[1];
1010 if (phi0 < 0) {
1011 phi0 = phi0 + 360.;
1012 }
1013 }
1014
1015 ds.accur = 0.005;
1016 ds.flag.prnt[0] = FALSE;
1017 ds.flag.prnt[1] = FALSE;
1018 ds.flag.prnt[2] = FALSE;
1019 ds.flag.prnt[3] = FALSE;
1020 ds.flag.prnt[4] = TRUE;
1021
1022 ds.flag.usrtau = FALSE;
1023 ds.flag.usrang = TRUE;
1024 ds.flag.spher = FALSE;
1025 ds.flag.general_source = FALSE;
1026 ds.flag.output_uum = FALSE;
1027
1028 ds.nlyr = static_cast<int>(p.nelem() - 1);
1029
1030 ds.flag.brdf_type = BRDF_NONE;
1031
1032 ds.flag.ibcnd = GENERAL_BC;
1033 ds.flag.usrang = TRUE;
1034
1035 if (emission) {
1036 ds.flag.planck = TRUE;
1037 } else {
1038 ds.flag.planck = FALSE;
1039 }
1040 ds.flag.onlyfl = FALSE;
1041 ds.flag.lamber = TRUE;
1042 ds.flag.quiet = FALSE;
1043 if (intensity_correction) {
1044 ds.flag.intensity_correction = TRUE;
1045 ds.flag.old_intensity_correction = FALSE;
1046 } else {
1047 ds.flag.intensity_correction = FALSE;
1048 ds.flag.old_intensity_correction = FALSE;
1049 }
1050
1051 ds.nstr = static_cast<int>(nstreams);
1052 ds.nphase = ds.nstr;
1053 ds.nmom = ds.nstr;
1054 //ds.ntau = ds.nlyr + 1; // With ds.flag.usrtau = FALSE; set by cdisort
1055 ds.numu = static_cast<int>(za_grid.nelem());
1056 ds.nphi = static_cast<int>(nphi);
1057 Index Nlegendre = nstreams + 1;
1058
1059 /* Allocate memory */
1060 c_disort_state_alloc(&ds);
1061 c_disort_out_alloc(&ds, &out);
1062
1063 // Looking direction of solar beam
1064 ds.bc.umu0 = umu0;
1065 ds.bc.phi0 = phi0;
1066
1067 // Intensity of bottom-boundary isotropic illumination
1068 ds.bc.fluor = 0.;
1069
1070 // fill up azimuth angle and temperature array
1071 for (Index i = 0; i < ds.nphi; i++) ds.phi[i] = aa_grid[i];
1072
1073 if (ds.flag.planck==TRUE){
1074 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1075 }
1076
1077 // Transform to mu, starting with negative values
1078 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);
1079
1080 //gas absorption
1081 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1082 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1083
1084 // Get particle bulk properties
1085 Index nang;
1086 Vector pfct_angs;
1087 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1088 Index nf_ssd = scat_data[0][0].f_grid.nelem();
1089 Tensor3 pha_bulk_par;
1090
1091 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1092 nang = Npfct;
1093 nlinspace(pfct_angs, 0, 180, nang);
1094
1095 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1096
1097 ext_bulk_par = 0.0;
1098 abs_bulk_par = 0.0;
1099 pha_bulk_par = 0.0;
1100
1101 Index iflat = 0;
1102
1103 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1104 const Index nse = scat_data[iss].nelem();
1105 ext_abs_pfun_from_tro(ext_bulk_par,
1106 abs_bulk_par,
1107 pha_bulk_par,
1108 scat_data[iss],
1109 iss,
1110 pnd(Range(iflat, nse), joker),
1111 cboxlims,
1112 t,
1113 pfct_angs);
1114 iflat += nse;
1115 }
1116 } else {
1117 get_angs(pfct_angs, scat_data, Npfct);
1118 nang = pfct_angs.nelem();
1119
1120 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1121
1123 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1124 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1125 }
1126
1127 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1128 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1129
1130 // Legendre polynomials of phase function
1131 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1132 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1133
1134 if (gas_scattering_do) {
1135 // gas scattering
1136
1137 // layer averaged particle scattering coefficient
1138 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1139 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1140
1141 // call gas_scattering_properties
1142 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1143 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1144 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1145
1147 sca_coeff_gas_layer,
1148 sca_coeff_gas_level,
1149 pmom_gas,
1150 f_grid,
1151 p,
1152 t,
1153 vmr,
1154 gas_scattering_agenda);
1155
1156 // call add_norm_phase_functions
1158 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1159
1160 // add gas_scat_ext to ext_bulk_par
1161 ext_bulk_par += sca_coeff_gas_level;
1162 }
1163
1164 // Optical depth of layers
1165 Matrix dtauc(nf, ds.nlyr);
1166 // Single scattering albedo of layers
1167 Matrix ssalb(nf, ds.nlyr);
1168 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1169
1170 //upper boundary conditions:
1171 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1172 // emission. Both are applied additively.
1173 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1174 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1175 // are additive we should use either the one or the other.
1176 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1177 // Generally, temis!=1 should be avoided since that technically implies a
1178 // reflective upper boundary (though it seems that this is not exploited in
1179 // DISORT1.2, which we so far use).
1180
1181 // Cosmic background
1182 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1183 ds.bc.fisot = 0;
1184
1185 // Top of the atmosphere temperature and emissivity
1186 ds.bc.ttemp = COSMIC_BG_TEMP;
1187 ds.bc.btemp = surface_skin_t;
1188 ds.bc.temis = 1.;
1189
1190 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1191 snprintf(ds.header, 128, "ARTS Calc f_index = %" PRId64, f_index);
1192
1193 std::memcpy(ds.dtauc,
1194 dtauc(f_index, joker).get_c_array(),
1195 sizeof(Numeric) * ds.nlyr);
1196 std::memcpy(ds.ssalb,
1197 ssalb(f_index, joker).get_c_array(),
1198 sizeof(Numeric) * ds.nlyr);
1199
1200 // Wavenumber in [1/cm]
1201 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
1202 ds.wvnmhi += ds.wvnmhi * 1e-7;
1203 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1204
1205 // set
1206 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1207
1208 // Set irradiance of incident solar beam at top boundary
1209 if (suns_do) {
1210 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1211 (100 * SPEED_OF_LIGHT)*scale_factor;
1212 }
1213 ds.bc.fbeam = fbeam;
1214
1215 std::memcpy(ds.pmom,
1216 pmom(f_index, joker, joker).get_c_array(),
1217 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1218
1219 enum class Status { FIRST_TRY, RETRY, SUCCESS };
1220 Status tries = Status::FIRST_TRY;
1221 const Numeric eps = 2e-4; //two times the value defined in cdisort.c:3653
1222 do {
1223 try {
1224 c_disort(&ds, &out);
1225 tries = Status::SUCCESS;
1226 } catch (const std::runtime_error& e) {
1227 //catch cases if solar zenith angle=quadrature angle
1228 if (tries == Status::FIRST_TRY) {
1229 // change angle
1230 if (umu0 < 1 - eps) {
1231 umu0 += eps;
1232 } else if (umu0 > 1 - eps) {
1233 umu0 -= eps;
1234 }
1235
1236 const Numeric shift =
1237 abs(Conversion::acosd(umu0) - Conversion::acosd(ds.bc.umu0));
1239 out1
1240 << "Solar zenith angle coincided with one of the quadrature angles\n"
1241 << "We needed to shift the solar sun angle by " << shift
1242 << "deg.\n";
1243
1244 ds.bc.umu0 = umu0;
1245 tries = Status::RETRY;
1246 } else
1247 throw e;
1248 }
1249 } while (tries != Status::SUCCESS);
1250
1251 for (Index i = 0; i < ds.nphi; i++) {
1252 for (Index j = 0; j < ds.numu; j++) {
1253 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1254 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, i, 0) =
1255 out.uu[j + ((ds.nlyr - k - cboxlims[0]) + i * (ds.nlyr + 1)) *
1256 ds.numu] /
1257 (ds.wvnmhi - ds.wvnmlo) / (100 * SPEED_OF_LIGHT);
1258 }
1259 // To avoid potential numerical problems at interpolation of the field,
1260 // we copy the surface field to underground altitudes
1261 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1262 cloudbox_field(f_index, k, 0, 0, j, i, 0) =
1263 cloudbox_field(f_index, k + 1, 0, 0, j, i, 0);
1264 }
1265 }
1266 }
1267 }
1268
1269 // Allocate aux data
1270 disort_aux.resize(disort_aux_vars.nelem());
1271 // Allocate and set (if possible here) iy_aux
1272 Index cnt=-1;
1273 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1274
1275
1276 if (disort_aux_vars[i] == "Layer optical thickness"){
1277 cnt += 1;
1278 Matrix deltatau(nf, ds.nlyr, 0);
1279 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1280 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1281 deltatau(f_index, k - 1 + ncboxremoved) =
1282 dtauc(f_index, ds.nlyr - k + ncboxremoved);
1283 }
1284 }
1285 disort_aux[cnt] = deltatau;
1286 }
1287 else if (disort_aux_vars[i] == "Single scattering albedo"){
1288 cnt+=1;
1289 Matrix snglsctalbedo(nf, ds.nlyr,0);
1290 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1291 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1292 snglsctalbedo(f_index, k - 1 + ncboxremoved) =
1293 ssalb(f_index, ds.nlyr - k + ncboxremoved);
1294 }
1295 }
1296 disort_aux[cnt]=snglsctalbedo;
1297 }
1298 else if (disort_aux_vars[i] == "Direct beam") {
1299 cnt += 1;
1300 Matrix directbeam(nf, ds.nlyr + 1, 0);
1301 if (suns_do) {
1302 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1303 directbeam(f_index, cboxlims[1] - cboxlims[0] + ncboxremoved) =
1304 suns[0].spectrum(f_index, 0)/PI;
1305
1306 for (Index k = cboxlims[1] - cboxlims[0]; k > 0; k--) {
1307 directbeam(f_index, k - 1 + ncboxremoved) =
1308 directbeam(f_index, k + ncboxremoved) *
1309 exp(-dtauc(f_index, ds.nlyr - k + ncboxremoved)/umu0);
1310 }
1311 }
1312 }
1313 disort_aux[cnt]=directbeam;
1314 } else {
1316 "The only allowed strings in *disort_aux_vars* are:\n"
1317 " \"Layer optical thickness\"\n"
1318 " \"Single scattering albedo\"\n"
1319 " \"Direct beam\"\n"
1320 "but you have selected: \"", disort_aux_vars[i], "\"\n");
1321 }
1322 }
1323
1324 /* Free allocated memory */
1325 c_disort_out_free(&ds, &out);
1326 c_disort_state_free(&ds);
1327}
1328
1330 Tensor5& spectral_irradiance_field,
1331 ArrayOfMatrix& disort_aux,
1332 ConstVectorView f_grid,
1333 ConstVectorView p_grid,
1334 ConstVectorView z_profile,
1335 const Numeric& z_surface,
1336 ConstVectorView t_profile,
1337 ConstMatrixView vmr_profiles,
1338 ConstMatrixView pnd_profiles,
1339 const ArrayOfArrayOfSingleScatteringData& scat_data,
1340 const ArrayOfSun& suns,
1341 const Agenda& propmat_clearsky_agenda,
1342 const Agenda& gas_scattering_agenda,
1343 const ArrayOfIndex& cloudbox_limits,
1344 const Numeric& surface_skin_t,
1345 const Vector& surface_scalar_reflectivity,
1346 ConstVectorView sun_rte_los,
1347 const Index& gas_scattering_do,
1348 const Index& suns_do,
1349 const ArrayOfString& disort_aux_vars,
1350 const Numeric& scale_factor,
1351 const Index& nstreams,
1352 const Index& Npfct,
1353 const Index& only_tro,
1354 const Index& quiet,
1355 const Index& emission,
1356 const Index& intensity_correction,
1357 const Verbosity& verbosity) {
1358 // Create an atmosphere starting at z_surface
1359 Vector p, z, t;
1360 Matrix vmr, pnd;
1361 ArrayOfIndex cboxlims;
1362 Index ncboxremoved;
1363 //
1364 reduced_1datm(p,
1365 z,
1366 t,
1367 vmr,
1368 pnd,
1369 cboxlims,
1370 ncboxremoved,
1371 p_grid,
1372 z_profile,
1373 z_surface,
1374 t_profile,
1375 vmr_profiles,
1376 pnd_profiles,
1377 cloudbox_limits);
1378
1379 disort_state ds;
1380 disort_output out;
1381
1382 if (quiet == 0)
1383 disort_verbosity = verbosity;
1384 else
1385 disort_verbosity = Verbosity(0, 0, 0);
1386
1387 const Index nf = f_grid.nelem();
1388
1389 // solar dependent properties if no sun is present
1390 //local zenith angle of sun
1391 Numeric umu0 = 0.;
1392 //local azimuth angle of sun
1393 Numeric phi0 = 0.;
1394 //Intensity of incident sun beam
1395 Numeric fbeam = 0.;
1396
1397 if (suns_do) {
1398 umu0 = Conversion::cosd(sun_rte_los[0]);
1399 phi0 = sun_rte_los[1];
1400 if (phi0 < 0) {
1401 phi0 = phi0 + 360.;
1402 }
1403 }
1404
1405 ds.accur = 0.005;
1406 ds.flag.prnt[0] = FALSE;
1407 ds.flag.prnt[1] = FALSE;
1408 ds.flag.prnt[2] = FALSE;
1409 ds.flag.prnt[3] = FALSE;
1410 ds.flag.prnt[4] = TRUE;
1411
1412 ds.flag.usrtau = FALSE;
1413 ds.flag.usrang = TRUE;
1414 ds.flag.spher = FALSE;
1415 ds.flag.general_source = FALSE;
1416 ds.flag.output_uum = FALSE;
1417
1418 ds.nlyr = static_cast<int>(p.nelem() - 1);
1419
1420 ds.flag.brdf_type = BRDF_NONE;
1421
1422 ds.flag.ibcnd = GENERAL_BC;
1423 ds.flag.usrang = FALSE;
1424
1425 if (emission) {
1426 ds.flag.planck = TRUE;
1427 } else {
1428 ds.flag.planck = FALSE;
1429 }
1430 ds.flag.onlyfl = TRUE;
1431 ds.flag.lamber = TRUE;
1432 ds.flag.quiet = FALSE;
1433 if (intensity_correction) {
1434 ds.flag.intensity_correction = TRUE;
1435 ds.flag.old_intensity_correction = FALSE;
1436 } else {
1437 ds.flag.intensity_correction = FALSE;
1438 ds.flag.old_intensity_correction = FALSE;
1439 }
1440
1441 ds.nstr = static_cast<int>(nstreams);
1442 ds.nphase = ds.nstr;
1443 ds.nmom = ds.nstr;
1444
1445 // set grid size of angular variable to dummy value
1446 ds.numu = 1;
1447 ds.nphi = 1;
1448
1449 //Set number of legendre terms
1450 Index Nlegendre = nstreams + 1;
1451
1452 /* Allocate memory */
1453 c_disort_state_alloc(&ds);
1454 c_disort_out_alloc(&ds, &out);
1455
1456 // Looking direction of solar beam
1457 ds.bc.umu0 = umu0;
1458 ds.bc.phi0 = phi0;
1459
1460 // Intensity of bottom-boundary isotropic illumination
1461 ds.bc.fluor = 0.;
1462
1463 // fill up temperature array
1464 if (ds.flag.planck==TRUE){
1465 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
1466 }
1467
1468 //gas absorption
1469 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
1470 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
1471
1472 // Get particle bulk properties
1473 Index nang;
1474 Vector pfct_angs;
1475 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
1476 Index nf_ssd = scat_data[0][0].f_grid.nelem();
1477 Tensor3 pha_bulk_par;
1478
1479 if (only_tro && (Npfct < 0 || Npfct > 3)) {
1480 nang = Npfct;
1481 nlinspace(pfct_angs, 0, 180, nang);
1482
1483 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1484
1485 ext_bulk_par = 0.0;
1486 abs_bulk_par = 0.0;
1487 pha_bulk_par = 0.0;
1488
1489 Index iflat = 0;
1490
1491 for (Index iss = 0; iss < scat_data.nelem(); iss++) {
1492 const Index nse = scat_data[iss].nelem();
1493 ext_abs_pfun_from_tro(ext_bulk_par,
1494 abs_bulk_par,
1495 pha_bulk_par,
1496 scat_data[iss],
1497 iss,
1498 pnd(Range(iflat, nse), joker),
1499 cboxlims,
1500 t,
1501 pfct_angs);
1502 iflat += nse;
1503 }
1504 } else {
1505 get_angs(pfct_angs, scat_data, Npfct);
1506 nang = pfct_angs.nelem();
1507
1508 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
1509
1511 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
1512 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
1513 }
1514
1515 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
1516 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
1517
1518 // Legendre polynomials of phase function
1519 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
1520 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
1521
1522 if (gas_scattering_do) {
1523 // gas scattering
1524
1525 // layer averaged particle scattering coefficient
1526 Matrix sca_bulk_par_layer(nf_ssd, ds.nlyr);
1527 get_scat_bulk_layer(sca_bulk_par_layer, ext_bulk_par, abs_bulk_par);
1528
1529 // call gas_scattering_properties
1530 Matrix sca_coeff_gas_layer(nf_ssd, ds.nlyr, 0.);
1531 Matrix sca_coeff_gas_level(nf_ssd, ds.nlyr + 1, 0.);
1532 Matrix pmom_gas(ds.nlyr, Nlegendre, 0.);
1533
1535 sca_coeff_gas_layer,
1536 sca_coeff_gas_level,
1537 pmom_gas,
1538 f_grid,
1539 p,
1540 t,
1541 vmr,
1542 gas_scattering_agenda);
1543
1544 // call add_norm_phase_functions
1546 pmom, sca_bulk_par_layer, pmom_gas, sca_coeff_gas_layer);
1547
1548 // add gas_scat_ext to ext_bulk_par
1549 ext_bulk_par += sca_coeff_gas_level;
1550 }
1551
1552 // Optical depth of layers
1553 Matrix dtauc(nf, ds.nlyr);
1554 // Single scattering albedo of layers
1555 Matrix ssalb(nf, ds.nlyr);
1556 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
1557
1558 //upper boundary conditions:
1559 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
1560 // emission. Both are applied additively.
1561 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
1562 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
1563 // are additive we should use either the one or the other.
1564 // Note: previous setup (using fisot) setting temis=0 should be avoided.
1565 // Generally, temis!=1 should be avoided since that technically implies a
1566 // reflective upper boundary (though it seems that this is not exploited in
1567 // DISORT1.2, which we so far use).
1568
1569 // Cosmic background
1570 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
1571 ds.bc.fisot = 0;
1572
1573 // Top of the atmosphere temperature and emissivity
1574 ds.bc.ttemp = COSMIC_BG_TEMP;
1575 ds.bc.btemp = surface_skin_t;
1576 ds.bc.temis = 1.;
1577
1578 Matrix spectral_direct_irradiance_field;
1579 Matrix dFdtau(nf, ds.nlyr+1,0);
1580 Matrix deltatau(nf, ds.nlyr,0);
1581 Matrix snglsctalbedo(nf, ds.nlyr,0);
1582
1583 if (suns_do){
1584 //Resize direct field
1585 spectral_direct_irradiance_field.resize(nf, ds.nlyr+1);
1586 spectral_direct_irradiance_field = 0;
1587 }
1588
1589 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
1590 snprintf(ds.header, 128, "ARTS Calc f_index = %" PRId64, f_index);
1591
1592 std::memcpy(ds.dtauc,
1593 dtauc(f_index, joker).get_c_array(),
1594 sizeof(Numeric) * ds.nlyr);
1595 std::memcpy(ds.ssalb,
1596 ssalb(f_index, joker).get_c_array(),
1597 sizeof(Numeric) * ds.nlyr);
1598
1599 // Wavenumber in [1/cm]
1600 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
1601 ds.wvnmhi += ds.wvnmhi * 1e-7;
1602 ds.wvnmlo -= ds.wvnmlo * 1e-7;
1603
1604 // set
1605 ds.bc.albedo = surface_scalar_reflectivity[f_index];
1606
1607 // Set irradiance of incident solar beam at top boundary
1608 if (suns_do) {
1609 fbeam = suns[0].spectrum(f_index, 0)*(ds.wvnmhi - ds.wvnmlo)*
1610 (100 * SPEED_OF_LIGHT)*scale_factor;
1611 }
1612 ds.bc.fbeam = fbeam;
1613
1614 std::memcpy(ds.pmom,
1615 pmom(f_index, joker, joker).get_c_array(),
1616 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
1617
1618 c_disort(&ds, &out);
1619
1620 //factor for converting it into spectral radiance units
1621 const Numeric conv_fac=(ds.wvnmhi - ds.wvnmlo) * (100 * SPEED_OF_LIGHT);
1622
1623 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
1624 if (suns_do){
1625 // downward direct flux
1626 spectral_direct_irradiance_field(f_index, k + ncboxremoved) =
1627 -out.rad[ds.nlyr - k - cboxlims[0]].rfldir/conv_fac;
1628
1629 // downward total flux
1630 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1631 -(out.rad[ds.nlyr - k - cboxlims[0]].rfldir +
1632 out.rad[ds.nlyr - k - cboxlims[0]].rfldn)/conv_fac;
1633
1634 } else {
1635 // downward total flux
1636 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 0) =
1637 -out.rad[ds.nlyr - k - cboxlims[0]].rfldn/conv_fac;
1638 }
1639
1640 // upward flux
1641 spectral_irradiance_field(f_index, k + ncboxremoved, 0, 0, 1) =
1642 out.rad[ds.nlyr - k - cboxlims[0]].flup/conv_fac;
1643
1644 // flux divergence in tau space
1645 dFdtau(f_index, k + ncboxremoved) =
1646 -out.rad[ds.nlyr - k - cboxlims[0]].dfdt;
1647
1648 // k is running over the number of levels but deltatau, ssalb is defined for layers,
1649 // therefore we need to exlude k==0 and remove one from the index.
1650 if (k>0){
1651 deltatau(f_index, k - 1 + ncboxremoved) = ds.dtauc[ds.nlyr - k - 1 - cboxlims[0]];
1652 snglsctalbedo(f_index, k - 1 + ncboxremoved) = ds.ssalb[ds.nlyr - k - 1 - cboxlims[0]];
1653 }
1654 }
1655
1656 // To avoid potential numerical problems at interpolation of the field,
1657 // we copy the surface field to underground altitudes
1658 for (Index k = ncboxremoved - 1; k >= 0; k--) {
1659 spectral_irradiance_field(f_index, k, 0, 0, joker) =
1660 spectral_irradiance_field(f_index, k + 1, 0, 0, joker);
1661
1662 if (suns_do) {
1663 spectral_direct_irradiance_field(f_index, k) =
1664 spectral_direct_irradiance_field(f_index, k + 1);
1665 }
1666 dFdtau(f_index, k) = dFdtau(f_index, k + 1);
1667 }
1668 }
1669
1670 // Allocate aux data
1671 disort_aux.resize(disort_aux_vars.nelem());
1672 // Allocate and set (if possible here) iy_aux
1673 Index cnt=-1;
1674 for (Index i = 0; i < disort_aux_vars.nelem(); i++) {
1675
1676
1677 if (disort_aux_vars[i] == "Layer optical thickness"){
1678 cnt+=1;
1679 disort_aux[cnt]=deltatau;
1680 }
1681 else if (disort_aux_vars[i] == "Single scattering albedo"){
1682 cnt+=1;
1683 disort_aux[cnt]=snglsctalbedo;
1684 }
1685 else if (disort_aux_vars[i] == "Direct downward spectral irradiance") {
1686 cnt += 1;
1687 disort_aux[cnt] = spectral_direct_irradiance_field;
1688 }
1689 else if (disort_aux_vars[i] == "dFdtau") {
1690 cnt += 1;
1691 disort_aux[cnt] = dFdtau;
1692 }
1693 else {
1695 "The only allowed strings in *disort_aux_vars* are:\n"
1696 " \"Layer optical thickness\"\n"
1697 " \"Single scattering albedo\"\n"
1698 " \"Direct downward spectral irradiance\"\n"
1699 " \"dFdtau\"\n"
1700 "but you have selected: \"", disort_aux_vars[i], "\"\n");
1701 }
1702 }
1703
1704
1705
1706 /* Free allocated memory */
1707 c_disort_out_free(&ds, &out);
1708 c_disort_state_free(&ds);
1709}
1710
1712 //Output
1713 VectorView albedo,
1714 Numeric& btemp,
1715 //Input
1716 const Agenda& surface_rtprop_agenda,
1717 ConstVectorView f_grid,
1718 ConstVectorView scat_za_grid,
1719 const Numeric& surf_alt,
1720 const Verbosity& verbosity) {
1721 // Here, we derive an average surface albedo of the setup as given by ARTS'
1722 // surface_rtprop_agenda to use with Disorts's proprietary Lambertian surface.
1723 // In this way, ARTS-Disort can approximately mimick all surface reflection
1724 // types that ARTS itself can handle.
1725 // Surface temperature as derived from surface_rtprop_agenda is also returned.
1726 //
1727 // We derive the reflection matrices over all incident and reflected polar
1728 // angle directions and derive their integrated value (or weighted average).
1729 // The surface_rtprop_agenda handles one reflected direction at a time
1730 // (rtp_los) and for the reflected directions we loop over all (upwelling)
1731 // angles as given by scat_za_grid. The derived reflection matrices already
1732 // represent the reflectivity (or power reflection coefficient) for the given
1733 // reflected direction including proper angle weighting. For integrating/
1734 // averaging over the reflected directions, we use the same approach, i.e.
1735 // weight each angle by its associated range as given by the half distances to
1736 // the neighboring grid angles (using ARTS' scat_za_grid means 0/180deg are
1737 // grid points, 90deg shouldn't be among them (resulting from even number
1738 // requirement for Disort and the (implicitly assumed?) requirement of a
1739 // horizon-symmetric za_grid)).
1740 //
1741 // We do all frequencies here at once (assuming this is the faster variant as
1742 // the agenda anyway (can) provide output for full f_grid at once and as we
1743 // have to apply the same inter/extrapolation to all the frequencies).
1744
1746
1747 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1748
1749 const Index nf = f_grid.nelem();
1750 Index frza = 0;
1751 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
1752 if (frza == scat_za_grid.nelem()) {
1753 ostringstream os;
1754 os << "No upwelling direction found in scat_za_grid.\n";
1755 throw runtime_error(os.str());
1756 }
1757 const Index nrza = scat_za_grid.nelem() - frza;
1758 //cout << nrza << " upwelling directions found, starting from element #"
1759 // << frza << " of scat_za_grid.\n";
1760 Matrix dir_refl_coeff(nrza, nf, 0.);
1761
1762 // Local input of surface_rtprop_agenda.
1763 Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
1764
1765 // first derive the (reflected-)direction dependent power reflection
1766 // coefficient
1767 for (Index rza = 0; rza < nrza; rza++) {
1768 // Local output of surface_rtprop_agenda.
1769 Numeric surface_skin_t;
1770 Matrix surface_los;
1771 Tensor4 surface_rmatrix;
1772 Matrix surface_emission;
1773
1774 Vector rtp_los(1, scat_za_grid[rza + frza]);
1775 out2 << "Doing reflected dir #" << rza << " at " << rtp_los[0] << " degs\n";
1776
1778 surface_skin_t,
1779 surface_emission,
1780 surface_los,
1781 surface_rmatrix,
1782 f_grid,
1783 rtp_pos,
1784 rtp_los,
1785 surface_rtprop_agenda);
1786 if (surface_los.nrows() > 1) {
1787 ostringstream os;
1788 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1789 << "return zero or one direction in *surface_los*.";
1790 throw runtime_error(os.str());
1791 }
1792
1793 if (rza == 0)
1794 btemp = surface_skin_t;
1795 else if (surface_skin_t != btemp) {
1796 ostringstream os;
1797 os << "Something went wrong.\n"
1798 << " *surface_rtprop_agenda* returned different surface_skin_t\n"
1799 << " for different LOS.\n";
1800 throw runtime_error(os.str());
1801 }
1802 // Nothing to do if no surface_los (which means blackbody surface)
1803 // as dir_refl_coeff already set to 0
1804 if (!surface_los.empty()) {
1805 for (Index f_index = 0; f_index < nf; f_index++)
1806 dir_refl_coeff(rza, f_index) =
1807 surface_rmatrix(joker, f_index, 0, 0).sum();
1808 }
1809 out2 << " directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1810 << "\n";
1811 }
1812
1813 if (btemp < 0. || btemp > 1000.) {
1814 ostringstream os;
1815 os << "Surface temperature has been derived as " << btemp << " K,\n"
1816 << "which is not considered a meaningful value.\n";
1817 throw runtime_error(os.str());
1818 }
1819
1820 // now integrate/average the (reflected-)direction dependent power reflection
1821 // coefficients
1822 //
1823 // starting with deriving the angles defining the angle ranges
1824 Vector surf_int_grid(nrza + 1);
1825 // the first angle grid point should be around (but above) 90deg and should
1826 // cover the angle range between the 90deg and half-way point towards the next
1827 // angle grid point. we probably also want to check, that we don't
1828 // 'extrapolate' too much.
1829 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1830 ostringstream os;
1831 os << "Looks like scat_za_grid contains the 90deg direction,\n"
1832 << "which it shouldn't for running Disort.\n";
1833 throw runtime_error(os.str());
1834 }
1835 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1836 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1837 const Numeric ok_extrapol = 0.5;
1838 if ((za_extrapol - ok_extrapol) > 1e-6) {
1839 ostringstream os;
1840 os << "Extrapolation range from shallowest scat_za_grid point\n"
1841 << "to horizon is too big.\n"
1842 << " Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1843 << ", which is " << za_extrapol - 0.5 << " too big.\n";
1844 throw runtime_error(os.str());
1845 }
1847 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1848 ostringstream os;
1849 os << "Looks like last point in scat_za_grid is not 180deg.\n";
1850 throw runtime_error(os.str());
1851 }
1852
1853 surf_int_grid[0] = 90.;
1854 surf_int_grid[nrza] = 180.;
1855 for (Index rza = 1; rza < nrza; rza++)
1856 surf_int_grid[rza] =
1857 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1858 surf_int_grid *= DEG2RAD;
1859
1860 // now calculating the actual weights and apply them
1861 for (Index rza = 0; rza < nrza; rza++) {
1862 //Numeric coslow = cos(2.*surf_int_grid[rza]);
1863 //Numeric coshigh = cos(2.*surf_int_grid[rza+1]);
1864 //Numeric w = 0.5*(coshigh-coslow);
1865 Numeric w =
1866 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1867 //cout << "at reflLOS[" << rza << "]=" << scat_za_grid[frza+rza] << ":\n";
1868 //cout << " angle weight derives as w = 0.5*(" << coshigh << "-"
1869 // << coslow << ") = " << w << "\n";
1870 //cout << " weighting directional reflection coefficient from "
1871 // << dir_refl_coeff(rza,joker);
1872 dir_refl_coeff(rza, joker) *= w;
1873 //cout << " to " << dir_refl_coeff(rza,joker) << "\n";
1874 }
1875
1876 // eventually sum up the weighted directional power reflection coefficients
1877 for (Index f_index = 0; f_index < nf; f_index++) {
1878 albedo[f_index] = dir_refl_coeff(joker, f_index).sum();
1879 out2 << "at f=" << f_grid[f_index] * 1e-9
1880 << " GHz, ending up with albedo=" << albedo[f_index] << "\n";
1881 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1882 ostringstream os;
1883 os << "Something went wrong: Albedo must be inside [0,1],\n"
1884 << " but is not at freq #" << f_index << " , where it is "
1885 << albedo[f_index] << ".\n";
1886 throw runtime_error(os.str());
1887 }
1888 }
1889}
1890
1892 //Output
1893 VectorView albedo,
1894 Numeric& btemp,
1895 //Input
1896 const Agenda& surface_rtprop_agenda,
1897 ConstVectorView f_grid,
1898 const Numeric& surf_alt,
1899 const Numeric& inc_angle) {
1900
1901 Vector rtp_pos(1, surf_alt);
1902 Vector rtp_los(1, 180-inc_angle);
1903
1904 Numeric surface_skin_t;
1905 Matrix surface_los;
1906 Tensor4 surface_rmatrix;
1907 Matrix surface_emission;
1908
1910 surface_skin_t,
1911 surface_emission,
1912 surface_los,
1913 surface_rmatrix,
1914 f_grid,
1915 rtp_pos,
1916 rtp_los,
1917 surface_rtprop_agenda);
1918
1919 if (surface_los.nrows() > 1) {
1920 ostringstream os;
1921 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1922 << "return zero or one direction in *surface_los*.";
1923 throw runtime_error(os.str());
1924 }
1925 if (surface_los.empty()) {
1926 albedo = 0; // Blackbody assumed if no surface_los
1927 } else {
1928 albedo[joker] = surface_rmatrix(0,joker,0,0);
1929 }
1930
1931 btemp = surface_skin_t;
1932 if (btemp < 0. || btemp > 1000.) {
1933 ostringstream os;
1934 os << "Surface temperature has been derived as " << btemp << " K,\n"
1935 << "which is not considered a meaningful value.\n";
1936 throw runtime_error(os.str());
1937 }
1938}
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
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:25058
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:25803
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:25987
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:1065
bool empty() const noexcept
Definition: matpackI.h:1086
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
Index nbooks() const noexcept
Definition: matpackIV.h:143
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:160
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3View class.
Definition: matpackIII.h:251
The Tensor3 class.
Definition: matpackIII.h:352
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:435
The Tensor5 class.
Definition: matpackV.h:524
The Tensor6 class.
Definition: matpackVI.h:1105
The Tensor7 class.
Definition: matpackVII.h:2407
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:674
The Vector class.
Definition: matpackI.h:910
Workspace class.
Definition: workspace_ng.h:53
#define DEBUG_ONLY(...)
Definition: debug.h:71
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
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:1711
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:858
constexpr Numeric SPEED_OF_LIGHT
Definition: disort.cc:59
#define MAX_WARNINGS
Definition: disort.cc:802
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: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 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:1329
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:1891
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:829
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 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:847
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:544
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
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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
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: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.