ARTS 2.5.4 (git: 4c0d3b4d)
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& only_tro,
751 const Index& quiet,
752 const Verbosity& verbosity) {
753
754 // Create an atmosphere starting at z_surface
755 Vector p, z, t;
756 Matrix vmr, pnd;
757 ArrayOfIndex cboxlims;
758 Index ncboxremoved;
759 //
761 z,
762 t,
763 vmr,
764 pnd,
765 cboxlims,
766 ncboxremoved,
767 p_grid,
768 z_profile,
769 z_surface,
770 t_profile,
771 vmr_profiles,
772 pnd_profiles,
773 cloudbox_limits);
774
775 disort_state ds;
776 disort_output out;
777
778 if (quiet == 0)
779 disort_verbosity = verbosity;
780 else
781 disort_verbosity = Verbosity(0, 0, 0);
782
783 const Index nf = f_grid.nelem();
784
785 ds.accur = 0.005;
786 ds.flag.prnt[0] = FALSE;
787 ds.flag.prnt[1] = FALSE;
788 ds.flag.prnt[2] = FALSE;
789 ds.flag.prnt[3] = FALSE;
790 ds.flag.prnt[4] = TRUE;
791
792 ds.flag.usrtau = FALSE;
793 ds.flag.usrang = TRUE;
794 ds.flag.spher = FALSE;
795 ds.flag.general_source = FALSE;
796 ds.flag.output_uum = FALSE;
797
798 ds.nlyr = static_cast<int>(p.nelem() - 1);
799
800 ds.flag.brdf_type = BRDF_NONE;
801
802 ds.flag.ibcnd = GENERAL_BC;
803 ds.flag.usrang = TRUE;
804 ds.flag.planck = TRUE;
805 ds.flag.onlyfl = FALSE;
806 ds.flag.lamber = TRUE;
807 ds.flag.quiet = FALSE;
808 ds.flag.intensity_correction = TRUE;
809 ds.flag.old_intensity_correction = TRUE;
810
811 ds.nstr = static_cast<int>(nstreams);
812 ds.nphase = ds.nstr;
813 ds.nmom = ds.nstr;
814 //ds.ntau = ds.nlyr + 1; // With ds.flag.usrtau = FALSE; set by cdisort
815 ds.numu = static_cast<int>(za_grid.nelem());
816 ds.nphi = 1;
817 Index Nlegendre = nstreams + 1;
818
819 /* Allocate memory */
820 c_disort_state_alloc(&ds);
821 c_disort_out_alloc(&ds, &out);
822
823 // Properties of solar beam, set to zero as they are not needed
824 ds.bc.fbeam = 0.;
825 ds.bc.umu0 = 0.;
826 ds.bc.phi0 = 0.;
827 ds.bc.fluor = 0.;
828
829 // Since we have no solar source there is no angular dependance
830 ds.phi[0] = 0.;
831
832 //upper boundary conditions:
833 // DISORT offers isotropic incoming radiance or emissivity-scaled planck
834 // emission. Both are applied additively.
835 // We want to have cosmic background radiation, for which ttemp=COSMIC_BG_TEMP
836 // and temis=1 should give identical results to fisot(COSMIC_BG_TEMP). As they
837 // are additive we should use either the one or the other.
838 // Note: previous setup (using fisot) setting temis=0 should be avoided.
839 // Generally, temis!=1 should be avoided since that technically implies a
840 // reflective upper boundary (though it seems that this is not exploited in
841 // DISORT1.2, which we so far use).
842
843 // Cosmic background
844 // we use temis*ttemp as upper boundary specification, hence CBR set to 0.
845 ds.bc.fisot = 0;
846
847 // Top of the atmosphere temperature and emissivity
848 ds.bc.ttemp = COSMIC_BG_TEMP;
849 ds.bc.btemp = surface_skin_t;
850 ds.bc.temis = 1.;
851
852
853 for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
854
855 // Absorption species
856 Matrix ext_bulk_gas(nf, ds.nlyr + 1);
857 get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
858
859 // Get particle bulk properties
860 Index nang;
861 Vector pfct_angs;
862 Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
863 Index nf_ssd = scat_data[0][0].f_grid.nelem();
864 Tensor3 pha_bulk_par;
865
866 if (only_tro && Npfct > 3) {
867 nang = Npfct;
868 nlinspace(pfct_angs, 0, 180, nang);
869
870 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
871
872 ext_bulk_par = 0.0;
873 abs_bulk_par = 0.0;
874 pha_bulk_par = 0.0;
875
876 Index iflat = 0;
877
878 for (Index iss=0; iss < scat_data.nelem(); iss++) {
879 const Index nse = scat_data[iss].nelem();
880 ext_abs_pfun_from_tro(ext_bulk_par,
881 abs_bulk_par,
882 pha_bulk_par,
883 scat_data[iss],
884 iss,
885 pnd(Range(iflat,nse),joker),
886 cboxlims,
887 t,
888 pfct_angs);
889 iflat += nse;
890 }
891 } else {
892 get_angs(pfct_angs, scat_data, Npfct);
893 nang = pfct_angs.nelem();
894
895 pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);
896
897 get_paroptprop(ext_bulk_par,
898 abs_bulk_par,
899 scat_data,
900 pnd,
901 t,
902 p,
903 cboxlims,
904 f_grid);
905 get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
906 }
907
908 // Optical depth of layers
909 Matrix dtauc(nf, ds.nlyr);
910 // Single scattering albedo of layers
911 Matrix ssalb(nf, ds.nlyr);
912 get_dtauc_ssalb(dtauc, ssalb, ext_bulk_gas, ext_bulk_par, abs_bulk_par, z);
913
914 // Transform to mu, starting with negative values
915 for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);
916
917 Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
918 get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);
919
920 // Legendre polynomials of phase function
921 Tensor3 pmom(nf_ssd, ds.nlyr, Nlegendre, 0.);
922 get_pmom(pmom, pfct_bulk_par, pfct_angs, Nlegendre);
923
924 for (Index f_index = 0; f_index < f_grid.nelem(); f_index++) {
925 sprintf(ds.header, "ARTS Calc f_index = %ld", f_index);
926
927 std::memcpy(ds.dtauc,
928 dtauc(f_index, joker).get_c_array(),
929 sizeof(Numeric) * ds.nlyr);
930 std::memcpy(ds.ssalb,
931 ssalb(f_index, joker).get_c_array(),
932 sizeof(Numeric) * ds.nlyr);
933
934 // Wavenumber in [1/cm]
935 ds.wvnmhi = ds.wvnmlo = (f_grid[f_index]) / (100. * SPEED_OF_LIGHT);
936 ds.wvnmhi += ds.wvnmhi * 1e-7;
937 ds.wvnmlo -= ds.wvnmlo * 1e-7;
938
939 ds.bc.albedo = surface_scalar_reflectivity[f_index];
940
941 std::memcpy(ds.pmom,
942 pmom(f_index, joker, joker).get_c_array(),
943 sizeof(Numeric) * pmom.nrows() * pmom.ncols());
944
945 c_disort(&ds, &out);
946
947 for (Index j = 0; j < ds.numu; j++) {
948 for (Index k = cboxlims[1] - cboxlims[0]; k >= 0; k--) {
949 cloudbox_field(f_index, k + ncboxremoved, 0, 0, j, 0, 0) =
950 out.uu[ds.numu * (ds.nlyr - k - cboxlims[0]) + j] /
951 (ds.wvnmhi - ds.wvnmlo) / (100 * SPEED_OF_LIGHT);
952 }
953 // To avoid potential numerical problems at interpolation of the field,
954 // we copy the surface field to underground altitudes
955 for (Index k = ncboxremoved - 1; k >= 0; k--) {
956 cloudbox_field(f_index, k, 0, 0, j, 0, 0) =
957 cloudbox_field(f_index, k + 1, 0, 0, j, 0, 0);
958 }
959 }
960 }
961
962 /* Free allocated memory */
963 c_disort_out_free(&ds, &out);
964 c_disort_state_free(&ds);
965}
966
968 //Output
969 VectorView albedo,
970 Numeric& btemp,
971 //Input
972 const Agenda& surface_rtprop_agenda,
973 ConstVectorView f_grid,
974 ConstVectorView scat_za_grid,
975 const Numeric& surf_alt,
976 const Verbosity& verbosity) {
977 // Here, we derive an average surface albedo of the setup as given by ARTS'
978 // surface_rtprop_agenda to use with Disorts's proprietary Lambertian surface.
979 // In this way, ARTS-Disort can approximately mimick all surface reflection
980 // types that ARTS itself can handle.
981 // Surface temperature as derived from surface_rtprop_agenda is also returned.
982 //
983 // We derive the reflection matrices over all incident and reflected polar
984 // angle directions and derive their integrated value (or weighted average).
985 // The surface_rtprop_agenda handles one reflected direction at a time
986 // (rtp_los) and for the reflected directions we loop over all (upwelling)
987 // angles as given by scat_za_grid. The derived reflection matrices already
988 // represent the reflectivity (or power reflection coefficient) for the given
989 // reflected direction including proper angle weighting. For integrating/
990 // averaging over the reflected directions, we use the same approach, i.e.
991 // weight each angle by its associated range as given by the half distances to
992 // the neighboring grid angles (using ARTS' scat_za_grid means 0/180deg are
993 // grid points, 90deg shouldn't be among them (resulting from even number
994 // requirement for Disort and the (implicitly assumed?) requirement of a
995 // horizon-symmetric za_grid)).
996 //
997 // We do all frequencies here at once (assuming this is the faster variant as
998 // the agenda anyway (can) provide output for full f_grid at once and as we
999 // have to apply the same inter/extrapolation to all the frequencies).
1000
1002
1003 chk_not_empty("surface_rtprop_agenda", surface_rtprop_agenda);
1004
1005 const Index nf = f_grid.nelem();
1006 Index frza = 0;
1007 while (frza < scat_za_grid.nelem() && scat_za_grid[frza] < 90.) frza++;
1008 if (frza == scat_za_grid.nelem()) {
1009 ostringstream os;
1010 os << "No upwelling direction found in scat_za_grid.\n";
1011 throw runtime_error(os.str());
1012 }
1013 const Index nrza = scat_za_grid.nelem() - frza;
1014 //cout << nrza << " upwelling directions found, starting from element #"
1015 // << frza << " of scat_za_grid.\n";
1016 Matrix dir_refl_coeff(nrza, nf, 0.);
1017
1018 // Local input of surface_rtprop_agenda.
1019 Vector rtp_pos(1, surf_alt); //atmosphere_dim is 1
1020
1021 // first derive the (reflected-)direction dependent power reflection
1022 // coefficient
1023 for (Index rza = 0; rza < nrza; rza++) {
1024 // Local output of surface_rtprop_agenda.
1025 Numeric surface_skin_t;
1026 Matrix surface_los;
1027 Tensor4 surface_rmatrix;
1028 Matrix surface_emission;
1029
1030 Vector rtp_los(1, scat_za_grid[rza + frza]);
1031 out2 << "Doing reflected dir #" << rza << " at " << rtp_los[0] << " degs\n";
1032
1034 surface_skin_t,
1035 surface_emission,
1036 surface_los,
1037 surface_rmatrix,
1038 f_grid,
1039 rtp_pos,
1040 rtp_los,
1041 surface_rtprop_agenda);
1042 if (surface_los.nrows() > 1) {
1043 ostringstream os;
1044 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1045 << "return zero or one direction in *surface_los*.";
1046 throw runtime_error(os.str());
1047 }
1048
1049 if (rza == 0)
1050 btemp = surface_skin_t;
1051 else if (surface_skin_t != btemp) {
1052 ostringstream os;
1053 os << "Something went wrong.\n"
1054 << " *surface_rtprop_agenda* returned different surface_skin_t\n"
1055 << " for different LOS.\n";
1056 throw runtime_error(os.str());
1057 }
1058 // Nothing to do if no surface_los (which means blackbody surface)
1059 // as dir_refl_coeff already set to 0
1060 if (!surface_los.empty()) {
1061 for (Index f_index = 0; f_index < nf; f_index++)
1062 dir_refl_coeff(rza, f_index) =
1063 surface_rmatrix(joker, f_index, 0, 0).sum();
1064 }
1065 out2 << " directional albedos[f_grid] = " << dir_refl_coeff(rza, joker)
1066 << "\n";
1067 }
1068
1069 if (btemp < 0. || btemp > 1000.) {
1070 ostringstream os;
1071 os << "Surface temperature has been derived as " << btemp << " K,\n"
1072 << "which is not considered a meaningful value.\n";
1073 throw runtime_error(os.str());
1074 }
1075
1076 // now integrate/average the (reflected-)direction dependent power reflection
1077 // coefficients
1078 //
1079 // starting with deriving the angles defining the angle ranges
1080 Vector surf_int_grid(nrza + 1);
1081 // the first angle grid point should be around (but above) 90deg and should
1082 // cover the angle range between the 90deg and half-way point towards the next
1083 // angle grid point. we probably also want to check, that we don't
1084 // 'extrapolate' too much.
1085 if (is_same_within_epsilon(scat_za_grid[frza], 90., 1e-6)) {
1086 ostringstream os;
1087 os << "Looks like scat_za_grid contains the 90deg direction,\n"
1088 << "which it shouldn't for running Disort.\n";
1089 throw runtime_error(os.str());
1090 }
1091 Numeric za_extrapol = (scat_za_grid[frza] - 90.) /
1092 (scat_za_grid[frza + 1] - scat_za_grid[frza]);
1093 const Numeric ok_extrapol = 0.5;
1094 if ((za_extrapol - ok_extrapol) > 1e-6) {
1095 ostringstream os;
1096 os << "Extrapolation range from shallowest scat_za_grid point\n"
1097 << "to horizon is too big.\n"
1098 << " Allowed extrapolation factor is 0.5.\n Yours is " << za_extrapol
1099 << ", which is " << za_extrapol - 0.5 << " too big.\n";
1100 throw runtime_error(os.str());
1101 }
1103 scat_za_grid[scat_za_grid.nelem() - 1], 180., 1e-6)) {
1104 ostringstream os;
1105 os << "Looks like last point in scat_za_grid is not 180deg.\n";
1106 throw runtime_error(os.str());
1107 }
1108
1109 surf_int_grid[0] = 90.;
1110 surf_int_grid[nrza] = 180.;
1111 for (Index rza = 1; rza < nrza; rza++)
1112 surf_int_grid[rza] =
1113 0.5 * (scat_za_grid[frza + rza - 1] + scat_za_grid[frza + rza]);
1114 surf_int_grid *= DEG2RAD;
1115
1116 // now calculating the actual weights and apply them
1117 for (Index rza = 0; rza < nrza; rza++) {
1118 //Numeric coslow = cos(2.*surf_int_grid[rza]);
1119 //Numeric coshigh = cos(2.*surf_int_grid[rza+1]);
1120 //Numeric w = 0.5*(coshigh-coslow);
1121 Numeric w =
1122 0.5 * (cos(2. * surf_int_grid[rza + 1]) - cos(2. * surf_int_grid[rza]));
1123 //cout << "at reflLOS[" << rza << "]=" << scat_za_grid[frza+rza] << ":\n";
1124 //cout << " angle weight derives as w = 0.5*(" << coshigh << "-"
1125 // << coslow << ") = " << w << "\n";
1126 //cout << " weighting directional reflection coefficient from "
1127 // << dir_refl_coeff(rza,joker);
1128 dir_refl_coeff(rza, joker) *= w;
1129 //cout << " to " << dir_refl_coeff(rza,joker) << "\n";
1130 }
1131
1132 // eventually sum up the weighted directional power reflection coefficients
1133 for (Index f_index = 0; f_index < nf; f_index++) {
1134 albedo[f_index] = dir_refl_coeff(joker, f_index).sum();
1135 out2 << "at f=" << f_grid[f_index] * 1e-9
1136 << " GHz, ending up with albedo=" << albedo[f_index] << "\n";
1137 if (albedo[f_index] < 0 || albedo[f_index] > 1.) {
1138 ostringstream os;
1139 os << "Something went wrong: Albedo must be inside [0,1],\n"
1140 << " but is not at freq #" << f_index << " , where it is "
1141 << albedo[f_index] << ".\n";
1142 throw runtime_error(os.str());
1143 }
1144 }
1145}
1146
1148 //Output
1149 VectorView albedo,
1150 Numeric& btemp,
1151 //Input
1152 const Agenda& surface_rtprop_agenda,
1153 ConstVectorView f_grid,
1154 const Numeric& surf_alt,
1155 const Numeric& inc_angle) {
1156
1157 Vector rtp_pos(1, surf_alt);
1158 Vector rtp_los(1, 180-inc_angle);
1159
1160 Numeric surface_skin_t;
1161 Matrix surface_los;
1162 Tensor4 surface_rmatrix;
1163 Matrix surface_emission;
1164
1166 surface_skin_t,
1167 surface_emission,
1168 surface_los,
1169 surface_rmatrix,
1170 f_grid,
1171 rtp_pos,
1172 rtp_los,
1173 surface_rtprop_agenda);
1174
1175 if (surface_los.nrows() > 1) {
1176 ostringstream os;
1177 os << "For this method, *surface_rtprop_agenda* must be set up to\n"
1178 << "return zero or one direction in *surface_los*.";
1179 throw runtime_error(os.str());
1180 }
1181 if (surface_los.empty()) {
1182 albedo = 0; // Blackbody assumed if no surface_los
1183 } else {
1184 albedo[joker] = surface_rmatrix(0,joker,0,0);
1185 }
1186
1187 btemp = surface_skin_t;
1188 if (btemp < 0. || btemp > 1000.) {
1189 ostringstream os;
1190 os << "Surface temperature has been derived as " << btemp << " K,\n"
1191 << "which is not considered a meaningful value.\n";
1192 throw runtime_error(os.str());
1193 }
1194}
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:23580
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:25148
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:51
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
A constant view of a Matrix.
Definition: matpackI.h:1050
bool empty() const noexcept
Definition: matpackI.h:1066
Index nrows() const noexcept
Definition: matpackI.h:1061
Index ncols() const noexcept
Definition: matpackI.h:1062
A constant view of a Tensor3.
Definition: matpackIII.h:130
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index nbooks() const noexcept
Definition: matpackIV.h:139
A constant view of a Vector.
Definition: matpackI.h:517
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
The MatrixView class.
Definition: matpackI.h:1169
The Matrix class.
Definition: matpackI.h:1270
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
The range class.
Definition: matpackI.h:166
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3View class.
Definition: matpackIII.h:244
The Tensor3 class.
Definition: matpackIII.h:344
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:427
The Tensor5 class.
Definition: matpackV.h:514
The Tensor6 class.
Definition: matpackVI.h:1097
The Tensor7 class.
Definition: matpackVII.h:2397
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:658
The Vector class.
Definition: matpackI.h:908
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)
surf_albedoCalc
Definition: disort.cc:967
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
#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
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:1147
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 &only_tro, const Index &quiet, const Verbosity &verbosity)
Calculate doit_i_feild with Disort.
Definition: disort.cc:733
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:545
#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
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:205
#define CREATE_OUT2
Definition: messages.h:206
#define CREATE_OUT0
Definition: messages.h:204
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr Numeric dl(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
void opt_prop_ScatSpecBulk(ArrayOfTensor5 &ext_mat, ArrayOfTensor4 &abs_vec, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor5 &ext_mat_se, const ArrayOfArrayOfTensor4 &abs_vec_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk extinction and absorption.
void pha_mat_ScatSpecBulk(ArrayOfTensor6 &pha_mat, ArrayOfIndex &ptype, const ArrayOfArrayOfTensor6 &pha_mat_se, const ArrayOfArrayOfIndex &ptypes_se, ConstMatrixView pnds, ConstMatrixView t_ok)
Scattering species bulk phase matrices.
void pha_mat_Bulk(Tensor6 &pha_mat, Index &ptype, const ArrayOfTensor6 &pha_mat_ss, const ArrayOfIndex &ptypes_ss)
Scattering species bulk phase matrix.
void opt_prop_Bulk(Tensor5 &ext_mat, Tensor4 &abs_vec, Index &ptype, const ArrayOfTensor5 &ext_mat_ss, const ArrayOfTensor4 &abs_vec_ss, const ArrayOfIndex &ptypes_ss)
one-line descript
void opt_prop_NScatElems(ArrayOfArrayOfTensor5 &ext_mat, ArrayOfArrayOfTensor4 &abs_vec, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &dir_array, const Index &f_index, const Index &t_interp_order)
Extinction and absorption from all scattering elements.
void pha_mat_NScatElems(ArrayOfArrayOfTensor6 &pha_mat, ArrayOfArrayOfIndex &ptypes, Matrix &t_ok, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &stokes_dim, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_index, const Index &t_interp_order)
Phase matrices from all scattering elements.
String PTypeToString(const PType &ptype)
Convert particle type enum value to String.
void ext_abs_pfun_from_tro(MatrixView ext_data, MatrixView abs_data, Tensor3View pfun_data, const ArrayOfSingleScatteringData &scat_data, const Index &iss, ConstMatrixView pnd_data, ArrayOfIndex &cloudbox_limits, ConstVectorView T_grid, ConstVectorView sa_grid)
Extinction, absorption and phase function for one scattering species, 1D and TRO.
@ PTYPE_GENERAL
Definition: optproperties.h:53
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
#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.