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