ARTS 2.5.0 (git: 9ee3ac6c)
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
27#include "disort.h"
28#include <cmath>
29#include <cstdlib>
30#include <iostream>
31#include <stdexcept>
32#include "agenda_class.h"
33#include "array.h"
34#include "auto_md.h"
35#include "check_input.h"
36
37extern "C" {
38#include "cdisort.h"
39}
40
41#include "disort.h"
42#include "interpolation.h"
43#include "logic.h"
44#include "math_funcs.h"
45#include "messages.h"
46#include "rte.h"
47#include "xml_io.h"
48
49extern const Numeric PI;
50extern const Numeric DEG2RAD;
51extern const Numeric PLANCK_CONST;
52extern const Numeric SPEED_OF_LIGHT;
53extern const Numeric COSMIC_BG_TEMP;
54
55void check_disort_input( // Input
56 const Index& cloudbox_on,
57 const Index& atmfields_checked,
58 const Index& atmgeom_checked,
59 const Index& cloudbox_checked,
60 const Index& scat_data_checked,
61 const Index& atmosphere_dim,
62 const Index& stokes_dim,
63 const ArrayOfIndex& cloudbox_limits,
65 ConstVectorView za_grid,
66 const Index& nstreams) {
67 if (!cloudbox_on) {
68 throw runtime_error(
69 "Cloudbox is off, no scattering calculations to be"
70 "performed.");
71 }
72
73 if (atmfields_checked != 1)
74 throw runtime_error(
75 "The atmospheric fields must be flagged to have "
76 "passed a consistency check (atmfields_checked=1).");
77 if (atmgeom_checked != 1)
78 throw runtime_error(
79 "The atmospheric geometry must be flagged to have "
80 "passed a consistency check (atmgeom_checked=1).");
81 if (cloudbox_checked != 1)
82 throw runtime_error(
83 "The cloudbox must be flagged to have "
84 "passed a consistency check (cloudbox_checked=1).");
85 if (scat_data_checked != 1)
86 throw runtime_error(
87 "The scat_data must be flagged to have "
88 "passed a consistency check (scat_data_checked=1).");
89
90 if (atmosphere_dim != 1)
91 throw runtime_error(
92 "For running DISORT, atmospheric dimensionality "
93 "must be 1.\n");
94
95 if (stokes_dim < 0 || stokes_dim > 1)
96 throw runtime_error(
97 "For running DISORT, the dimension of stokes vector "
98 "must be 1.\n");
99
100 if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
101 throw runtime_error(
102 "*cloudbox_limits* is a vector which contains the"
103 "upper and lower limit of the cloud for all "
104 "atmospheric dimensions. So its dimension must"
105 "be 2 x *atmosphere_dim*");
106
107 if (cloudbox_limits[0] != 0) {
108 ostringstream os;
109 os << "DISORT calculations currently only possible with "
110 << "lower cloudbox limit\n"
111 << "at 0th atmospheric level "
112 << "(assumes surface there, ignoring z_surface).\n";
113 throw runtime_error(os.str());
114 }
115
116 if (scat_data.empty())
117 throw runtime_error(
118 "No single scattering data present.\n"
119 "See documentation of WSV *scat_data* for options to "
120 "make single scattering data available.\n");
121
122 // DISORT requires even number of streams:
123 // nstreams is total number of directions, up- and downwelling, and the up-
124 // and downwelling directions need to be symmetrically distributed, i.e. same
125 // number of directions in both hemispheres is required. horizontal direction
126 // (90deg) can not be covered in a plane-parallel atmosphere.
127 if (nstreams / 2 * 2 != nstreams) {
128 ostringstream os;
129 os << "DISORT requires an even number of streams, but yours is " << nstreams
130 << ".\n";
131 throw runtime_error(os.str());
132 }
133
134 // Zenith angle grid.
135 Index nza = za_grid.nelem();
136
137 // za_grid here is only relevant to provide an i_field from which the
138 // sensor los angles can be interpolated by yCalc; it does not the determine
139 // the accuracy of the DISORT output itself at these angles. So we can only
140 // apply a very rough test here, whether the grid is appropriate. However, we
141 // set the threshold fairly high since calculation costs for a higher number
142 // of angles are negligible.
143 if (nza < 20) {
144 ostringstream os;
145 os << "We require size of za_grid to be >= 20, to ensure a\n"
146 << "reasonable interpolation of the calculated cloudbox field.\n"
147 << "Note that for DISORT additional computation costs for\n"
148 << "larger numbers of angles are negligible.";
149 throw runtime_error(os.str());
150 }
151
152 if (za_grid[0] != 0. || za_grid[nza - 1] != 180.)
153 throw runtime_error("The range of *za_grid* must [0 180].");
154
155 if (!is_increasing(za_grid))
156 throw runtime_error("*za_grid* must be increasing.");
157
158 Index i = 1;
159 while (za_grid[i] <= 90) {
160 if (za_grid[i] == 90)
161 throw runtime_error("*za_grid* is not allowed to contain the value 90");
162 i++;
163 }
164
165 // DISORT can only handle randomly oriented particles.
166 bool all_totrand = true;
167 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
168 for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++)
169 if (scat_data[i_ss][i_se].ptype != PTYPE_TOTAL_RND) all_totrand = false;
170 if (!all_totrand) {
171 ostringstream os;
172 os << "DISORT can only handle scattering elements of type "
173 << PTYPE_TOTAL_RND << " (" << PTypeToString(PTYPE_TOTAL_RND) << "),\n"
174 << "but at least one element of other type (" << PTYPE_AZIMUTH_RND << "="
175 << PTypeToString(PTYPE_AZIMUTH_RND) << " or " << PTYPE_GENERAL << "="
176 << PTypeToString(PTYPE_GENERAL) << ") is present.\n";
177 throw runtime_error(os.str());
178 }
179}
180
181void init_ifield( // Output
182 Tensor7& cloudbox_field,
183 // Input
184 const Vector& f_grid,
185 const ArrayOfIndex& cloudbox_limits,
186 const Index& nang,
187 const Index& stokes_dim) {
188 const Index Nf = f_grid.nelem();
189 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
190 //const Index Nza = za_grid.nelem();
191
192 // Resize and initialize radiation field in the cloudbox
193 cloudbox_field.resize(Nf, Np_cloud, 1, 1, nang, 1, stokes_dim);
194 cloudbox_field = NAN;
195}
196
197void get_disortsurf_props( // Output
198 Vector& albedo,
199 Numeric& btemp,
200 // Input
201 ConstVectorView f_grid,
202 const Numeric& surface_skin_t,
203 ConstVectorView surface_scalar_reflectivity) {
204 // temperature of surface
205 if (surface_skin_t < 0. || surface_skin_t > 1000.) {
206 ostringstream os;
207 os << "Surface temperature has been set or derived as " << btemp << " K,\n"
208 << "which is not considered a meaningful value.\n"
209 << "For surface method 'L', *surface_skin_t* needs to\n"
210 << "be set and passed explicitly. Maybe you didn't do this?";
211 throw runtime_error(os.str());
212 }
213 btemp = surface_skin_t;
214
215 // surface albedo
216 if (surface_scalar_reflectivity.nelem() != f_grid.nelem() &&
217 surface_scalar_reflectivity.nelem() != 1) {
218 ostringstream os;
219 os << "The number of elements in *surface_scalar_reflectivity*\n"
220 << "should match length of *f_grid* or be 1."
221 << "\n length of *f_grid* : " << f_grid.nelem()
222 << "\n length of *surface_scalar_reflectivity* : "
223 << surface_scalar_reflectivity.nelem() << "\n";
224 throw runtime_error(os.str());
225 }
226
227 if (min(surface_scalar_reflectivity) < 0 ||
228 max(surface_scalar_reflectivity) > 1) {
229 throw runtime_error(
230 "All values in *surface_scalar_reflectivity*"
231 " must be inside [0,1].");
232 }
233
234 if (surface_scalar_reflectivity.nelem() > 1)
235 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
236 albedo[f_index] = surface_scalar_reflectivity[f_index];
237 else
238 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++)
239 albedo[f_index] = surface_scalar_reflectivity[0];
240}
241
243 MatrixView ext_bulk_gas,
244 const Agenda& propmat_clearsky_agenda,
245 ConstVectorView t_profile,
246 ConstMatrixView vmr_profiles,
247 ConstVectorView p_grid,
248 ConstVectorView f_grid) {
249 const Index Np = p_grid.nelem();
250
251 ARTS_ASSERT(ext_bulk_gas.nrows() == f_grid.nelem());
252 ARTS_ASSERT(ext_bulk_gas.ncols() == Np);
253
254 // Initialization
255 ext_bulk_gas = 0.;
256
257 // making gas property output containers and input dummies
258 const EnergyLevelMap rtp_nlte_dummy;
259 const Vector rtp_mag_dummy(3, 0);
260 const Vector ppath_los_dummy;
261 StokesVector nlte_dummy;
262 ArrayOfStokesVector partial_nlte_dummy;
263 ArrayOfPropagationMatrix partial_dummy;
264
265 PropagationMatrix propmat_clearsky_local;
266 for (Index ip = 0; ip < Np; ip++) {
268 propmat_clearsky_local,
269 nlte_dummy,
270 partial_dummy,
271 partial_nlte_dummy,
273 f_grid,
274 rtp_mag_dummy,
275 ppath_los_dummy,
276 p_grid[ip],
277 t_profile[ip],
278 rtp_nlte_dummy,
279 vmr_profiles(joker, ip),
280 propmat_clearsky_agenda);
281 ext_bulk_gas(joker, ip) += propmat_clearsky_local.Kjj();
282 }
283}
284
285void get_paroptprop(MatrixView ext_bulk_par,
286 MatrixView abs_bulk_par,
287 const ArrayOfArrayOfSingleScatteringData& scat_data,
288 ConstMatrixView pnd_profiles,
289 ConstVectorView t_profile,
291 const ArrayOfIndex& cloudbox_limits,
292 ConstVectorView f_grid) {
293 const Index Np_cloud = pnd_profiles.ncols();
294 DEBUG_ONLY(const Index Np = p_grid.nelem());
295 const Index nf = f_grid.nelem();
296
297 ARTS_ASSERT(ext_bulk_par.nrows() == nf);
298 ARTS_ASSERT(abs_bulk_par.nrows() == nf);
299 ARTS_ASSERT(ext_bulk_par.ncols() == Np);
300 ARTS_ASSERT(abs_bulk_par.ncols() == Np);
301
302 // Initialization
303 ext_bulk_par = 0.;
304 abs_bulk_par = 0.;
305
306 // preparing input data
307 Vector T_array = t_profile[Range(cloudbox_limits[0], Np_cloud)];
308 Matrix dir_array(1, 2, 0.); // just a dummy. only tot_random allowed, ie.
309 // optprop are independent of direction.
310
311 // making particle property output containers
312 ArrayOfArrayOfTensor5 ext_mat_Nse;
313 ArrayOfArrayOfTensor4 abs_vec_Nse;
314 ArrayOfArrayOfIndex ptypes_Nse;
315 Matrix t_ok;
316 ArrayOfTensor5 ext_mat_ssbulk;
317 ArrayOfTensor4 abs_vec_ssbulk;
318 ArrayOfIndex ptype_ssbulk;
319 Tensor5 ext_mat_bulk; //nf,nT,ndir,nst,nst
320 Tensor4 abs_vec_bulk;
321 Index ptype_bulk;
322
323 // calculate particle optical properties
324 opt_prop_NScatElems(ext_mat_Nse,
325 abs_vec_Nse,
326 ptypes_Nse,
327 t_ok,
328 scat_data,
329 1,
330 T_array,
331 dir_array,
332 -1);
333 opt_prop_ScatSpecBulk(ext_mat_ssbulk,
334 abs_vec_ssbulk,
335 ptype_ssbulk,
336 ext_mat_Nse,
337 abs_vec_Nse,
338 ptypes_Nse,
339 pnd_profiles,
340 t_ok);
341 opt_prop_Bulk(ext_mat_bulk,
342 abs_vec_bulk,
343 ptype_bulk,
344 ext_mat_ssbulk,
345 abs_vec_ssbulk,
346 ptype_ssbulk);
347
348 Index f_this = 0;
349 bool pf = (abs_vec_bulk.nbooks() != 1);
350 for (Index ip = 0; ip < Np_cloud; ip++)
351 for (Index f_index = 0; f_index < nf; f_index++) {
352 if (pf) f_this = f_index;
353 ext_bulk_par(f_index, ip + cloudbox_limits[0]) =
354 ext_mat_bulk(f_this, ip, 0, 0, 0);
355 abs_bulk_par(f_index, ip + cloudbox_limits[0]) =
356 abs_vec_bulk(f_this, ip, 0, 0);
357 }
358}
359
361 MatrixView ssalb,
362 ConstMatrixView ext_bulk_gas,
363 ConstMatrixView ext_bulk_par,
364 ConstMatrixView abs_bulk_par,
365 ConstVectorView z_profile) {
366 const Index nf = ext_bulk_gas.nrows();
367 const Index Np = ext_bulk_gas.ncols();
368
369 ARTS_ASSERT(dtauc.nrows() == nf);
370 ARTS_ASSERT(ssalb.nrows() == nf);
371 ARTS_ASSERT(dtauc.ncols() == Np - 1);
372 ARTS_ASSERT(ssalb.ncols() == Np - 1);
373
374 // Initialization
375 dtauc = 0.;
376 ssalb = 0.;
377
378 for (Index ip = 0; ip < Np - 1; ip++)
379 // Do layer averaging and derive single scattering albedo & optical depth
380 for (Index f_index = 0; f_index < nf; f_index++) {
381 Numeric ext =
382 0.5 * (ext_bulk_gas(f_index, ip) + ext_bulk_par(f_index, ip) +
383 ext_bulk_gas(f_index, ip + 1) + ext_bulk_par(f_index, ip + 1));
384 if (ext != 0) {
385 Numeric abs =
386 0.5 *
387 (ext_bulk_gas(f_index, ip) + abs_bulk_par(f_index, ip) +
388 ext_bulk_gas(f_index, ip + 1) + abs_bulk_par(f_index, ip + 1));
389 ssalb(f_index, Np - 2 - ip) = (ext - abs) / ext;
390 }
391
392 dtauc(f_index, Np - 2 - ip) = ext * (z_profile[ip + 1] - z_profile[ip]);
393 }
394}
395
396void get_angs(Vector& pfct_angs,
397 const ArrayOfArrayOfSingleScatteringData& scat_data,
398 const Index& Npfct) {
399 const Index min_nang = 3;
400 Index nang = Npfct;
401
402 if (Npfct < 0) {
403 Index this_ss = 0, this_se = 0;
404 // determine nang and pfct_angs from scat_data with finest za_grid
405 for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
406 for (Index i_se = scat_data[i_ss].nelem() - 1; i_se >= 0; i_se--)
407 // considering scat elems within one species mostly sorted from small to
408 // large sizes with large sizes corresponding to large za_grid. that is,
409 // starting searching from back should trigger if-clause (and variable
410 // update) less often.
411 if (nang < scat_data[i_ss][i_se].za_grid.nelem()) {
412 nang = scat_data[i_ss][i_se].za_grid.nelem();
413 this_ss = i_ss;
414 this_se = i_se;
415 }
416 pfct_angs = scat_data[this_ss][this_se].za_grid;
417 } else if (Npfct < min_nang) {
418 ostringstream os;
419 os << "Number of requested angular grid points (Npfct=" << Npfct
420 << ") is insufficient.\n"
421 << "At least " << min_nang << " points required.\n";
422 throw runtime_error(os.str());
423 } else {
424 nlinspace(pfct_angs, 0, 180, nang);
425 }
426}
427
428void get_parZ(Tensor3& pha_bulk_par,
429 const ArrayOfArrayOfSingleScatteringData& scat_data,
430 ConstMatrixView pnd_profiles,
431 ConstVectorView t_profile,
432 ConstVectorView pfct_angs,
433 const ArrayOfIndex& cloudbox_limits) {
434 const Index Np_cloud = pnd_profiles.ncols();
435 const Index nang = pfct_angs.nelem();
436
437 // Initialization
438 pha_bulk_par = 0.;
439
440 // preparing input data
441 Vector T_array = t_profile[Range(cloudbox_limits[0], Np_cloud)];
442 Matrix idir_array(1, 2, 0.); // we want pfct on sca ang grid, hence set
443 // pdir(*,0) to sca ang, all other to 0.
444 Matrix pdir_array(nang, 2, 0.);
445 pdir_array(joker, 0) = pfct_angs;
446
447 // making particle property output containers
448 ArrayOfArrayOfTensor6 pha_mat_Nse;
449 ArrayOfArrayOfIndex ptypes_Nse;
450 Matrix t_ok;
451 ArrayOfTensor6 pha_mat_ssbulk;
452 ArrayOfIndex ptype_ssbulk;
453 Tensor6 pha_mat_bulk; //nf,nT,npdir,nidir,nst,nst
454 Index ptype_bulk;
455
456 // calculate phase matrix
457 // FIXME: might be optimized by instead just executing pha_mat_1ScatElem where
458 // ext_bulk_par or pnd_field are non-zero.
459 pha_mat_NScatElems(pha_mat_Nse,
460 ptypes_Nse,
461 t_ok,
462 scat_data,
463 1,
464 T_array,
465 pdir_array,
466 idir_array,
467 -1);
468 pha_mat_ScatSpecBulk(pha_mat_ssbulk,
469 ptype_ssbulk,
470 pha_mat_Nse,
471 ptypes_Nse,
472 pnd_profiles,
473 t_ok);
474 pha_mat_Bulk(pha_mat_bulk, ptype_bulk, pha_mat_ssbulk, ptype_ssbulk);
475
476 pha_bulk_par(joker, Range(cloudbox_limits[0], Np_cloud), joker) =
477 pha_mat_bulk(joker, joker, joker, 0, 0, 0);
478}
479
480void get_pfct(Tensor3& pfct_bulk_par,
481 ConstTensor3View& pha_bulk_par,
482 ConstMatrixView ext_bulk_par,
483 ConstMatrixView abs_bulk_par,
484 const ArrayOfIndex& cloudbox_limits) {
485 const Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
486 const Index Np = pha_bulk_par.nrows();
487 const Index nf = pha_bulk_par.npages();
488 Index nang = pha_bulk_par.ncols();
489
490 ARTS_ASSERT(pfct_bulk_par.npages() == nf);
491 ARTS_ASSERT(pfct_bulk_par.nrows() == Np - 1);
492 ARTS_ASSERT(pfct_bulk_par.ncols() == nang);
493
494 // Initialization
495 pfct_bulk_par = 0.;
496
497 for (Index ip = cloudbox_limits[0]; ip < Np_cloud - 1; ip++)
498 for (Index f_index = 0; f_index < nf; f_index++) {
499 // Calculate layer averaged scattering (omitting factor 0.5 here as
500 // omitting it for layer averaged Z below, too).
501 Numeric sca =
502 (ext_bulk_par(f_index, ip) + ext_bulk_par(f_index, ip + 1)) -
503 (abs_bulk_par(f_index, ip) + abs_bulk_par(f_index, ip + 1));
504 if (sca != 0.) {
505 // Calculate layer averaged Z (omitting factor 0.5) and rescale from
506 // Z (Csca) to P (4Pi)
507 for (Index ia = 0; ia < nang; ia++)
508 pfct_bulk_par(f_index, Np - 2 - ip, ia) +=
509 pha_bulk_par(f_index, ip, ia) + pha_bulk_par(f_index, ip + 1, ia);
510 pfct_bulk_par(f_index, Np - 2 - ip, joker) *= 4 * PI / sca;
511 }
512 }
513}
514
516 ConstTensor3View pfct_bulk_par,
517 ConstVectorView pfct_angs,
518 const Index& Nlegendre) {
519 const Index nf = pfct_bulk_par.npages();
520 const Index nlyr = pfct_bulk_par.nrows();
521 const Index nang = pfct_bulk_par.ncols();
522
523 ARTS_ASSERT(nang == pfct_angs.nelem());
524
525 ARTS_ASSERT(pmom.npages() == nf);
526 ARTS_ASSERT(pmom.nrows() == nlyr);
527 ARTS_ASSERT(pmom.ncols() == Nlegendre);
528
529 Numeric pfct_threshold = 0.1;
530
531 // Initialization
532 pmom = 0.;
533
534 // we need the cosine of the pfct angles
535 Vector u(nang), adu(nang - 1);
536 Tensor3 px(nang - 1, Nlegendre, 2, 0.);
537 u[0] = cos(pfct_angs[0] * PI / 180.);
538 px(joker, 0, joker) = 1.;
539 for (Index ia = 1; ia < nang; ia++) {
540 u[ia] = cos(pfct_angs[ia] * PI / 180.);
541 adu[ia - 1] = abs(u[ia] - u[ia - 1]);
542 px(ia - 1, 1, 0) = u[ia - 1];
543 px(ia - 1, 1, 1) = u[ia];
544 for (Index l = 2; l < Nlegendre; l++) {
545 Numeric dl = (double)l;
546 px(ia - 1, l, 0) = (2 * dl - 1) / dl * u[ia - 1] * px(ia - 1, l - 1, 0) -
547 (dl - 1) / dl * px(ia - 1, l - 2, 0);
548 px(ia - 1, l, 1) = (2 * dl - 1) / dl * u[ia] * px(ia - 1, l - 1, 1) -
549 (dl - 1) / dl * px(ia - 1, l - 2, 1);
550 }
551 }
552
553 for (Index il = 0; il < nlyr; il++)
554 if (pfct_bulk_par(joker, il, 0).sum() != 0.)
555 for (Index f_index = 0; f_index < nf; f_index++) {
556 if (pfct_bulk_par(f_index, il, 0) != 0) {
557 Vector pfct = pfct_bulk_par(f_index, il, joker);
558
559 // Check if phase function is properly normalized
560 Numeric pint = 0.;
561 for (Index ia = 0; ia < nang - 1; ia++)
562 pint += 0.5 * adu[ia] * (pfct[ia] + pfct[ia + 1]);
563
564 if (abs(pint / 2. - 1.) > pfct_threshold) {
565 ostringstream os;
566 os << "Phase function normalization deviates from expected value by\n"
567 << 1e2 * pint / 2. - 1e2 << "(allowed: " << pfct_threshold * 1e2
568 << "%).\n"
569 << "Occurs at layer #" << il << " and frequency #" << f_index
570 << ".\n"
571 << "Something is wrong with your scattering data. Check!\n";
572 throw runtime_error(os.str());
573 }
574
575 // for the rest, rescale pfct to norm 2
576 pfct *= 2. / pint;
577
578 pmom(f_index, il, 0) = 1.;
579 for (Index ia = 0; ia < nang - 1; ia++) {
580 //for (Index l=0; l<Nlegendre; l++)
581 for (Index l = 1; l < Nlegendre; l++)
582 pmom(f_index, il, l) +=
583 0.25 * adu[ia] *
584 (px(ia, l, 0) * pfct[ia] + px(ia, l, 1) * pfct[ia + 1]);
585 }
586 }
587 }
588}
589
590// Use a thread_local variable to communicate the Verbosity to the
591// Disort error and warning functions. Ugly workaround, to avoid
592// passing a Verbosity argument throughout the whole cdisort code.
593// We want to avoid changes to the original code to keep it maintainable
594// in respect to upstream updates.
596
597#define MAX_WARNINGS 100
598
600void c_errmsg(const char* messag, int type) {
601 Verbosity verbosity = disort_verbosity;
602 static int warning_limit = FALSE, num_warnings = 0;
603
604 if (type == DS_ERROR) {
606 out0 << " ******* ERROR >>>>>> " << messag << "\n";
607 arts_exit(1);
608 }
609
610 if (warning_limit) return;
611
612 if (++num_warnings <= MAX_WARNINGS) {
614 out1 << " ******* WARNING >>>>>> " << messag << "\n";
615 } else {
617 out1 << " >>>>>> TOO MANY WARNING MESSAGES -- They will no longer be "
618 "printed <<<<<<<\n\n";
619 warning_limit = TRUE;
620 }
621
622 return;
623}
624
625#undef MAX_WARNINGS
626
628int c_write_bad_var(int quiet, const char* varnam) {
629 const int maxmsg = 50;
630 static int nummsg = 0;
631
632 nummsg++;
633 if (quiet != QUIET) {
634 Verbosity verbosity = disort_verbosity;
636 out1 << " **** Input variable " << varnam << " in error ****\n";
637 if (nummsg == maxmsg) {
638 c_errmsg("Too many input errors. Aborting...", DS_ERROR);
639 }
640 }
641
642 return TRUE;
643}
644
646int c_write_too_small_dim(int quiet, const char* dimnam, int minval) {
647 if (quiet != QUIET) {
648 Verbosity verbosity = disort_verbosity;
650 out1 << " **** Symbolic dimension " << dimnam
651 << " should be increased to at least " << minval << " ****\n";
652 }
653
654 return TRUE;
655}
656
658 Vector& z,
659 Vector& t,
660 Matrix& vmr,
661 Matrix& pnd,
662 ArrayOfIndex& cboxlims,
663 Index& ncboxremoved,
664 ConstVectorView p_grid,
665 ConstVectorView z_profile,
666 const Numeric& z_surface,
667 ConstVectorView t_profile,
668 ConstMatrixView vmr_profiles,
669 ConstMatrixView pnd_profiles,
670 const ArrayOfIndex& cloudbox_limits) {
671 // Surface at p_grid[0] and we just need to copy the original data
672 if (abs(z_surface - z_profile[0]) < 1e-3) {
673 p = p_grid;
674 z = z_profile;
675 t = t_profile;
676 vmr = vmr_profiles;
677 pnd = pnd_profiles;
678 cboxlims = cloudbox_limits;
679 ncboxremoved = 0;
680 }
681 // Surface above p_grid[0]
682 else {
683 // Some counters
684 Index np = p_grid.nelem(), ifirst = 0;
685 // Determine where to start with respect to z_profile
686 for (; z_surface >= z_profile[ifirst + 1]; ifirst++) {
687 }
688 np -= ifirst;
689 // Start by copying from ifirst to end
690 Range ind(ifirst, np);
691 p = p_grid[ind];
692 z = z_profile[ind];
693 t = t_profile[ind];
694 vmr = vmr_profiles(joker, ind);
695 // Insert surface altitude
696 z[0] = z_surface;
697 // Prepare interpolation
698 ArrayOfGridPos gp(1);
699 gridpos(gp[0], z_profile, z_surface);
700 Vector itw(2);
701 interpweights(itw, gp[0]);
702 // t and vmr
703 t[0] = interp(itw, t, gp[0]);
704 for (int i = 0; i < vmr.nrows(); i++) {
705 vmr(i, 0) = interp(itw, vmr(i, joker), gp[0]);
706 }
707 // p (we need a matrix version of iwt to use the function *itw2p*)
708 Matrix itw2(1, 2);
709 itw2(0, 0) = itw[0];
710 itw2(0, 1) = itw[1];
711 itw2p(p[0], p, gp, itw2);
712 // pnd_field and cloudbox limits need special treatment
713 cboxlims = cloudbox_limits;
714 if (ifirst < cloudbox_limits[0]) { // Surface below cloudbox
715 cboxlims[0] -= ifirst;
716 cboxlims[1] -= ifirst;
717 pnd = pnd_profiles;
718 ncboxremoved = 0;
719 } else { // Surface inside cloudbox
720 ncboxremoved = ifirst - cboxlims[0];
721 cboxlims[0] = 0;
722 cboxlims[1] = cloudbox_limits[1] - cloudbox_limits[0] - ncboxremoved;
723 ind = Range(ncboxremoved, cboxlims[1] + 1);
724 pnd = pnd_profiles(joker, ind);
725 gp[0].idx -= cloudbox_limits[0] + ncboxremoved;
726 for (int i = 0; i < pnd.nrows(); i++) {
727 pnd(i, 0) = interp(itw, pnd(i, joker), gp[0]);
728 }
729 }
730 }
731}
732
734 Tensor7& cloudbox_field,
735 ConstVectorView f_grid,
736 ConstVectorView p_grid,
737 ConstVectorView z_profile,
738 const Numeric& z_surface,
739 ConstVectorView t_profile,
740 ConstMatrixView vmr_profiles,
741 ConstMatrixView pnd_profiles,
742 const ArrayOfArrayOfSingleScatteringData& scat_data,
743 const Agenda& propmat_clearsky_agenda,
744 const ArrayOfIndex& cloudbox_limits,
745 const Numeric& surface_skin_t,
746 const Vector& surface_scalar_reflectivity,
747 ConstVectorView za_grid,
748 const Index& nstreams,
749 const Index& Npfct,
750 const Index& quiet,
751 const Verbosity& verbosity) {
752 // Create an atmosphere starting at z_surface
753 Vector p, z, t;
754 Matrix vmr, pnd;
755 ArrayOfIndex cboxlims;
756 Index ncboxremoved;
757 //
759 z,
760 t,
761 vmr,
762 pnd,
763 cboxlims,
764 ncboxremoved,
765 p_grid,
766 z_profile,
767 z_surface,
768 t_profile,
769 vmr_profiles,
770 pnd_profiles,
771 cloudbox_limits);
772
773 disort_state ds;
774 disort_output out;
775
776 if (quiet == 0)
777 disort_verbosity = verbosity;
778 else
779 disort_verbosity = Verbosity(0, 0, 0);
780
781 const Index nf = f_grid.nelem();
782
783 ds.accur = 0.005;
784 ds.flag.prnt[0] = FALSE;
785 ds.flag.prnt[1] = FALSE;
786 ds.flag.prnt[2] = FALSE;
787 ds.flag.prnt[3] = FALSE;
788 ds.flag.prnt[4] = TRUE;
789
790 ds.flag.usrtau = FALSE;
791 ds.flag.usrang = TRUE;
792 ds.flag.spher = FALSE;
793 ds.flag.general_source = FALSE;
794 ds.flag.output_uum = FALSE;
795
796 ds.nlyr = static_cast<int>(p.nelem() - 1);
797
798 ds.flag.brdf_type = BRDF_NONE;
799
800 ds.flag.ibcnd = GENERAL_BC;
801 ds.flag.usrang = TRUE;
802 ds.flag.planck = TRUE;
803 ds.flag.onlyfl = FALSE;
804 ds.flag.lamber = TRUE;
805 ds.flag.quiet = FALSE;
806 ds.flag.intensity_correction = TRUE;
807 ds.flag.old_intensity_correction = TRUE;
808
809 ds.nstr = static_cast<int>(nstreams);
810 ds.nphase = ds.nstr;
811 ds.nmom = ds.nstr;
812 //ds.ntau = ds.nlyr + 1; // With ds.flag.usrtau = FALSE; set by cdisort
813 ds.numu = static_cast<int>(za_grid.nelem());
814 ds.nphi = 1;
815 Index Nlegendre = nstreams + 1;
816
817 /* Allocate memory */
818 c_disort_state_alloc(&ds);
819 c_disort_out_alloc(&ds, &out);
820
821 // Properties of solar beam, set to zero as they are not needed
822 ds.bc.fbeam = 0.;
823 ds.bc.umu0 = 0.;
824 ds.bc.phi0 = 0.;
825 ds.bc.fluor = 0.;
826
827 // Since we have no solar source there is no angular dependance
828 ds.phi[0] = 0.;
829
830 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
831
832 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
833 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
834 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
836 ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
837
838 // Optical depth of layers
839 Matrix dtauc(nf, ds.nlyr);
840 // Single scattering albedo of layers
841 Matrix ssalb(nf, ds.nlyr);
842 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
843
844 // Transform to mu, starting with negative values
845 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);
846
847 //upper boundary conditions:
848 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
849 // emission. Both are applied additively.
850 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
851 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
852 // are additive we should use either the one or the other.
853 // Note: previous setup (using fisot) setting temis=0 should be avoided.
854 // Generally, temis!=1 should be avoided since that technically implies a
855 // reflective upper boundary (though it seems that this is not exploited in
856 // DISORT1.2, which we so far use).
857
858 // Cosmic background
859 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
860 ds.bc.fisot = 0;
861
862 // Top of the atmosphere temperature and emissivity
863 ds.bc.ttemp = COSMIC_BG_TEMP;
864 ds.bc.btemp = surface_skin_t;
865 ds.bc.temis = 1.;
866
867 Vector pfct_angs;
868 get_angs(pfct_angs, scat_data, Npfct);
869 Index nang = pfct_angs.nelem();
870
871 Index nf_ssd = scat_data[0][0].f_grid.nelem();
872 Tensor3 pha_bulk_par(nf_ssd, ds.nlyr + 1, nang);
873 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
874 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
875 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
876
877 // Legendre polynomials of phase function
878 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
879 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
880
881 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
882 sprintf(ds.header, "ARTS Calc f_index = %ld", f_index);
883
884 std::memcpy(ds.dtauc,
885 dtauc(f_index, joker).get_c_array(),
886 sizeof(Numeric) * ds.nlyr);
887 std::memcpy(ds.ssalb,
888 ssalb(f_index, joker).get_c_array(),
889 sizeof(Numeric) * ds.nlyr);
890
891 // Wavenumber in [1/cm]
892 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
893 ds.wvnmhi += ds.wvnmhi * 1e-7;
894 ds.wvnmlo -= ds.wvnmlo * 1e-7;
895
896 ds.bc.albedo = surface_scalar_reflectivity[f_index];
897
898 std::memcpy(ds.pmom,
899 pmom(f_index, joker, joker).get_c_array(),
900 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
901
902 c_disort(&ds, &out);
903
904 for (Index j = 0; j < ds.numu; j++) {
905 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
906 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, 0, 0) =
907 out.uu[ds.numu * (ds.nlyr - k - cboxlims[0]) + j] /
908 (ds.wvnmhi - ds.wvnmlo) / (100 * SPEED_OF_LIGHT);
909 }
910 // To avoid potential numerical problems at interpolation of the field,
911 // we copy the surface field to underground altitudes
912 for (Index k = ncboxremoved - 1; k >= 0; k--) {
913 cloudbox_field(f_index, k, 0, 0, j, 0, 0) =
914 cloudbox_field(f_index, k + 1, 0, 0, j, 0, 0);
915 }
916 }
917 }
918
919 /* Free allocated memory */
920 c_disort_out_free(&ds, &out);
921 c_disort_state_free(&ds);
922}
923
925 //Output
926 VectorView albedo,
927 Numeric& btemp,
928 //Input
929 const Agenda& surface_rtprop_agenda,
930 ConstVectorView f_grid,
931 ConstVectorView scat_za_grid,
932 const Numeric& surf_alt,
933 const Verbosity& verbosity) {
934 // Here, we derive an average surface albedo of the setup as given by ARTS'
935 // surface_rtprop_agenda to use with Disorts's proprietary Lambertian surface.
936 // In this way, ARTS-Disort can approximately mimick all surface reflection
937 // types that ARTS itself can handle.
938 // Surface temperature as derived from surface_rtprop_agenda is also returned.
939 //
940 // We derive the reflection matrices over all incident and reflected polar
941 // angle directions and derive their integrated value (or weighted average).
942 // The surface_rtprop_agenda handles one reflected direction at a time
943 // (rtp_los) and for the reflected directions we loop over all (upwelling)
944 // angles as given by scat_za_grid. The derived reflection matrices already
945 // represent the reflectivity (or power reflection coefficient) for the given
946 // reflected direction including proper angle weighting. For integrating/
947 // averaging over the reflected directions, we use the same approach, i.e.
948 // weight each angle by its associated range as given by the half distances to
949 // the neighboring grid angles (using ARTS' scat_za_grid means 0/180deg are
950 // grid points, 90deg shouldn't be among them (resulting from even number
951 // requirement for Disort and the (implicitly assumed?) requirement of a
952 // horizon-symmetric za_grid)).
953 //
954 // We do all frequencies here at once (assuming this is the faster variant as
955 // the agenda anyway (can) provide output for full f_grid at once and as we
956 // have to apply the same inter/extrapolation to all the frequencies).
957
959
960 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
961
962 const Index nf = f_grid.nelem();
963 Index frza = 0;
964 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
965 if (frza == scat_za_grid.nelem()) {
966 ostringstream os;
967 os << "No upwelling direction found in scat_za_grid.\n";
968 throw runtime_error(os.str());
969 }
970 const Index nrza = scat_za_grid.nelem() - frza;
971 //cout << nrza << " upwelling directions found, starting from element #"
972 // << frza << " of scat_za_grid.\n";
973 Matrix dir_refl_coeff(nrza, nf, 0.);
974
975 // Local input of surface_rtprop_agenda.
976 Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
977
978 // first derive the (reflected-)direction dependent power reflection
979 // coefficient
980 for (Index rza = 0; rza < nrza; rza++) {
981 // Local output of surface_rtprop_agenda.
982 Numeric surface_skin_t;
983 Matrix surface_los;
984 Tensor4 surface_rmatrix;
985 Matrix surface_emission;
986
987 Vector rtp_los(1, scat_za_grid[rza + frza]);
988 out2 << "Doing reflected dir #" << rza << " at " << rtp_los[0] << " degs\n";
989
991 surface_skin_t,
992 surface_emission,
993 surface_los,
994 surface_rmatrix,
995 f_grid,
996 rtp_pos,
997 rtp_los,
998 surface_rtprop_agenda);
999 //cout << "surf_los has " << surface_los.ncols() << " columns and "
1000 // << surface_los.nrows() << " rows.\n";
1001 ARTS_ASSERT(surface_los.ncols() == 1 || surface_los.nrows() == 0);
1002 if (rza == 0)
1003 btemp = surface_skin_t;
1004 else if (surface_skin_t != btemp) {
1005 ostringstream os;
1006 os << "Something went wrong.\n"
1007 << " *surface_rtprop_agenda* returned different surface_skin_t\n"
1008 << " for different LOS.\n";
1009 throw runtime_error(os.str());
1010 }
1011 if (surface_los.nrows() > 0) {
1012 for (Index f_index = 0; f_index < nf; f_index++)
1013 dir_refl_coeff(rza, f_index) =
1014 surface_rmatrix(joker, f_index, 0, 0).sum();
1015 }
1016 out2 << " directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1017 << "\n";
1018 }
1019
1020 if (btemp < 0. || btemp > 1000.) {
1021 ostringstream os;
1022 os << "Surface temperature has been derived as " << btemp << " K,\n"
1023 << "which is not considered a meaningful value.\n";
1024 throw runtime_error(os.str());
1025 }
1026
1027 // now integrate/average the (reflected-)direction dependent power reflection
1028 // coefficients
1029 //
1030 // starting with deriving the angles defining the angle ranges
1031 Vector surf_int_grid(nrza + 1);
1032 // the first angle grid point should be around (but above) 90deg and should
1033 // cover the angle range between the 90deg and half-way point towards the next
1034 // angle grid point. we probably also want to check, that we don't
1035 // 'extrapolate' too much.
1036 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1037 ostringstream os;
1038 os << "Looks like scat_za_grid contains the 90deg direction,\n"
1039 << "which it shouldn't for running Disort.\n";
1040 throw runtime_error(os.str());
1041 }
1042 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1043 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1044 const Numeric ok_extrapol = 0.5;
1045 if ((za_extrapol - ok_extrapol) > 1e-6) {
1046 ostringstream os;
1047 os << "Extrapolation range from shallowest scat_za_grid point\n"
1048 << "to horizon is too big.\n"
1049 << " Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1050 << ", which is " << za_extrapol - 0.5 << " too big.\n";
1051 throw runtime_error(os.str());
1052 }
1054 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1055 ostringstream os;
1056 os << "Looks like last point in scat_za_grid is not 180deg.\n";
1057 throw runtime_error(os.str());
1058 }
1059
1060 surf_int_grid[0] = 90.;
1061 surf_int_grid[nrza] = 180.;
1062 for (Index rza = 1; rza < nrza; rza++)
1063 surf_int_grid[rza] =
1064 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1065 surf_int_grid *= DEG2RAD;
1066
1067 // now calculating the actual weights and apply them
1068 for (Index rza = 0; rza < nrza; rza++) {
1069 //Numeric coslow = cos(2.*surf_int_grid[rza]);
1070 //Numeric coshigh = cos(2.*surf_int_grid[rza+1]);
1071 //Numeric w = 0.5*(coshigh-coslow);
1072 Numeric w =
1073 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1074 //cout << "at reflLOS[" << rza << "]=" << scat_za_grid[frza+rza] << ":\n";
1075 //cout << " angle weight derives as w = 0.5*(" << coshigh << "-"
1076 // << coslow << ") = " << w << "\n";
1077 //cout << " weighting directional reflection coefficient from "
1078 // << dir_refl_coeff(rza,joker);
1079 dir_refl_coeff(rza, joker) *= w;
1080 //cout << " to " << dir_refl_coeff(rza,joker) << "\n";
1081 }
1082
1083 // eventually sum up the weighted directional power reflection coefficients
1084 for (Index f_index = 0; f_index < nf; f_index++) {
1085 albedo[f_index] = dir_refl_coeff(joker, f_index).sum();
1086 out2 << "at f=" << f_grid[f_index] * 1e-9
1087 << " GHz, ending up with albedo=" << albedo[f_index] << "\n";
1088 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1089 ostringstream os;
1090 os << "Something went wrong: Albedo must be inside [0,1],\n"
1091 << " but is not at freq #" << f_index << " , where it is "
1092 << albedo[f_index] << ".\n";
1093 throw runtime_error(os.str());
1094 }
1095 }
1096}
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
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 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:22824
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:24392
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:44
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
A constant view of a Matrix.
Definition: matpackI.h:1014
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The MatrixView class.
Definition: matpackI.h:1125
The Matrix class.
Definition: matpackI.h:1225
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:165
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3View class.
Definition: matpackIII.h:239
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4 class.
Definition: matpackIV.h:421
The Tensor5 class.
Definition: matpackV.h:506
The Tensor6 class.
Definition: matpackVI.h:1088
The Tensor7 class.
Definition: matpackVII.h:2382
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5482
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
Workspace class.
Definition: workspace_ng.h:40
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
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)
Surface albed.
Definition: disort.cc:924
void c_errmsg(const char *messag, int type)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:600
const Numeric COSMIC_BG_TEMP
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:480
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:657
void run_cdisort(Workspace &ws, Tensor7 &cloudbox_field, 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 Agenda &propmat_clearsky_agenda, const ArrayOfIndex &cloudbox_limits, const Numeric &surface_skin_t, const Vector &surface_scalar_reflectivity, ConstVectorView za_grid, const Index &nstreams, const Index &Npfct, const Index &quiet, const Verbosity &verbosity)
Calculate doit_i_feild with Disort.
Definition: disort.cc:733
#define MAX_WARNINGS
Definition: disort.cc:597
void init_ifield(Tensor7 &cloudbox_field, const Vector &f_grid, const ArrayOfIndex &cloudbox_limits, const Index &nang, const Index &stokes_dim)
init_ifield.
Definition: disort.cc:181
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:428
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:197
int c_write_bad_var(int quiet, const char *varnam)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:628
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:360
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:55
thread_local Verbosity disort_verbosity
Definition: disort.cc:595
void get_angs(Vector &pfct_angs, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &Npfct)
get_angs.
Definition: disort.cc:396
const Numeric PI
void get_pmom(Tensor3View pmom, ConstTensor3View pfct_bulk_par, ConstVectorView pfct_angs, const Index &Nlegendre)
get_pmom
Definition: disort.cc:515
const Numeric SPEED_OF_LIGHT
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:242
const Numeric DEG2RAD
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:285
const Numeric PLANCK_CONST
int c_write_too_small_dim(int quiet, const char *dimnam, int minval)
Verbosity enabled replacement for the original cdisort function.
Definition: disort.cc:646
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:560
#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:215
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:351
Header file for logic.cc.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:205
#define CREATE_OUT2
Definition: messages.h:206
#define CREATE_OUT0
Definition: messages.h:204
Index nelem(const Lines &l)
Number of lines.
char Type type
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.
@ PTYPE_GENERAL
Definition: optproperties.h:53
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
#define u
#define w
Declaration of functions in rte.cc.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
This file contains basic functions to handle XML data files.