ARTS  2.3.1285(git:92a29ea9-dirty)
m_abs_lookup.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>
2 
3  This program is free software; you can redistribute it and/or
4  modify it under the terms of the GNU General Public License as
5  published by the Free Software Foundation; either version 2 of the
6  License, or (at your option) any 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 
26 #include <algorithm>
27 #include <limits>
28 #include <map>
29 
30 #include "absorption.h"
31 #include "agenda_class.h"
32 #include "arts.h"
33 #include "arts_omp.h"
34 #include "auto_md.h"
35 #include "check_input.h"
36 #include "cloudbox.h"
37 #include "gas_abs_lookup.h"
38 #include "global_data.h"
39 #include "interpolation_poly.h"
40 #include "math_funcs.h"
41 #include "matpackV.h"
42 #include "messages.h"
43 #include "physics_funcs.h"
44 #include "rng.h"
45 
46 extern const Index GFIELD4_FIELD_NAMES;
47 extern const Index GFIELD4_P_GRID;
48 
49 /* Workspace method: Doxygen documentation will be auto-generated */
50 void abs_lookupInit(GasAbsLookup& x, const Verbosity& verbosity) {
51  ArtsOut2 out2(verbosity);
52  // Nothing to do here.
53  // That means, we rely on the default constructor.
54 
55  x = GasAbsLookup();
56  out2 << " Created an empty gas absorption lookup table.\n";
57 }
58 
59 /* Workspace method: Doxygen documentation will be auto-generated */
60 void abs_lookupCalc( // Workspace reference:
61  Workspace& ws,
62  // WS Output:
63  GasAbsLookup& abs_lookup,
64  Index& abs_lookup_is_adapted,
65  // WS Input:
66  const ArrayOfArrayOfSpeciesTag& abs_species,
67  const ArrayOfArrayOfSpeciesTag& abs_nls,
68  const Vector& f_grid,
69  const Vector& abs_p,
70  const Matrix& abs_vmrs,
71  const Vector& abs_t,
72  const Vector& abs_t_pert,
73  const Vector& abs_nls_pert,
74  const Agenda& abs_xsec_agenda,
75  // Verbosity object:
76  const Verbosity& verbosity) {
79 
80  // We need this temporary variable to make a local copy of all VMRs,
81  // where we then perturb the H2O profile as needed
82  Matrix these_all_vmrs = abs_vmrs;
83 
84  // We will be calling an absorption agenda one species at a
85  // time. This is better than doing all simultaneously, because is
86  // saves memory and allows for consistent treatment of nonlinear
87  // species. But it means we need local copies of species, line list,
88  // and line shapes for agenda communication.
89 
90  // 1. Output of absorption calculations:
91 
92  // Absorption coefficients:
93  Matrix these_abs_coef;
94 
95  // Absorption cross sections per tag group.
96  ArrayOfMatrix abs_xsec_per_species, src_xsec_per_species;
97  ArrayOfArrayOfMatrix dabs_xsec_per_species_dx, dsrc_xsec_per_species_dx;
98 
99  // 2. Determine various important sizes:
100  const Index n_species = abs_species.nelem(); // Number of abs species
101  const Index n_nls = abs_nls.nelem(); // Number of nonlinear species
102  const Index n_f_grid = f_grid.nelem(); // Number of frequency grid points
103  const Index n_p_grid = abs_p.nelem(); // Number of presure grid points
104  const Index n_t_pert = abs_t_pert.nelem(); // Number of temp. perturbations
105  const Index n_nls_pert = abs_nls_pert.nelem(); // Number of VMR pert. for NLS
106 
107  // 3. Input to absorption calculations:
108 
109  // Absorption vmrs and temperature:
110  Matrix this_vmr(1, n_p_grid);
111  Vector abs_h2o(n_p_grid);
112  Vector this_t; // Has same dimension, but is
113  // initialized by assignment later.
114  const EnergyLevelMap this_nlte_dummy;
115 
116  // Species list, lines, and line shapes, all with only 1 element:
117  ArrayOfArrayOfSpeciesTag this_species(1);
118 
119  // List of active species for agenda call. Will always be filled with only
120  // one species.
121  ArrayOfIndex abs_species_active(1);
122 
123  // Local copy of nls_pert and t_pert:
124  Vector these_nls_pert; // Is resized every time when used
125  Vector these_t_pert; // Is resized later on
126 
127  // 4. Checks of input parameter correctness:
128 
129  const Index h2o_index = find_first_species_tg(
130  abs_species, species_index_from_species_name("H2O"));
131 
132  if (h2o_index < 0) {
133  // If there are nonlinear species, then at least one species must be
134  // H2O. We will use that to set h2o_abs, and to perturb in the case
135  // of nonlinear species.
136  if (n_nls > 0) {
137  ostringstream os;
138  os << "If you have nonlinear species, at least one\n"
139  << "species must be a H2O species.";
140  throw runtime_error(os.str());
141  } else {
142  out2 << " You have no H2O species. Absorption continua will not work.\n"
143  << " You should get a runtime error if you try them anyway.\n";
144  }
145  }
146 
147  // abs_species, f_grid, and p_grid should not be empty:
148  if (0 == n_species || 0 == n_f_grid || 0 == n_p_grid) {
149  ostringstream os;
150  os << "One of the required input variables is empty:\n"
151  << "abs_species.nelem() = " << n_species << ",\n"
152  << "f_grid.nelem() = " << n_f_grid << ",\n"
153  << "abs_p.nelem() = " << n_p_grid << ".";
154  throw runtime_error(os.str());
155  }
156 
157  // Set up the index array abs_nls from the tag array
158  // abs_nls. Give an error message if these
159  // tags are not included in abs_species.
160  ArrayOfIndex abs_nls_idx;
161  for (Index i = 0; i < n_nls; ++i) {
162  Index s;
163  for (s = 0; s < n_species; ++s) {
164  if (abs_nls[i] == abs_species[s]) {
165  abs_nls_idx.push_back(s);
166  break;
167  }
168  }
169  if (s == n_species) {
170  ostringstream os;
171  os << "Did not find *abs_nls* tag group \""
172  << get_tag_group_name(abs_nls[i]) << "\" in *abs_species*.";
173  throw runtime_error(os.str());
174  }
175  }
176 
177  // Furthermore abs_nls_idx must not contain duplicate values:
178  if (!is_unique(abs_nls_idx)) {
179  ostringstream os;
180  os << "The variable *abs_nls* must not have duplicate species.\n"
181  << "Value of *abs_nls*: " << abs_nls_idx;
182  throw runtime_error(os.str());
183  }
184 
185  // VMR matrix must match species list and pressure levels:
186  chk_size("abs_vmrs", abs_vmrs, n_species, n_p_grid);
187 
188  // Temperature vector must match number of pressure levels:
189  chk_size("abs_t", abs_t, n_p_grid);
190 
191  // abs_nls_pert should only be not-empty if we have nonlinear species:
192  if (((0 == n_nls) && (0 != n_nls_pert)) ||
193  ((0 != n_nls) && (0 == n_nls_pert))) {
194  ostringstream os;
195  os << "You have to set both abs_nls and abs_nls_pert, or none.";
196  throw runtime_error(os.str());
197  }
198 
199  // 4.a Set up a logical array for the nonlinear species.
200  ArrayOfIndex non_linear(n_species, 0);
201  for (Index s = 0; s < n_nls; ++s) {
202  non_linear[abs_nls_idx[s]] = 1;
203  }
204 
205  // 5. Set general lookup table properties:
206  abs_lookup.species = abs_species; // Species list
207  abs_lookup.nonlinear_species =
208  abs_nls_idx; // Nonlinear species (e.g., H2O, O2)
209  abs_lookup.f_grid = f_grid; // Frequency grid
210  abs_lookup.p_grid = abs_p; // Pressure grid
211  abs_lookup.vmrs_ref = abs_vmrs;
212  abs_lookup.t_ref = abs_t;
213  abs_lookup.t_pert = abs_t_pert;
214  abs_lookup.nls_pert = abs_nls_pert;
215 
216  // 5.a. Set log_p_grid:
217  abs_lookup.log_p_grid.resize(n_p_grid);
218  transform(abs_lookup.log_p_grid, log, abs_lookup.p_grid);
219 
220  // 6. Create abs_lookup.xsec with the right dimensions:
221  {
222  Index a, b, c, d;
223 
224  if (0 == n_t_pert)
225  a = 1;
226  else
227  a = n_t_pert;
228 
229  b = n_species + n_nls * (n_nls_pert - 1);
230 
231  c = n_f_grid;
232 
233  d = n_p_grid;
234 
235  abs_lookup.xsec.resize(a, b, c, d);
236  abs_lookup.xsec = NAN;
237  }
238 
239  // 6.a. Set up these_t_pert. This is done so that we can use the
240  // same loop over the perturbations, independent of
241  // whether we have temperature perturbations or not.
242  if (0 != n_t_pert) {
243  out2 << " With temperature perturbations.\n";
244  these_t_pert.resize(n_t_pert);
245  these_t_pert = abs_t_pert;
246  } else {
247  out2 << " No temperature perturbations.\n";
248  these_t_pert.resize(1);
249  these_t_pert = 0;
250  }
251 
252  const Index these_t_pert_nelem = these_t_pert.nelem();
253 
254  // 7. Now we have to fill abs_lookup.xsec with the right values!
255 
256  String fail_msg;
257  bool failed = false;
258 
259  // We have to make a local copy of the Workspace and the agenda because
260  // only non-reference types can be declared firstprivate in OpenMP
261  Workspace l_ws(ws);
262  Agenda l_abs_xsec_agenda(abs_xsec_agenda);
263 
264  // Loop species:
265  for (Index i = 0, spec = 0; i < n_species; ++i) {
266  // Skipping Zeeman and free_electrons species.
267  // (Mixed tag groups between those and other species are not allowed.)
268  if (is_zeeman(abs_species[i]) ||
269  abs_species[i][0].Type() == SpeciesTag::TYPE_FREE_ELECTRONS ||
270  abs_species[i][0].Type() == SpeciesTag::TYPE_PARTICLES) {
271  spec++;
272  continue;
273  }
274 
275  // spec is the index for the second dimension of abs_lookup.xsec.
276 
277  // Prepare absorption agenda input for this species:
278  out2 << " Doing species " << i + 1 << " of " << n_species << ": "
279  << abs_species[i] << ".\n";
280 
281  // Set active species:
282  abs_species_active[0] = i;
283 
284  // Get a dummy list of tag groups with only the current element:
285  this_species[0].resize(abs_species[i].nelem());
286  this_species[0] = abs_species[i];
287 
288  // Set up these_nls_pert. This is done so that we can use the
289  // same loop over the perturbations, independent of
290  // whether we have nonlinear species or not.
291  if (non_linear[i]) {
292  out2 << " This is a species with H2O VMR perturbations.\n";
293  these_nls_pert.resize(n_nls_pert);
294  these_nls_pert = abs_nls_pert;
295  } else {
296  these_nls_pert.resize(1);
297  these_nls_pert = 1;
298  }
299 
300  // Loop these_nls_pert:
301  for (Index s = 0; s < these_nls_pert.nelem(); ++s, ++spec) {
302  // Remember, spec is the index for the second dimension of
303  // abs_lookup.xsec
304 
305  if (non_linear[i]) {
306  out2 << " Doing H2O VMR variant " << s + 1 << " of " << n_nls_pert
307  << ": " << abs_nls_pert[s] << ".\n";
308  }
309 
310  // Make a local copy of the VMRs, and manipulate the H2O VMR within it.
311  // Note: We do not need a runtime error check that h2o_index is ok here,
312  // because earlier on we throw an error if there is no H2O species although we
313  // need it. So, if h2o_indes is -1, we here simply assume that there
314  // should not be a perturbation
315  if (h2o_index >= 0) {
316  these_all_vmrs(h2o_index, joker) = abs_vmrs(h2o_index, joker);
317  these_all_vmrs(h2o_index, joker) *=
318  these_nls_pert[s]; // Add perturbation
319  }
320 
321  // VMR for this species (still needed by interfact to continua):
322  // FIXME: This variable may go away eventually, when the continuum
323  // part no longer needs it.
324  this_vmr(0, joker) = these_all_vmrs(i, joker);
325 
326  // For abs_h2o, we can always add the perturbations (it will
327  // not make a difference if the species itself is also H2O).
328  // Attention, we need to treat here also the case that there
329  // is no H2O species. We will then set abs_h2o to
330  // -1. Absorption routines that do not really need abs_h2o
331  // will still run.
332  //
333  // FIXME: abs_h2o is currently still needed by the continuum part.
334  // Should go away eventually.
335  if (h2o_index == -1) {
336  // The case without H2O species.
337  abs_h2o.resize(1);
338  abs_h2o = -1;
339  } else {
340  // The normal case.
341  abs_h2o = these_all_vmrs(h2o_index, joker);
342  }
343 
344  // Loop temperature perturbations
345  // ------------------------------
346 
347  // We use a parallel for loop for this.
348 
349  // There is something strange here: abs_lookup seems to be
350  // "shared" by default, although I have set default(none). I
351  // suspect that the reason for this behavior is that
352  // abs_lookup is a return by reference parameter of this
353  // function. Anyway, shared is the correct setting for
354  // abs_lookup, so there is no problem.
355 
356 #pragma omp parallel for if ( \
357  !arts_omp_in_parallel() && \
358  these_t_pert_nelem >= \
359  arts_omp_get_max_threads()) private(this_t, \
360  abs_xsec_per_species, \
361  src_xsec_per_species, \
362  dabs_xsec_per_species_dx, \
363  dsrc_xsec_per_species_dx) \
364  firstprivate(l_ws, l_abs_xsec_agenda)
365  for (Index j = 0; j < these_t_pert_nelem; ++j) {
366  // Skip remaining iterations if an error occurred
367  if (failed) continue;
368 
369  // The try block here is necessary to correctly handle
370  // exceptions inside the parallel region.
371  try {
372  if (0 != n_t_pert) {
373  // We first prepare the output in a string here,
374  // so that we can write it to out3 with a single
375  // operation. This avoids messy output from
376  // multiple threads.
377  ostringstream os;
378 
379  os << " Doing temperature variant " << j + 1 << " of " << n_t_pert
380  << ": " << these_t_pert[j] << ".\n";
381 
382  out3 << os.str();
383  }
384 
385  // Create perturbed temperature profile:
386  this_t = abs_lookup.t_ref;
387  this_t += these_t_pert[j];
388 
389  // Call agenda to calculate absorption:
391  abs_xsec_per_species,
392  src_xsec_per_species,
393  dabs_xsec_per_species_dx,
394  dsrc_xsec_per_species_dx,
395  abs_species,
397  abs_species_active,
398  f_grid,
399  abs_p,
400  this_t,
401  this_nlte_dummy,
402  these_all_vmrs,
403  l_abs_xsec_agenda);
404 
405  // Store in the right place:
406  // Loop through all altitudes
407  for (Index p = 0; p < n_p_grid; ++p) {
408  abs_lookup.xsec(j, spec, Range(joker), p) =
409  abs_xsec_per_species[i](Range(joker), p);
410 
411  // There used to be a division by the number density
412  // n here. This is no longer necessary, since
413  // abs_xsec_per_species now contains true absorption
414  // cross sections.
415 
416  // IMPORTANT: There was a bug in my old Matlab
417  // function "create_lookup.m" to generate the lookup
418  // table. (The number density was always the
419  // reference one, and did not change for different
420  // temperatures.) Patricks Atmlab function
421  // "arts_abstable_from_arts1.m" did *not* have this bug.
422 
423  // Calculate the number density for the given pressure and
424  // temperature:
425  // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
426  // const Numeric n = number_density( abs_lookup.p_grid[p],
427  // this_t[p] );
428  // abs_lookup.xsec( j, spec, Range(joker), p ) /= n;
429  }
430  } // end of try block
431  catch (const std::runtime_error& e) {
432 #pragma omp critical(abs_lookupCalc_fail)
433  {
434  fail_msg = e.what();
435  failed = true;
436  }
437  }
438  } // end of parallel for loop
439 
440  if (failed) throw runtime_error(fail_msg);
441  }
442  }
443 
444  // 6. Initialize fgp_default.
445  abs_lookup.fgp_default.resize(f_grid.nelem());
446  gridpos_poly(abs_lookup.fgp_default, abs_lookup.f_grid, abs_lookup.f_grid, 0);
447 
448  // Set the abs_lookup_is_adapted flag. After all, the table fits the
449  // current frequency grid and species selection.
450  abs_lookup_is_adapted = 1;
451 }
452 
454 
472  const ArrayOfArrayOfSpeciesTag& abs_species,
473  const Verbosity& verbosity) {
474  CREATE_OUT3;
475 
476  cont.resize(0);
477 
478  // This is quite complicated, unfortunately. The approach here
479  // is copied from abs_xsec_per_speciesAddConts. For explanation,
480  // see there.
481 
482  // Loop tag groups:
483  for (Index i = 0; i < abs_species.nelem(); ++i) {
484  // Loop tags in tag group
485  for (Index s = 0; s < abs_species[i].nelem(); ++s) {
486  // Check for continuum tags
487  if (abs_species[i][s].Type() == SpeciesTag::TYPE_PREDEF ||
488  abs_species[i][s].Type() == SpeciesTag::TYPE_CIA) {
489  const String thisname = abs_species[i][s].Name();
490  // Ok, now we know this is a continuum tag.
491  out3 << " Continuum tag: " << thisname;
492 
493  // Check whether we want nonlinear treatment for
494  // this or not. We have three classes of continua:
495  // 1. Those that we know do not require it
496  // 2. Those that we know require h2o_abs
497  // 3. Those for which we do not know
498 
499  // The list here is not at all perfect, just a first
500  // guess. If you need particular models, you better
501  // check that the right thing is done with your model.
502 
503  // 1. Continua known to not use h2o_abs
504  // We take also H2O itself here, since this is
505  // handled separately
506  if (species_index_from_species_name("H2O") ==
507  abs_species[i][s].Species() ||
508  "N2-" == thisname.substr(0, 3) || "CO2-" == thisname.substr(0, 4) ||
509  "O2-CIA" == thisname.substr(0, 6) ||
510  "O2-v0v" == thisname.substr(0, 6) ||
511  "O2-v1v" == thisname.substr(0, 6) ||
512  "H2-CIA" == thisname.substr(0, 6) ||
513  "He-CIA" == thisname.substr(0, 6) ||
514  "CH4-CIA" == thisname.substr(0, 7) ||
515  "liquidcloud-MPM93" == thisname.substr(0, 17) ||
516  "liquidcloud-ELL07" == thisname.substr(0, 17)) {
517  out3 << " --> not added.\n";
518  break;
519  }
520 
521  // 2. Continua known to use h2o_abs
522  if ("O2-" == thisname.substr(0, 3)) {
523  cont.push_back(i);
524  out3 << " --> added to abs_nls.\n";
525  break;
526  }
527 
528  // 3. abs_species tags that are NOT allowed in LUT
529  // calculations
530  if ("icecloud-" == thisname.substr(0, 9) ||
531  "rain-" == thisname.substr(0, 5)) {
532  ostringstream os;
533  os << "Tag " << thisname << " not allowed in absorption "
534  << "lookup tables.";
535  throw runtime_error(os.str());
536  }
537 
538  // If we get here, then the tag was neither in the
539  // posivitive nor in the negative list. We through a
540  // runtime error.
541  out3 << " --> unknown.\n";
542  ostringstream os;
543  os << "Unknown whether tag " << thisname
544  << " is a nonlinear species (i.e. uses h2o_abs) or not.\n"
545  << "Cannot set abs_nls automatically.";
546  throw runtime_error(os.str());
547  }
548  }
549  }
550 }
551 
553 
562  const ArrayOfArrayOfSpeciesTag& abs_species,
563  const Verbosity& verbosity) {
564  CREATE_OUT2;
565 
566  abs_nls.resize(0);
567 
568  // Add all H2O species as non-linear:
569  Index next_h2o = 0;
570  while (-1 !=
571  (next_h2o = find_next_species_tg(
572  abs_species, species_index_from_species_name("H2O"), next_h2o))) {
573  abs_nls.push_back(abs_species[next_h2o]);
574  ++next_h2o;
575  }
576 
577  // Certain continuum models also depend on abs_h2o. There is a
578  // helper function that contains a list of these.
579  ArrayOfIndex cont;
580  find_nonlinear_continua(cont, abs_species, verbosity);
581 
582  // Add these to abs_nls:
583  for (Index i = 0; i < cont.nelem(); ++i) {
584  abs_nls.push_back(abs_species[cont[i]]);
585  }
586 
587  out2 << " Species marked for nonlinear treatment:\n";
588  for (Index i = 0; i < abs_nls.nelem(); ++i) {
589  out2 << " ";
590  for (Index j = 0; j < abs_nls[i].nelem(); ++j) {
591  if (j != 0) out2 << ", ";
592  out2 << abs_nls[i][j].Name();
593  }
594  out2 << "\n";
595  }
596 }
597 
599 
612 void choose_abs_t_pert(Vector& abs_t_pert,
613  ConstVectorView abs_t,
614  ConstVectorView tmin,
615  ConstVectorView tmax,
616  const Numeric& step,
617  const Index& p_interp_order,
618  const Index& t_interp_order,
619  const Verbosity& verbosity) {
620  CREATE_OUT2;
621  CREATE_OUT3;
622 
623  // The code to find out the range for perturbation is a bit
624  // complicated. The problem is that, since we use higher order
625  // interpolation for p, we may require temperatures well outside the
626  // local min/max range at the reference level. We solve this by
627  // really looking at the min/max range for those levels required by
628  // the p_interp_order.
629 
630  Numeric mindev = 1e9;
631  Numeric maxdev = -1e9;
632 
633  Vector the_grid(0, abs_t.nelem(), 1);
634  for (Index i = 0; i < the_grid.nelem(); ++i) {
635  GridPosPoly gp;
636  gridpos_poly(gp, the_grid, (Numeric)i, p_interp_order);
637 
638  for (Index j = 0; j < gp.idx.nelem(); ++j) {
639  // Our pressure grid for the lookup table may be coarser than
640  // the original one for the batch cases. This may lead to max/min
641  // values in the original data that exceed those we assumed
642  // above. We add some extra margin here to account for
643  // that. (The margin is +-10K)
644 
645  Numeric delta_min = tmin[i] - abs_t[gp.idx[j]] - 10;
646  Numeric delta_max = tmax[i] - abs_t[gp.idx[j]] + 10;
647 
648  if (delta_min < mindev) mindev = delta_min;
649  if (delta_max > maxdev) maxdev = delta_max;
650  }
651  }
652 
653  out3 << " abs_t_pert: mindev/maxdev : " << mindev << " / " << maxdev << "\n";
654 
655  // We divide the interval between mindev and maxdev, so that the
656  // steps are of size *step* or smaller. But we also need at least
657  // *t_interp_order*+1 points.
658  Index div = t_interp_order;
659  Numeric effective_step;
660  do {
661  effective_step = (maxdev - mindev) / (Numeric)div;
662  ++div;
663  } while (effective_step > step);
664 
665  abs_t_pert = Vector(mindev, div, effective_step);
666 
667  out2 << " abs_t_pert: " << abs_t_pert[0] << " K to "
668  << abs_t_pert[abs_t_pert.nelem() - 1] << " K in steps of "
669  << effective_step << " K (" << abs_t_pert.nelem() << " grid points)\n";
670 }
671 
673 
686 void choose_abs_nls_pert(Vector& abs_nls_pert,
687  ConstVectorView refprof,
688  ConstVectorView minprof,
689  ConstVectorView maxprof,
690  const Numeric& step,
691  const Index& p_interp_order,
692  const Index& nls_interp_order,
693  const Verbosity& verbosity) {
694  CREATE_OUT2;
695  CREATE_OUT3;
696 
697  // The code to find out the range for perturbation is a bit
698  // complicated. The problem is that, since we use higher order
699  // interpolation for p, we may require humidities well outside the
700  // local min/max range at the reference level. We solve this by
701  // really looking at the min/max range for those levels required by
702  // the p_interp_order.
703 
704  Numeric mindev = 0;
705  Numeric maxdev = -1e9;
706 
707  // mindev is set to zero from the start, since we always want to
708  // include 0.
709 
710  Vector the_grid(0, refprof.nelem(), 1);
711  for (Index i = 0; i < the_grid.nelem(); ++i) {
712  // cout << "!!!!!! i = " << i << "\n";
713  // cout << " min/ref/max = " << minprof[i] << " / "
714  // << refprof[i] << " / "
715  // << maxprof[i] << "\n";
716 
717  GridPosPoly gp;
718  gridpos_poly(gp, the_grid, (Numeric)i, p_interp_order);
719 
720  for (Index j = 0; j < gp.idx.nelem(); ++j) {
721  // cout << "!!!!!! j = " << j << "\n";
722  // cout << " ref[j] = " << refprof[gp.idx[j]] << " ";
723  // cout << " minfrac[j] = " << minprof[i] / refprof[gp.idx[j]] << " ";
724  // cout << " maxfrac[j] = " << maxprof[i] / refprof[gp.idx[j]] << " \n";
725 
726  // Our pressure grid for the lookup table may be coarser than
727  // the original one for the batch cases. This may lead to max/min
728  // values in the original data that exceed those we assumed
729  // above. We add some extra margin to the max value here to account for
730  // that. (The margin is a factor of 2.)
731 
732  Numeric delta_min = minprof[i] / refprof[gp.idx[j]];
733  Numeric delta_max = 2 * maxprof[i] / refprof[gp.idx[j]];
734 
735  if (delta_min < mindev) mindev = delta_min;
736  // do not update maxdev, when delta_max is infinity (this results from
737  // refprof being 0)
738  if (!std::isinf(delta_max) && (delta_max > maxdev)) maxdev = delta_max;
739  }
740  }
741 
742  out3 << " abs_nls_pert: mindev/maxdev : " << mindev << " / " << maxdev
743  << "\n";
744 
745  bool allownegative = false;
746  if (mindev < 0) {
747  out2
748  << " Warning: I am getting a negative fractional distance to the H2O\n"
749  << " reference profile. Some of your H2O fields may contain negative values.\n"
750  << " Will allow negative values also for abs_nls_pert.\n";
751  allownegative = true;
752  }
753 
754  if (!allownegative) {
755  mindev = 0;
756  out3 << " Adjusted mindev : " << mindev << "\n";
757  }
758 
759  if (std::isinf(maxdev)) {
760  ostringstream os;
761  os << "Perturbation upper limit is infinity (likely due to the reference\n"
762  << "profile being 0 at at least one pressure level). Can not work\n"
763  << "with that.";
764  throw runtime_error(os.str());
765  }
766 
767  // We divide the interval between mindev and maxdev, so that the
768  // steps are of size *step* or smaller. But we also need at least
769  // *nls_interp_order*+1 points.
770  Index div = nls_interp_order;
771  Numeric effective_step;
772  do {
773  effective_step = (maxdev - mindev) / (Numeric)div;
774  ++div;
775  } while (effective_step > step);
776 
777  abs_nls_pert = Vector(mindev, div, effective_step);
778 
779  // If there are negative values, we also add 0. The reason for this
780  // is that 0 is a turning point.
781  if (allownegative) {
782  VectorInsertGridPoints(abs_nls_pert, abs_nls_pert, {0}, verbosity);
783  out2
784  << " I am including also 0 in the abs_nls_pert, because it is a turning \n"
785  << " point. Consider to use a higher abs_nls_interp_order, for example 4.\n";
786  }
787 
788  out2 << " abs_nls_pert: " << abs_nls_pert[0] << " to "
789  << abs_nls_pert[abs_nls_pert.nelem() - 1]
790  << " (fractional units) in steps of "
791  << abs_nls_pert[1] - abs_nls_pert[0] << " (" << abs_nls_pert.nelem()
792  << " grid points)\n";
793 }
794 
795 /* Workspace method: Doxygen documentation will be auto-generated */
796 void abs_lookupSetup( // WS Output:
797  Vector& abs_p,
798  Vector& abs_t,
799  Vector& abs_t_pert,
800  Matrix& abs_vmrs,
801  ArrayOfArrayOfSpeciesTag& abs_nls,
802  Vector& abs_nls_pert,
803  // WS Input:
804  const Index& atmosphere_dim,
805  const Vector& p_grid,
806  // const Vector& lat_grid,
807  // const Vector& lon_grid,
808  const Tensor3& t_field,
809  const Tensor4& vmr_field,
810  const Index& atmfields_checked,
811  const ArrayOfArrayOfSpeciesTag& abs_species,
812  const Index& abs_p_interp_order,
813  const Index& abs_t_interp_order,
814  const Index& abs_nls_interp_order,
815  // Control Parameters:
816  const Numeric& p_step10,
817  const Numeric& t_step,
818  const Numeric& h2o_step,
819  const Verbosity& verbosity) {
820  // Checks on input parameters:
821 
822  if (atmfields_checked != 1)
823  throw runtime_error(
824  "The atmospheric fields must be flagged to have "
825  "passed a consistency check (atmfields_checked=1).");
826 
827  // We don't actually need lat_grid and lon_grid, but we have them as
828  // input variables, so that we can use the standard functions to
829  // check atmospheric fields and grids. A bit cheesy, but I don't
830  // want to program all the checks explicitly.
831 
832  // Check grids (outcommented the ones that have been done by
833  // atmfields_checkedCalc already):
834  //chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
835 
836  if (p_grid.nelem() < 2) {
837  ostringstream os;
838  os << "You need at least two pressure levels.";
839  throw runtime_error(os.str());
840  }
841 
842  // Check T field:
843  //chk_atm_field("t_field", t_field, atmosphere_dim,
844  // p_grid, lat_grid, lon_grid);
845 
846  // Check VMR field (and abs_species):
847  //chk_atm_field("vmr_field", vmr_field, atmosphere_dim,
848  // abs_species.nelem(), p_grid, lat_grid, lon_grid);
849 
850  // Check the keyword arguments:
851  if (p_step10 <= 0 || t_step <= 0 || h2o_step <= 0) {
852  ostringstream os;
853  os << "The keyword arguments p_step, t_step, and h2o_step must be >0.";
854  throw runtime_error(os.str());
855  }
856 
857  // Ok, all input parameters seem to be reasonable.
858 
859  // For consistency with other code around arts (e.g., correlation
860  // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
861  // convert it here to the natural log:
862  const Numeric p_step = log(pow(10.0, p_step10));
863 
864  // We will need the log of the pressure grid:
865  Vector log_p_grid(p_grid.nelem());
866  transform(log_p_grid, log, p_grid);
867 
868  // const Numeric epsilon = 0.01 * p_step; // This is the epsilon that
869  // // we use for comparing p grid spacings.
870 
871  // Construct abs_p
872  // ---------------
873 
874  ArrayOfNumeric log_abs_p_a; // We take log_abs_p_a as an array of
875  // Numeric, so that we can easily
876  // build it up by appending new elements to the end.
877 
878  // Check whether there are pressure levels that are further apart
879  // (in log(p)) than p_step, and insert additional levels if
880  // necessary:
881 
882  log_abs_p_a.push_back(log_p_grid[0]);
883 
884  for (Index i = 1; i < log_p_grid.nelem(); ++i) {
885  const Numeric dp =
886  log_p_grid[i - 1] - log_p_grid[i]; // The grid is descending.
887 
888  const Numeric dp_by_p_step = dp / p_step;
889  // cout << "dp_by_p_step: " << dp_by_p_step << "\n";
890 
891  // How many times does p_step fit into dp?
892  const Index n = (Index)ceil(dp_by_p_step);
893  // n is the number of intervals that we want to have in the
894  // new grid. The number of additional points to insert is
895  // n-1. But we have to insert the original point as well.
896  // cout << n << "\n";
897 
898  const Numeric ddp = dp / (Numeric)n;
899  // cout << "ddp: " << ddp << "\n";
900 
901  for (Index j = 1; j <= n; ++j)
902  log_abs_p_a.push_back(log_p_grid[i - 1] - (Numeric)j * ddp);
903  }
904 
905  // Copy to a proper vector, we need this also later for
906  // interpolation:
907  Vector log_abs_p(log_abs_p_a.nelem());
908  for (Index i = 0; i < log_abs_p_a.nelem(); ++i) log_abs_p[i] = log_abs_p_a[i];
909 
910  // Copy the new grid to abs_p, removing the log:
911  abs_p.resize(log_abs_p.nelem());
912  transform(abs_p, exp, log_abs_p);
913 
914  // Check that abs_p has enough points for the interpolation order
915  if (abs_p.nelem() < abs_p_interp_order + 1) {
916  ostringstream os;
917  os << "Your pressure grid does not have enough levels for the desired interpolation order.";
918  throw runtime_error(os.str());
919  }
920 
921  // We will also have to interpolate T and VMR profiles to the new
922  // pressure grid. We interpolate in log(p), as usual in ARTS.
923 
924  // Grid positions:
925  ArrayOfGridPos gp(log_abs_p.nelem());
926  gridpos(gp, log_p_grid, log_abs_p);
927 
928  // Interpolation weights:
929  Matrix itw(gp.nelem(), 2);
930  interpweights(itw, gp);
931 
932  // In the 1D case the lookup table is just a lookup table in
933  // pressure. We treat this simple case first.
934  if (1 == atmosphere_dim) {
935  // Reference temperature,
936  // interpolate abs_t from t_field:
937  abs_t.resize(log_abs_p.nelem());
938  interp(abs_t, itw, t_field(joker, 0, 0), gp);
939 
940  // Temperature perturbations:
941  abs_t_pert.resize(0);
942 
943  // Reference VMR profiles,
944  // interpolate abs_vmrs from vmr_field:
945  abs_vmrs.resize(abs_species.nelem(), log_abs_p.nelem());
946  for (Index i = 0; i < abs_species.nelem(); ++i)
947  interp(abs_vmrs(i, joker), itw, vmr_field(i, joker, 0, 0), gp);
948 
949  // Species for which H2O VMR is perturbed:
950  abs_nls.resize(0);
951 
952  // H2O VMR perturbations:
953  abs_nls_pert.resize(0);
954  } else {
955  // 2D or 3D case. We have to set up T and nonlinear species variations.
956 
957  // Make an intelligent choice for the nonlinear species.
958  choose_abs_nls(abs_nls, abs_species, verbosity);
959 
960  // Now comes a part where we analyse the atmospheric fields.
961  // We need to find the max, min, and mean profile for
962  // temperature and VMRs.
963  // We do this on the original p grid, not on the refined p
964  // grid, to be more efficient.
965 
966  // Temperature:
967  Vector tmin(p_grid.nelem());
968  Vector tmax(p_grid.nelem());
969  Vector tmean(p_grid.nelem());
970 
971  for (Index i = 0; i < p_grid.nelem(); ++i) {
972  tmin[i] = min(t_field(i, joker, joker));
973  tmax[i] = max(t_field(i, joker, joker));
974  tmean[i] = mean(t_field(i, joker, joker));
975  }
976 
977  // cout << "Tmin: " << tmin << "\n";
978  // cout << "Tmax: " << tmax << "\n";
979  // cout << "Tmean: " << tmean << "\n";
980 
981  // Calculate mean profiles of all species. (We need all for abs_vmrs
982  // not only H2O.)
983  Matrix vmrmean(abs_species.nelem(), p_grid.nelem());
984  for (Index s = 0; s < abs_species.nelem(); ++s)
985  for (Index i = 0; i < p_grid.nelem(); ++i) {
986  vmrmean(s, i) = mean(vmr_field(s, i, joker, joker));
987  }
988 
989  // If there are NLS, determine H2O statistics:
990 
991  // We have to define these here, outside the if block, because
992  // we need the values later.
993  Vector h2omin(p_grid.nelem());
994  Vector h2omax(p_grid.nelem());
995  const Index h2o_index = find_first_species_tg(
996  abs_species, species_index_from_species_name("H2O"));
997  // We need this inside the if clauses for nonlinear species
998  // treatment. The function returns "-1" if there is no H2O
999  // species. There is a check for that in the next if block, with
1000  // an appropriate runtime error.
1001 
1002  if (0 < abs_nls.nelem()) {
1003  // Check if there really is a H2O species.
1004  if (h2o_index < 0) {
1005  ostringstream os;
1006  os << "Some of your species require nonlinear treatment,\n"
1007  << "but you have no H2O species.";
1008  throw runtime_error(os.str());
1009  }
1010 
1011  for (Index i = 0; i < p_grid.nelem(); ++i) {
1012  h2omin[i] = min(vmr_field(h2o_index, i, joker, joker));
1013  h2omax[i] = max(vmr_field(h2o_index, i, joker, joker));
1014  }
1015 
1016  // cout << "H2Omin: " << h2omin << "\n";
1017  // cout << "H2Omax: " << h2omax << "\n";
1018  // cout << "H2Omean: " << vmrmean(h2o_index,joker) << "\n";
1019  }
1020 
1021  // Interpolate in pressure, set abs_t, abs_vmr...
1022 
1023  // Reference temperature,
1024  // interpolate abs_t from tmean:
1025  abs_t.resize(log_abs_p.nelem());
1026  interp(abs_t, itw, tmean, gp);
1027 
1028  // Temperature perturbations:
1029  choose_abs_t_pert(abs_t_pert,
1030  tmean,
1031  tmin,
1032  tmax,
1033  t_step,
1034  abs_p_interp_order,
1035  abs_t_interp_order,
1036  verbosity);
1037  // cout << "abs_t_pert: " << abs_t_pert << "\n";
1038 
1039  // Reference VMR profiles,
1040  // interpolate abs_vmrs from vmrmean:
1041  abs_vmrs.resize(abs_species.nelem(), log_abs_p.nelem());
1042  for (Index i = 0; i < abs_species.nelem(); ++i)
1043  interp(abs_vmrs(i, joker), itw, vmrmean(i, joker), gp);
1044 
1045  if (0 < abs_nls.nelem()) {
1046  // Construct abs_nls_pert:
1047  choose_abs_nls_pert(abs_nls_pert,
1048  vmrmean(h2o_index, joker),
1049  h2omin,
1050  h2omax,
1051  h2o_step,
1052  abs_p_interp_order,
1053  abs_nls_interp_order,
1054  verbosity);
1055  } else {
1056  // Empty abs_nls_pert:
1057  abs_nls_pert.resize(0);
1058  }
1059  // cout << "abs_nls_pert: " << abs_nls_pert << "\n";
1060  }
1061 }
1062 
1063 /* Workspace method: Doxygen documentation will be auto-generated */
1064 void abs_lookupSetupBatch( // WS Output:
1065  Vector& abs_p,
1066  Vector& abs_t,
1067  Vector& abs_t_pert,
1068  Matrix& abs_vmrs,
1069  ArrayOfArrayOfSpeciesTag& abs_nls,
1070  Vector& abs_nls_pert,
1071  // WS Input:
1072  const ArrayOfArrayOfSpeciesTag& abs_species,
1073  const ArrayOfGriddedField4& batch_fields,
1074  const Index& abs_p_interp_order,
1075  const Index& abs_t_interp_order,
1076  const Index& abs_nls_interp_order,
1077  const Index& atmosphere_dim,
1078  // Control Parameters:
1079  const Numeric& p_step10,
1080  const Numeric& t_step,
1081  const Numeric& h2o_step,
1082  const Vector& extremes,
1083  const Index& robust,
1084  const Index& check_gridnames,
1085  const Verbosity& verbosity) {
1086  CREATE_OUT1;
1087  CREATE_OUT2;
1088  CREATE_OUT3;
1089 
1090  // For consistency with other code around arts (e.g., correlation
1091  // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
1092  // convert it here to the natural log:
1093  const Numeric p_step = log(pow(10.0, p_step10));
1094 
1095  // Derive which abs_species is H2O (required for nonlinear species handling)
1096  // returns -1 if no H2O present
1097  const Index h2o_index = find_first_species_tg(
1098  abs_species, species_index_from_species_name("H2O"));
1099  // cout << "The " << h2o_index+1 << ". species in abs_species is H2O\n";
1100  // cout << "That is, H2O is expected to be the " << indoff+h2o_index
1101  // << ". column of the atmospheric fields\n";
1102 
1103  ArrayOfIndex batch_index(abs_species.nelem());
1104  Index T_index = -1;
1105  Index z_index = -1;
1106 
1107  ArrayOfString species_names(abs_species.nelem());
1108  for (Index i = 0; i < abs_species.nelem(); ++i)
1109  species_names[i] = get_species_name(abs_species[i]);
1110 
1111  const ArrayOfString field_names = batch_fields[0].get_string_grid(0);
1112 
1113  String species_type;
1114  String species_name;
1115  const String delim = "-";
1116 
1117  // Check that the field names in batch_fields are good. (We check
1118  // only the first field in the batch.)
1119  {
1120  ostringstream os;
1121  bool bail = false;
1122 
1123  // One of the abs_species has to be H2O
1124  // (ideally that should be checked later, when we know whether there are
1125  // nonlinear species present. however, currently batch mode REQUIRES H2O to
1126  // be present, so we check that immediately)
1127  if (h2o_index < 0) {
1128  os << "One of the atmospheric fields must be H2O.\n";
1129  bail = true;
1130  }
1131 
1132  // First simply check dimensions of abs_species and field_names
1133  if (field_names.nelem() < 3) {
1134  os << "Atmospheric states in batch_fields must have at\n"
1135  << "least three fields: T, z, and at least one absorption species.";
1136  throw runtime_error(os.str());
1137  }
1138 
1139  if (abs_species.nelem() < 1) {
1140  os << "At least one absorption species needs to be defined "
1141  << "in abs_species.";
1142  throw runtime_error(os.str());
1143  }
1144 
1145  // Check that all required fields are present.
1146  const Index nf = field_names.nelem();
1147  bool found;
1148 
1149  // Looking for temperature field:
1150  found = false;
1151  for (Index fi = 0; fi < nf; ++fi) {
1152  parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1153  if (species_type == "T") {
1154  if (found) {
1155  os << "Only one temperature ('T') field allowed, "
1156  << "but found at least 2.\n";
1157  bail = true;
1158  } else {
1159  found = true;
1160  T_index = fi;
1161  }
1162  }
1163  }
1164  if (!found) {
1165  os << "One temperature ('T') field required, but none found.\n";
1166  bail = true;
1167  }
1168 
1169  // Looking for altitude field:
1170  found = false;
1171  for (Index fi = 0; fi < nf; ++fi) {
1172  parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1173  if (species_type == "z") {
1174  if (found) {
1175  os << "Only one altitude ('z') field allowed, "
1176  << "but found at least 2.\n";
1177  bail = true;
1178  } else {
1179  found = true;
1180  z_index = fi;
1181  }
1182  }
1183  }
1184  if (!found) {
1185  os << "One altitude ('z') field required, but none found.\n";
1186  bail = true;
1187  }
1188 
1189  // Now going over all abs_species elements and match the corresponding
1190  // field_name (we always take the first occurence of a matching field). We
1191  // don't care that batch_fields contains further fields, too. And we can't
1192  // expect the fields being in the same order as the abs_species.
1193  // At the same time, we keep the indices of the fields corresponding to the
1194  // abs_species elements, because later we will need the field data.
1195  Index fi;
1196  for (Index j = 0; j < abs_species.nelem(); ++j) {
1197  found = false;
1198  fi = 0;
1199  while (!found && fi < nf) {
1200  parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1201  // do we have an abs_species type field?
1202  if (species_type == "abs_species") {
1203  parse_atmcompact_speciesname(species_name, field_names[fi], delim);
1204  if (species_name == species_names[j]) {
1205  found = true;
1206  batch_index[j] = fi;
1207  }
1208  }
1209  fi++;
1210  }
1211  if (!found) {
1212  os << "No field for absorption species '" << species_names[j]
1213  << "'found.\n";
1214  bail = true;
1215  }
1216  }
1217 
1218  os << "Your field names are:\n" << field_names;
1219 
1220  if (bail) throw runtime_error(os.str());
1221  }
1222 
1223  // FIXME: Adjustment of min/max values for Jacobian perturbations is still missing.
1224 
1225  // Make an intelligent choice for the nonlinear species.
1226  choose_abs_nls(abs_nls, abs_species, verbosity);
1227 
1228  // Find out maximum and minimum pressure and check that pressure grid is decreasing.
1229  Numeric maxp = batch_fields[0].get_numeric_grid(GFIELD4_P_GRID)[0];
1230  Numeric minp = batch_fields[0].get_numeric_grid(
1231  GFIELD4_P_GRID)[batch_fields[0].get_numeric_grid(GFIELD4_P_GRID).nelem() -
1232  1];
1233 
1234  ArrayOfIndex valid_field_indices;
1235  for (Index i = 0; i < batch_fields.nelem(); ++i) {
1236  // Local variables for atmfields_check.
1237  Index atmfields_checked;
1238  Tensor4 t4_dummy;
1239  Tensor3 t3_dummy;
1240 
1241  Vector p_grid;
1242  Vector lat_grid;
1243  Vector lon_grid;
1244  Tensor3 t_field;
1245  Tensor3 z_field;
1246  Tensor4 vmr_field;
1247  Tensor4 particle_bulkprop_field;
1248  ArrayOfString particle_bulkprop_names;
1249  GriddedField4 atm_fields_compact;
1250  SpeciesAuxData partition_functions;
1251  Index abs_f_interp_order;
1252 
1253  // Extract fields from atmfield and check their validity.
1254  // This closes the loophole when only calculating lookup tables.
1255  atm_fields_compact = batch_fields[i];
1256 
1258  lat_grid,
1259  lon_grid,
1260  t_field,
1261  z_field,
1262  vmr_field,
1263  particle_bulkprop_field,
1264  particle_bulkprop_names,
1265  abs_species,
1266  atm_fields_compact,
1267  atmosphere_dim,
1268  "-",
1269  0,
1270  check_gridnames,
1271  verbosity);
1272 
1273  try {
1274  atmfields_checkedCalc(atmfields_checked,
1275  atmosphere_dim,
1276  p_grid,
1277  lat_grid,
1278  lon_grid,
1279  abs_species,
1280  t_field,
1281  vmr_field,
1282  t3_dummy,
1283  t3_dummy,
1284  t3_dummy,
1285  t3_dummy,
1286  t3_dummy,
1287  t3_dummy,
1288  partition_functions,
1289  abs_f_interp_order,
1290  0,
1291  0,
1292  verbosity);
1293  } catch (const std::exception& e) {
1294  // If `robust`, skip field and continue, ...
1295  if (robust) {
1296  out1 << " WARNING! Skipped invalid atmfield "
1297  << "at batch_atmfield index " << i << ".\n"
1298  << "The runtime error produced was:\n"
1299  << e.what() << "\n";
1300  continue;
1301  }
1302  // ... else throw an error.
1303  else {
1304  stringstream err;
1305  err << "Invalid atmfield at batch_atmfield index " << i << ".\n"
1306  << "The runtime error produced was:\n"
1307  << e.what() << "\n";
1308  throw std::runtime_error(err.str());
1309  }
1310  };
1311  valid_field_indices.push_back(i); // Append index to list of valid fields.
1312 
1313  if (maxp < batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)[0])
1314  maxp = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)[0];
1315  if (minp >
1316  batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)
1317  [batch_fields[i].get_numeric_grid(GFIELD4_P_GRID).nelem() - 1])
1318  minp = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)
1319  [batch_fields[i].get_numeric_grid(GFIELD4_P_GRID).nelem() - 1];
1320  }
1321  // cout << " minp/maxp: " << minp << " / " << maxp << "\n";
1322 
1323  // Information on the number of skipped atmospheres.
1324  if (batch_fields.nelem() > valid_field_indices.nelem()) {
1325  out1 << " " << batch_fields.nelem() - valid_field_indices.nelem()
1326  << " out of " << batch_fields.nelem() << " atmospheres ignored.\n";
1327  }
1328 
1329  // Throw error if no atmfield passed the check.
1330  if (valid_field_indices.nelem() < 1) {
1331  stringstream err;
1332  err << "You need at least one valid profile.\n"
1333  << "It seems that no atmfield passed the checks!\n";
1334  throw std::runtime_error(err.str());
1335  }
1336 
1337  if (maxp == minp) {
1338  ostringstream os;
1339  os << "You need at least two pressure levels.";
1340  throw runtime_error(os.str());
1341  }
1342 
1343  // We construct the pressure grid as follows:
1344  // - Everything is done in log(p).
1345  // - Start with maxp and go down in steps of p_step until we are <= minp.
1346  // - Adjust the final pressure value to be exactly minp, otherwise
1347  // we have problems in getting min, max, and mean values for this
1348  // point later.
1349  Index np = (Index)ceil((log(maxp) - log(minp)) / p_step) + 1;
1350  // The +1 above has to be there because we must include also both
1351  // grid end points.
1352 
1353  // If np is too small for the interpolation order, we increase it:
1354  if (np < abs_p_interp_order + 1) np = abs_p_interp_order + 1;
1355 
1356  Vector log_abs_p(log(maxp), np, -p_step);
1357  log_abs_p[np - 1] = log(minp);
1358 
1359  abs_p.resize(np);
1360  transform(abs_p, exp, log_abs_p);
1361  out2 << " abs_p: " << abs_p[0] << " Pa to " << abs_p[np - 1]
1362  << " Pa in log10 steps of " << p_step10 << " (" << np
1363  << " grid points)\n";
1364 
1365  // Now we have to determine the statistics of T and VMRs, we need
1366  // profiles of min, max, and mean of these, on the abs_p grid.
1367 
1368  // Index n_variables = batch_fields[0].data.nbooks();
1369  // we will do statistics for data fields excluding the scat_species fields
1370  Index n_variables = 2 + abs_species.nelem();
1371 
1372  // The first dimension of datamin, datamax, and datamean is the
1373  // variable (T,Z,H2O,O3,...). The second dimension is pressure. We
1374  // assume all elements of the batch have the same variables.
1375 
1376  Matrix datamin(n_variables, np, numeric_limits<Numeric>::max());
1377  Matrix datamax(n_variables, np, numeric_limits<Numeric>::min());
1378  // The limits here are from the header file <limits>
1379  Matrix datamean(n_variables, np, 0);
1380  Vector mean_norm(np, 0); // Divide by this to get mean.
1381 
1382  // We will loop over all batch cases to analyze the statistics and
1383  // calculate reference profiles. As a little side project, we will
1384  // also calculate the absolute min and max of temperature and
1385  // humidity. This is handy as input to abs_lookupSetupWide.
1386  Numeric mint = +1e99;
1387  Numeric maxt = -1e99;
1388  Numeric minh2o = +1e99;
1389  Numeric maxh2o = -1e99;
1390 
1391  // Loop over valid atmfields.
1392  for (Index vi = 0; vi < valid_field_indices.nelem(); ++vi) {
1393  // Get internal field index.
1394  Index i = valid_field_indices[vi];
1395 
1396  // Check that really each case has the same variables (see
1397  // comment above.)
1398  if (batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES) !=
1399  batch_fields[0].get_string_grid(GFIELD4_FIELD_NAMES))
1400  throw runtime_error(
1401  "All batch atmospheres must contain the same field names.");
1402 
1403  for (Index j = 0;
1404  j < batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES).nelem();
1405  ++j)
1406  if (batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES)[j] !=
1407  batch_fields[0].get_string_grid(GFIELD4_FIELD_NAMES)[j])
1408  throw runtime_error(
1409  "All batch atmospheres must contain the same individual field names.");
1410 
1411  // Create convenient handles:
1412  const Vector& p_grid = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID);
1413  const Tensor4& data = batch_fields[i].data;
1414 
1415  // Update our global max/min values for T and H2O:
1416 
1417  // We have to loop over pressure, latitudes, and longitudes
1418  // here. The dimensions of data are:
1419  // data[N_fields][N_p][N_lat][N_lon]
1420 
1421  for (Index ip = 0; ip < data.npages(); ++ip)
1422  for (Index ilat = 0; ilat < data.nrows(); ++ilat)
1423  for (Index ilon = 0; ilon < data.ncols(); ++ilon) {
1424  // Field T_index is temperature:
1425  if (data(T_index, ip, ilat, ilon) < mint)
1426  mint = data(T_index, ip, ilat, ilon);
1427  if (data(T_index, ip, ilat, ilon) > maxt)
1428  maxt = data(T_index, ip, ilat, ilon);
1429  // Field batch_index[h2o_index] is H2O:
1430  if (data(batch_index[h2o_index], ip, ilat, ilon) < minh2o) {
1431  minh2o = data(batch_index[h2o_index], ip, ilat, ilon);
1432  }
1433  if (data(batch_index[h2o_index], ip, ilat, ilon) > maxh2o) {
1434  maxh2o = data(batch_index[h2o_index], ip, ilat, ilon);
1435  }
1436  }
1437 
1438  // Interpolate the current batch fields to the abs_p grid. We
1439  // have to do this for each batch case, since the grids could
1440  // all be different.
1441 
1442  Vector log_p_grid(p_grid.nelem());
1443  transform(log_p_grid, log, p_grid);
1444 
1445  // There is a catch here: We can only use the part of abs_p that
1446  // is inside the current pressure grid p_grid, otherwise we
1447  // would have to extrapolate.
1448  // The eps_bottom and eps_top account for the little bit of
1449  // extrapolation that is allowed.
1450 
1451  const Numeric eps_bottom = (log_p_grid[0] - log_p_grid[1]) / 2.1;
1452  Index first_p = 0;
1453  while (log_abs_p[first_p] > log_p_grid[0] + eps_bottom) ++first_p;
1454 
1455  const Numeric eps_top = (log_p_grid[log_p_grid.nelem() - 2] -
1456  log_p_grid[log_p_grid.nelem() - 1]) /
1457  2.1;
1458  Index last_p = log_abs_p.nelem() - 1;
1459  while (log_abs_p[last_p] < log_p_grid[log_p_grid.nelem() - 1] - eps_top)
1460  --last_p;
1461 
1462  const Index extent_p = last_p - first_p + 1;
1463 
1464  // This was too complicated to get right:
1465  // const Index first_p = (Index) round ( (log_abs_p[0] - log_p_grid[0]) / p_step);
1466  // const Index extent_p = (Index) round ( (log_abs_p[first_p] - log_p_grid[log_p_grid.nelem()-1]) / p_step) + 1;
1467 
1468  ConstVectorView active_log_abs_p = log_abs_p[Range(first_p, extent_p)];
1469 
1470  // cout << "first_p / last_p / extent_p : " << first_p << " / " << last_p << " / " << extent_p << "\n";
1471  // cout << "log_p_grid: " << log_p_grid << "\n";
1472  // cout << "log_abs_p: " << log_abs_p << "\n";
1473  // cout << "active_log_abs_p: " << active_log_abs_p << "\n";
1474  // cout << "=============================================================\n";
1475  // arts_exit();
1476 
1477  // Grid positions:
1478  ArrayOfGridPos gp(active_log_abs_p.nelem());
1479  gridpos(gp, log_p_grid, active_log_abs_p);
1480  // gridpos(gp, log_p_grid, active_log_abs_p, 100);
1481  // We allow much more extrapolation here than normal (0.5 is
1482  // normal). If we do not do this, then we get problems for
1483  // p_grids that are much finer than abs_p.
1484 
1485  // Interpolation weights:
1486  Matrix itw(gp.nelem(), 2);
1487  interpweights(itw, gp);
1488 
1489  // We have to loop over fields, latitudes, and longitudes here, doing the
1490  // pressure interpolation for all. The dimensions of data are:
1491  // data[N_fields][N_p][N_lat][N_lon]
1492  // For data_interp we reduce data by the particle related fields, i.e.,
1493  // the ones in between the first two and the last N=abs_species.nelem().
1494  // Tensor4 data_interp(data.nbooks(),
1495  Tensor4 data_interp(
1496  n_variables, active_log_abs_p.nelem(), data.nrows(), data.ncols());
1497 
1498  for (Index lo = 0; lo < data.ncols(); ++lo)
1499  for (Index la = 0; la < data.nrows(); ++la) {
1500  // for (Index fi=0; fi<data.nbooks(); ++fi)
1501  // we have to handle T/z and abs_species parts separately since they
1502  // can have scat_species in between, which we want to ignore
1503  interp(data_interp(0, joker, la, lo),
1504  itw,
1505  data(T_index, joker, la, lo),
1506  gp);
1507  interp(data_interp(1, joker, la, lo),
1508  itw,
1509  data(z_index, joker, la, lo),
1510  gp);
1511  for (Index fi = 0; fi < abs_species.nelem(); ++fi)
1512  interp(data_interp(fi + 2, joker, la, lo),
1513  itw,
1514  data(batch_index[fi], joker, la, lo),
1515  gp);
1516  }
1517 
1518  // Now update our datamin, datamax, and datamean variables.
1519  // We need the min and max only for the T and H2O profile,
1520  // not for others. But we need the mean for all. We are just
1521  // hopping over the case that we do not need below. This is not
1522  // very clean, but efficient. And it avoids handling all the
1523  // different cases explicitly.
1524  for (Index lo = 0; lo < data_interp.ncols(); ++lo)
1525  for (Index la = 0; la < data_interp.nrows(); ++la) {
1526  for (Index fi = 0; fi < data_interp.nbooks(); ++fi) {
1527  if (1 != fi) // We skip the z field, which we do not need
1528  for (Index pr = 0; pr < data_interp.npages(); ++pr) {
1529  // Min and max only needed for T and H2o
1530  if (0 == fi || (h2o_index + 2) == fi) {
1531  if (data_interp(fi, pr, la, lo) < datamin(fi, first_p + pr))
1532  datamin(fi, first_p + pr) = data_interp(fi, pr, la, lo);
1533  if (data_interp(fi, pr, la, lo) > datamax(fi, first_p + pr))
1534  datamax(fi, first_p + pr) = data_interp(fi, pr, la, lo);
1535  }
1536 
1537  datamean(fi, first_p + pr) += data_interp(fi, pr, la, lo);
1538  }
1539  }
1540 
1541  // The mean_norm is actually a bit tricky. It depends on
1542  // pressure, since different numbers of cases contribute
1543  // to the mean for different pressures. At the very eges
1544  // of the grid, typically only a single case contributes.
1545 
1546  mean_norm[Range(first_p, extent_p)] += 1;
1547  }
1548  }
1549 
1550  out2 << " Global statistics:\n"
1551  << " min(p) / max(p) [Pa]: " << minp << " / " << maxp << "\n"
1552  << " min(T) / max(T) [K]: " << mint << " / " << maxt << "\n"
1553  << " min(H2O) / max(H2O) [VMR]: " << minh2o << " / " << maxh2o << "\n";
1554 
1555  // Divide mean by mean_norm to get the mean:
1556  assert(np == mean_norm.nelem());
1557  for (Index fi = 0; fi < datamean.nrows(); ++fi)
1558  if (1 != fi) // We skip the z field, which we do not need
1559  for (Index pi = 0; pi < np; ++pi) {
1560  // We do this in an explicit loop here, since we have to
1561  // check whether there really were data points to calculate
1562  // the mean at each level.
1563  if (0 < mean_norm[pi])
1564  datamean(fi, pi) /= mean_norm[pi];
1565  else {
1566  ostringstream os;
1567  os << "No data at pressure level " << pi + 1 << " of " << np << " ("
1568  << abs_p[pi] << " hPa).";
1569  throw runtime_error(os.str());
1570  }
1571  }
1572 
1573  // If we do higher order pressure interpolation, then we should
1574  // smooth the reference profiles with a boxcar filter of width
1575  // p_interp_order+1. Otherwise we get numerical problems if there
1576  // are any sharp spikes in the reference profiles.
1577  assert(log_abs_p.nelem() == np);
1578  Matrix smooth_datamean(datamean.nrows(), datamean.ncols(), 0);
1579  for (Index i = 0; i < np; ++i) {
1580  GridPosPoly gp;
1581  gridpos_poly(gp, log_abs_p, log_abs_p[i], abs_p_interp_order);
1582 
1583  // We do this in practice by using the indices returned by
1584  // gridpos_poly. We simply take a mean over all points that
1585  // would be used in the interpolation.
1586 
1587  for (Index fi = 0; fi < datamean.nrows(); ++fi)
1588  if (1 != fi) // We skip the z field, which we do not need
1589  {
1590  for (Index j = 0; j < gp.idx.nelem(); ++j) {
1591  smooth_datamean(fi, i) += datamean(fi, gp.idx[j]);
1592  }
1593  smooth_datamean(fi, i) /= (Numeric)gp.idx.nelem();
1594  }
1595  // cout << "H2O-raw / H2O-smooth: "
1596  // << datamean(h2o_index+2,i) << " / "
1597  // << smooth_datamean(h2o_index+2,i) << "\n";
1598  }
1599 
1600  // There is another complication: If the (smoothed) mean for the H2O
1601  // reference profile is 0, then we have to adjust both mean and max
1602  // value to a non-zero number, otherwise the reference profile will
1603  // be zero, and we will get numerical problems later on when we
1604  // divide by the reference value. So, we set it here to 1e-9.
1605  for (Index i = 0; i < np; ++i) {
1606  // Assert that really H2O has index batch_index[h2o_index] VMR field list
1607  // and h2o_index in abs_species
1609  species_type, field_names[batch_index[h2o_index]], delim);
1610  assert(species_type == "abs_species");
1612  species_name, field_names[batch_index[h2o_index]], delim);
1613  assert("H2O" == species_name);
1614  assert("H2O" == species_names[h2o_index]);
1615 
1616  // Find mean and max H2O for this level:
1617  Numeric& mean_h2o = smooth_datamean(h2o_index + 2, i);
1618  Numeric& max_h2o = datamax(h2o_index + 2, i);
1619  if (mean_h2o <= 0) {
1620  mean_h2o = 1e-9;
1621  max_h2o = 1e-9;
1622  out3 << " H2O profile contained zero values, adjusted to 1e-9.\n";
1623  }
1624  }
1625 
1626  // Set abs_t:
1627  abs_t.resize(np);
1628  abs_t = smooth_datamean(0, joker);
1629  // cout << "abs_t: " << abs_t << "\n\n";
1630  // cout << "tmin: " << datamin(0,joker) << "\n\n";
1631  // cout << "tmax: " << datamax(0,joker) << "\n";
1632 
1633  // Set abs_vmrs:
1634  assert(abs_species.nelem() == smooth_datamean.nrows() - 2);
1635  abs_vmrs.resize(abs_species.nelem(), np);
1636  abs_vmrs = smooth_datamean(Range(2, abs_species.nelem()), joker);
1637  // cout << "\n\nabs_vmrs: " << abs_vmrs << "\n\n";
1638 
1639  // Construct abs_t_pert:
1640  ConstVectorView tmin = datamin(0, joker);
1641  ConstVectorView tmax = datamax(0, joker);
1642  choose_abs_t_pert(abs_t_pert,
1643  abs_t,
1644  tmin,
1645  tmax,
1646  t_step,
1647  abs_p_interp_order,
1648  abs_t_interp_order,
1649  verbosity);
1650  // cout << "abs_t_pert: " << abs_t_pert << "\n";
1651 
1652  // Construct abs_nls_pert:
1653  ConstVectorView h2omin = datamin(h2o_index + 2, joker);
1654  ConstVectorView h2omax = datamax(h2o_index + 2, joker);
1655  choose_abs_nls_pert(abs_nls_pert,
1656  abs_vmrs(h2o_index, joker),
1657  h2omin,
1658  h2omax,
1659  h2o_step,
1660  abs_p_interp_order,
1661  abs_nls_interp_order,
1662  verbosity);
1663  // cout << "abs_nls_pert: " << abs_nls_pert << "\n";
1664 
1665  // Append the explicitly given user extreme values, if necessary:
1666  if (0 != extremes.nelem()) {
1667  // There must be 4 values in this case: t_min, t_max, h2o_min, h2o_max
1668  if (4 != extremes.nelem()) {
1669  ostringstream os;
1670  os << "There must be exactly 4 elements in extremes:\n"
1671  << "min(abs_t_pert), max(abs_t_pert), min(abs_nls_pert), max(abs_nls_pert)";
1672  throw runtime_error(os.str());
1673  }
1674 
1675  // t_min:
1676  if (extremes[0] < abs_t_pert[0]) {
1677  Vector dummy = abs_t_pert;
1678  abs_t_pert.resize(abs_t_pert.nelem() + 1);
1679  abs_t_pert[0] = extremes[0];
1680  abs_t_pert[Range(1, dummy.nelem())] = dummy;
1681  out2 << " Added min extreme value for abs_t_pert: " << abs_t_pert[0]
1682  << "\n";
1683  }
1684 
1685  // t_max:
1686  if (extremes[1] > abs_t_pert[abs_t_pert.nelem() - 1]) {
1687  Vector dummy = abs_t_pert;
1688  abs_t_pert.resize(abs_t_pert.nelem() + 1);
1689  abs_t_pert[Range(0, dummy.nelem())] = dummy;
1690  abs_t_pert[abs_t_pert.nelem() - 1] = extremes[1];
1691  out2 << " Added max extreme value for abs_t_pert: "
1692  << abs_t_pert[abs_t_pert.nelem() - 1] << "\n";
1693  }
1694 
1695  // nls_min:
1696  if (extremes[2] < abs_nls_pert[0]) {
1697  Vector dummy = abs_nls_pert;
1698  abs_nls_pert.resize(abs_nls_pert.nelem() + 1);
1699  abs_nls_pert[0] = extremes[2];
1700  abs_nls_pert[Range(1, dummy.nelem())] = dummy;
1701  out2 << " Added min extreme value for abs_nls_pert: " << abs_nls_pert[0]
1702  << "\n";
1703  }
1704 
1705  // nls_max:
1706  if (extremes[3] > abs_nls_pert[abs_nls_pert.nelem() - 1]) {
1707  Vector dummy = abs_nls_pert;
1708  abs_nls_pert.resize(abs_nls_pert.nelem() + 1);
1709  abs_nls_pert[Range(0, dummy.nelem())] = dummy;
1710  abs_nls_pert[abs_nls_pert.nelem() - 1] = extremes[3];
1711  out2 << " Added max extreme value for abs_nls_pert: "
1712  << abs_nls_pert[abs_nls_pert.nelem() - 1] << "\n";
1713  }
1714  }
1715 }
1716 
1717 void abs_lookupSetupWide( // WS Output:
1718  Vector& abs_p,
1719  Vector& abs_t,
1720  Vector& abs_t_pert,
1721  Matrix& abs_vmrs,
1722  ArrayOfArrayOfSpeciesTag& abs_nls,
1723  Vector& abs_nls_pert,
1724  // WS Input:
1725  const ArrayOfArrayOfSpeciesTag& abs_species,
1726  const Index& abs_p_interp_order,
1727  const Index& abs_t_interp_order,
1728  const Index& abs_nls_interp_order,
1729  // Control Parameters:
1730  const Numeric& p_min,
1731  const Numeric& p_max,
1732  const Numeric& p_step10,
1733  const Numeric& t_min,
1734  const Numeric& t_max,
1735  const Numeric& h2o_min,
1736  const Numeric& h2o_max,
1737  const Verbosity& verbosity) {
1738  CREATE_OUT2;
1739 
1740  // For consistency with other code around arts (e.g., correlation
1741  // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
1742  // convert it here to the natural log:
1743  const Numeric p_step = log(pow(10.0, p_step10));
1744 
1745  // Make an intelligent choice for the nonlinear species.
1746  choose_abs_nls(abs_nls, abs_species, verbosity);
1747 
1748  // 1. Fix pressure grid abs_p
1749  // --------------------------
1750 
1751  // We construct the pressure grid as follows:
1752  // - Everything is done in log(p).
1753  // - Start with p_max and go down in steps of p_step until we are <= p_min.
1754 
1755  Index np = (Index)ceil((log(p_max) - log(p_min)) / p_step) + 1;
1756  // The +1 above has to be there because we must include also both
1757  // grid end points.
1758 
1759  // If np is too small for the interpolation order, we increase it:
1760  if (np < abs_p_interp_order + 1) np = abs_p_interp_order + 1;
1761 
1762  Vector log_abs_p(log(p_max), np, -p_step);
1763 
1764  abs_p.resize(np);
1765  transform(abs_p, exp, log_abs_p);
1766  out2 << " abs_p: " << abs_p[0] << " Pa to " << abs_p[np - 1]
1767  << " Pa in log10 steps of " << p_step10 << " (" << np
1768  << " grid points)\n";
1769 
1770  // 2. Fix reference temperature profile abs_t and temperature perturbations
1771  // ------------------------------------------------------------------------
1772 
1773  // We simply take a constant reference profile.
1774 
1775  Numeric t_ref = (t_min + t_max) / 2;
1776 
1777  abs_t.resize(np);
1778  abs_t = t_ref; // Assign same value to all elements.
1779 
1780  // We have to make vectors out of t_min and t_max, so we can use
1781  // them in the choose_abs_t_pert function call:
1782  Vector min_prof(np), max_prof(np);
1783  min_prof = t_min;
1784  max_prof = t_max;
1785 
1786  // Chose temperature perturbations:
1787  choose_abs_t_pert(abs_t_pert,
1788  abs_t,
1789  min_prof,
1790  max_prof,
1791  20,
1792  abs_p_interp_order,
1793  abs_t_interp_order,
1794  verbosity);
1795 
1796  // 3. Fix reference H2O profile and abs_nls_pert
1797  // ---------------------------------------------
1798 
1799  // We take a constant reference profile of 1000ppm (=1e-3) for H2O
1800  Numeric const h2o_ref = 1e-3;
1801 
1802  // And 1 ppt (1e-9) as default for all VMRs
1803  Numeric const other_ref = 1e-9;
1804 
1805  // We have to assign this value to all pressures of the H2O profile,
1806  // and 0 to all other profiles.
1807 
1808  // abs_vmrs has dimension [n_species, np].
1809  abs_vmrs.resize(abs_species.nelem(), np);
1810  abs_vmrs = other_ref;
1811 
1812  // We look for O2 and N2, and assign constant values to them.
1813  // The values are from Wallace&Hobbs, 2nd edition.
1814 
1815  const Index o2_index =
1817  if (o2_index >= 0) {
1818  abs_vmrs(o2_index, joker) = 0.2095;
1819  }
1820 
1821  const Index n2_index =
1823  if (n2_index >= 0) {
1824  abs_vmrs(n2_index, joker) = 0.7808;
1825  }
1826 
1827  // Which species is H2O?
1828  const Index h2o_index = find_first_species_tg(
1829  abs_species, species_index_from_species_name("H2O"));
1830 
1831  // The function returns "-1" if there is no H2O
1832  // species.
1833  if (0 < abs_nls.nelem()) {
1834  if (h2o_index < 0) {
1835  ostringstream os;
1836  os << "Some of your species require nonlinear treatment,\n"
1837  << "but you have no H2O species.";
1838  throw runtime_error(os.str());
1839  }
1840 
1841  // Assign constant reference value to all H2O levels:
1842  abs_vmrs(h2o_index, joker) = h2o_ref;
1843 
1844  // We have to make vectors out of h2o_min and h2o_max, so we can use
1845  // them in the choose_abs_nls_pert function call.
1846  // We re-use the vectors min_prof and max_prof that we have
1847  // defined above.
1848  min_prof = h2o_min;
1849  max_prof = h2o_max;
1850 
1851  // Construct abs_nls_pert:
1852  choose_abs_nls_pert(abs_nls_pert,
1853  abs_vmrs(h2o_index, joker),
1854  min_prof,
1855  max_prof,
1856  1e99,
1857  abs_p_interp_order,
1858  abs_nls_interp_order,
1859  verbosity);
1860  } else {
1861  CREATE_OUT1;
1862  out1 << " WARNING:\n"
1863  << " You have no species that require H2O variations.\n"
1864  << " This case might work, but it has never been tested.\n"
1865  << " Please test it, then remove this warning.\n";
1866  }
1867 }
1868 
1869 /* Workspace method: Doxygen documentation will be auto-generated */
1870 void abs_speciesAdd( // WS Output:
1871  ArrayOfArrayOfSpeciesTag& abs_species,
1872  Index& propmat_clearsky_agenda_checked,
1873  Index& abs_xsec_agenda_checked,
1874  // Control Parameters:
1875  const ArrayOfString& names,
1876  const Verbosity& verbosity) {
1877  CREATE_OUT3;
1878 
1879  // Invalidate agenda check flags
1880  propmat_clearsky_agenda_checked = false;
1881  abs_xsec_agenda_checked = false;
1882 
1883  // Size of initial array
1884  Index n_gs = abs_species.nelem();
1885 
1886  // Temporary ArrayOfSpeciesTag
1888 
1889  // Each element of the array of Strings names defines one tag
1890  // group. Let's work through them one by one.
1891  for (Index i = 0; i < names.nelem(); ++i) {
1892  array_species_tag_from_string(temp, names[i]);
1893  abs_species.push_back(temp);
1894  }
1895 
1896  check_abs_species(abs_species);
1897 
1898  // Print list of tag groups to the most verbose output stream:
1899  out3 << " Added tag groups:";
1900  for (Index i = n_gs; i < abs_species.nelem(); ++i) {
1901  out3 << "\n " << i << ":";
1902  for (Index s = 0; s < abs_species[i].nelem(); ++s) {
1903  out3 << " " << abs_species[i][s].Name();
1904  }
1905  }
1906  out3 << '\n';
1907 }
1908 
1909 /* Workspace method: Doxygen documentation will be auto-generated */
1910 void abs_speciesAdd2( // WS Output:
1911  Workspace& ws,
1912  ArrayOfArrayOfSpeciesTag& abs_species,
1914  Agenda& jacobian_agenda,
1915  Index& propmat_clearsky_agenda_checked,
1916  Index& abs_xsec_agenda_checked,
1917  // WS Input:
1918  const Index& atmosphere_dim,
1919  const Vector& p_grid,
1920  const Vector& lat_grid,
1921  const Vector& lon_grid,
1922  // WS Generic Input:
1923  const Vector& rq_p_grid,
1924  const Vector& rq_lat_grid,
1925  const Vector& rq_lon_grid,
1926  // Control Parameters:
1927  const String& species,
1928  const String& mode,
1929  const Verbosity& verbosity) {
1930  CREATE_OUT3;
1931 
1932  // Invalidate agenda check flags
1933  propmat_clearsky_agenda_checked = false;
1934  abs_xsec_agenda_checked = false;
1935 
1936  // Add species to *abs_species*
1937  ArrayOfSpeciesTag tags;
1938  array_species_tag_from_string(tags, species);
1939  abs_species.push_back(tags);
1940 
1941  check_abs_species(abs_species);
1942 
1943  // Print list of added tag group to the most verbose output stream:
1944  out3 << " Appended tag group:";
1945  out3 << "\n " << abs_species.nelem() - 1 << ":";
1946  for (Index s = 0; s < tags.nelem(); ++s) {
1947  out3 << " " << tags[s].Name();
1948  }
1949  out3 << '\n';
1950 
1951  // Do retrieval part
1953  jq,
1954  jacobian_agenda,
1955  atmosphere_dim,
1956  p_grid,
1957  lat_grid,
1958  lon_grid,
1959  rq_p_grid,
1960  rq_lat_grid,
1961  rq_lon_grid,
1962  species,
1963  mode,
1964  1,
1965  verbosity);
1966 }
1967 
1968 /* Workspace method: Doxygen documentation will be auto-generated */
1970  abs_species.resize(0);
1971 }
1972 
1973 /* Workspace method: Doxygen documentation will be auto-generated */
1974 void abs_speciesSet( // WS Output:
1975  ArrayOfArrayOfSpeciesTag& abs_species,
1976  Index& abs_xsec_agenda_checked,
1977  Index& propmat_clearsky_agenda_checked,
1978  // Control Parameters:
1979  const ArrayOfString& names,
1980  const Verbosity& verbosity) {
1981  CREATE_OUT3;
1982 
1983  // Invalidate agenda check flags
1984  propmat_clearsky_agenda_checked = false;
1985  abs_xsec_agenda_checked = false;
1986 
1987  abs_species.resize(names.nelem());
1988 
1989  //cout << "Names: " << names << "\n";
1990 
1991  // Each element of the array of Strings names defines one tag
1992  // group. Let's work through them one by one.
1993  for (Index i = 0; i < names.nelem(); ++i) {
1994  // This part has now been moved to array_species_tag_from_string.
1995  // Call this function.
1996  array_species_tag_from_string(abs_species[i], names[i]);
1997  }
1998 
1999  check_abs_species(abs_species);
2000 
2001  // Print list of tag groups to the most verbose output stream:
2002  out3 << " Defined tag groups: ";
2003  for (Index i = 0; i < abs_species.nelem(); ++i) {
2004  out3 << "\n " << i << ":";
2005  for (Index s = 0; s < abs_species[i].nelem(); ++s)
2006  out3 << " " << abs_species[i][s].Name();
2007  }
2008  out3 << '\n';
2009 }
2010 
2011 /* Workspace method: Doxygen documentation will be auto-generated */
2013  Index& abs_lookup_is_adapted,
2014  const ArrayOfArrayOfSpeciesTag& abs_species,
2015  const Vector& f_grid,
2016  const Verbosity& verbosity) {
2017  abs_lookup.Adapt(abs_species, f_grid, verbosity);
2018  abs_lookup_is_adapted = 1;
2019 }
2020 
2021 /* Workspace method: Doxygen documentation will be auto-generated */
2023  ArrayOfPropagationMatrix& propmat_clearsky,
2024  ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
2025  const GasAbsLookup& abs_lookup,
2026  const Index& abs_lookup_is_adapted,
2027  const Index& abs_p_interp_order,
2028  const Index& abs_t_interp_order,
2029  const Index& abs_nls_interp_order,
2030  const Index& abs_f_interp_order,
2031  const Vector& f_grid,
2032  const Numeric& a_pressure,
2033  const Numeric& a_temperature,
2034  const Vector& a_vmr_list,
2035  const ArrayOfRetrievalQuantity& jacobian_quantities,
2036  const Numeric& extpolfac,
2037  const Verbosity& verbosity) {
2038  CREATE_OUT3;
2039 
2040  // Variables needed by abs_lookup.Extract:
2041  Matrix abs_scalar_gas, dabs_scalar_gas_df, dabs_scalar_gas_dt;
2042 
2043  // Check if the table has been adapted:
2044  if (1 != abs_lookup_is_adapted)
2045  throw runtime_error(
2046  "Gas absorption lookup table must be adapted,\n"
2047  "use method abs_lookupAdapt.");
2048 
2049  const bool do_jac = supports_lookup(jacobian_quantities);
2050  const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
2051  const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
2052  const Numeric df = frequency_perturbation(jacobian_quantities);
2053  const Numeric dt = temperature_perturbation(jacobian_quantities);
2054  const ArrayOfIndex jacobian_quantities_position =
2055  equivalent_propmattype_indexes(jacobian_quantities);
2056 
2057  // The combination of doing frequency jacobian together with an
2058  // absorption lookup table is quite dangerous. If the frequency
2059  // interpolation order for the table is zero, the Jacobian will be
2060  // zero, and the cause for this may be difficult for a user to
2061  // find. So we do not allow this combination.
2062  if (do_freq_jac and (1 > abs_f_interp_order))
2063  throw std::runtime_error("Wind/frequency Jacobian is not possible without at least first\n"
2064  "order frequency interpolation in the lookup table. Please use\n"
2065  "abs_f_interp_order>0 or remove wind/frequency Jacobian.");
2066 
2067  // The function we are going to call here is one of the few helper
2068  // functions that adjust the size of their output argument
2069  // automatically.
2070  abs_lookup.Extract(abs_scalar_gas,
2071  abs_p_interp_order,
2072  abs_t_interp_order,
2073  abs_nls_interp_order,
2074  abs_f_interp_order,
2075  a_pressure,
2076  a_temperature,
2077  a_vmr_list,
2078  f_grid,
2079  extpolfac);
2080  if (do_freq_jac) {
2081  Vector dfreq = f_grid;
2082  dfreq += df;
2083  abs_lookup.Extract(dabs_scalar_gas_df,
2084  abs_p_interp_order,
2085  abs_t_interp_order,
2086  abs_nls_interp_order,
2087  abs_f_interp_order,
2088  a_pressure,
2089  a_temperature,
2090  a_vmr_list,
2091  dfreq,
2092  extpolfac);
2093  }
2094  if (do_temp_jac) {
2095  const Numeric dtemp = a_temperature + dt;
2096  abs_lookup.Extract(dabs_scalar_gas_dt,
2097  abs_p_interp_order,
2098  abs_t_interp_order,
2099  abs_nls_interp_order,
2100  abs_f_interp_order,
2101  a_pressure,
2102  dtemp,
2103  a_vmr_list,
2104  f_grid,
2105  extpolfac);
2106  }
2107 
2108  // Now add to the right place in the absorption matrix.
2109 
2110  if (not do_jac) {
2111  for (Index ii = 0; ii < propmat_clearsky.nelem(); ii++) {
2112  propmat_clearsky[ii].Kjj() += abs_scalar_gas(ii, joker);
2113  }
2114  } else {
2115  for (Index isp = 0; isp < propmat_clearsky.nelem(); isp++) {
2116  propmat_clearsky[isp].Kjj() += abs_scalar_gas(isp, joker);
2117 
2118  for (Index iv = 0; iv < propmat_clearsky[isp].NumberOfFrequencies();
2119  iv++) {
2120  for (Index iq = 0; iq < jacobian_quantities_position.nelem(); iq++) {
2121  if (jacobian_quantities[jacobian_quantities_position[iq]] ==
2123  dpropmat_clearsky_dx[iq].Kjj()[iv] +=
2124  (dabs_scalar_gas_dt(isp, iv) - abs_scalar_gas(isp, iv)) / dt;
2125  } else if (is_frequency_parameter(
2126  jacobian_quantities
2127  [jacobian_quantities_position[iq]])) {
2128  dpropmat_clearsky_dx[iq].Kjj()[iv] +=
2129  (dabs_scalar_gas_df(isp, iv) - abs_scalar_gas(isp, iv)) / df;
2130  } else if (jacobian_quantities[jacobian_quantities_position[iq]] ==
2132  if (jacobian_quantities[jacobian_quantities_position[iq]]
2133  .QuantumIdentity()
2134  .Species() not_eq abs_lookup.GetSpeciesIndex(isp))
2135  continue; // FIXME?: Ignoring isotopologue???
2136 
2137  // WARNING: If CIA in list, this scales wrong by factor 2
2138  dpropmat_clearsky_dx[iq].Kjj()[iv] +=
2139  abs_scalar_gas(isp, iv) / a_vmr_list[isp];
2140  }
2141  }
2142  }
2143  }
2144  }
2145 }
2146 
2147 /* Workspace method: Doxygen documentation will be auto-generated */
2149  // WS Output:
2150  Tensor7& propmat_clearsky_field,
2151  Tensor6& nlte_source_field,
2152  // WS Input:
2153  const Index& atmfields_checked,
2154  const Vector& f_grid,
2155  const Index& stokes_dim,
2156  const Vector& p_grid,
2157  const Vector& lat_grid,
2158  const Vector& lon_grid,
2159  const Tensor3& t_field,
2160  const Tensor4& vmr_field,
2161  const EnergyLevelMap& nlte_field,
2162  const Tensor3& mag_u_field,
2163  const Tensor3& mag_v_field,
2164  const Tensor3& mag_w_field,
2165  const Agenda& abs_agenda,
2166  // WS Generic Input:
2167  const Vector& doppler,
2168  const Vector& los,
2169  const Verbosity& verbosity) {
2170  CREATE_OUT2;
2171  CREATE_OUT3;
2172 
2173  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2174  if (atmfields_checked != 1)
2175  throw runtime_error(
2176  "The atmospheric fields must be flagged to have "
2177  "passed a consistency check (atmfields_checked=1).");
2178 
2179  ArrayOfPropagationMatrix partial_abs;
2180  ArrayOfStokesVector partial_nlte_source,
2181  nlte_partial_source; // FIXME: This is not stored!
2183  ArrayOfStokesVector nlte;
2184  Vector a_vmr_list;
2185  EnergyLevelMap a_nlte_list;
2186 
2187  // Get the number of species from the leading dimension of vmr_field:
2188  const Index n_species = vmr_field.nbooks();
2189 
2190  // Number of frequencies:
2191  const Index n_frequencies = f_grid.nelem();
2192 
2193  // Number of pressure levels:
2194  const Index n_pressures = p_grid.nelem();
2195 
2196  // Number of latitude grid points (must be at least one):
2197  const Index n_latitudes = max(Index(1), lat_grid.nelem());
2198 
2199  // Number of longitude grid points (must be at least one):
2200  const Index n_longitudes = max(Index(1), lon_grid.nelem());
2201 
2202  // Check that doppler is empty or matches p_grid
2203  if (0 != doppler.nelem() && p_grid.nelem() != doppler.nelem()) {
2204  ostringstream os;
2205  os << "Variable doppler must either be empty, or match the dimension of "
2206  << "p_grid.";
2207  throw runtime_error(os.str());
2208  }
2209 
2210  // Resize output field.
2211  // The dimension in lat and lon must be at least one, even if these
2212  // grids are empty.
2213  out2 << " Creating propmat field with dimensions:\n"
2214  << " " << n_species << " gas species,\n"
2215  << " " << n_frequencies << " frequencies,\n"
2216  << " " << stokes_dim << " stokes dimension,\n"
2217  << " " << stokes_dim << " stokes dimension,\n"
2218  << " " << n_pressures << " pressures,\n"
2219  << " " << n_latitudes << " latitudes,\n"
2220  << " " << n_longitudes << " longitudes.\n";
2221 
2222  propmat_clearsky_field.resize(n_species,
2223  n_frequencies,
2224  stokes_dim,
2225  stokes_dim,
2226  n_pressures,
2227  n_latitudes,
2228  n_longitudes);
2229  if (nlte_field.Data().empty()) {
2230  out2 << " Creating source field with dimensions:\n"
2231  << " " << n_species << " gas species,\n"
2232  << " " << n_frequencies << " frequencies,\n"
2233  << " " << stokes_dim << " stokes dimension,\n"
2234  << " " << n_pressures << " pressures,\n"
2235  << " " << n_latitudes << " latitudes,\n"
2236  << " " << n_longitudes << " longitudes.\n";
2237  nlte_source_field.resize(n_species,
2238  n_frequencies,
2239  stokes_dim,
2240  n_pressures,
2241  n_latitudes,
2242  n_longitudes);
2243  } else {
2244  out2 << " Creating source field with dimensions:\n"
2245  << " " << 0 << " gas species,\n"
2246  << " " << 0 << " frequencies,\n"
2247  << " " << 0 << " stokes dimension,\n"
2248  << " " << 0 << " pressures,\n"
2249  << " " << 0 << " latitudes,\n"
2250  << " " << 0 << " longitudes.\n";
2251  nlte_source_field.resize(0, 0, 0, 0, 0, 0);
2252  }
2253 
2254  // We have to make a local copy of the Workspace and the agendas because
2255  // only non-reference types can be declared firstprivate in OpenMP
2256  Workspace l_ws(ws);
2257  Agenda l_abs_agenda(abs_agenda);
2258 
2259  String fail_msg;
2260  bool failed = false;
2261 
2262  // Make local copy of f_grid, so that we can apply Dopler if we want.
2263  Vector this_f_grid = f_grid;
2264 
2265  // Now we have to loop all points in the atmosphere:
2266  if (n_pressures)
2267 #pragma omp parallel for if (!arts_omp_in_parallel() && \
2268  n_pressures >= arts_omp_get_max_threads()) \
2269  firstprivate(l_ws, l_abs_agenda, this_f_grid) private(abs, \
2270  nlte, \
2271  partial_abs, \
2272  partial_nlte_source, \
2273  nlte_partial_source, \
2274  a_vmr_list)
2275  for (Index ipr = 0; ipr < n_pressures; ++ipr) // Pressure: ipr
2276  {
2277  // Skip remaining iterations if an error occurred
2278  if (failed) continue;
2279 
2280  // The try block here is necessary to correctly handle
2281  // exceptions inside the parallel region.
2282  try {
2283  Numeric a_pressure = p_grid[ipr];
2284 
2285  if (0 != doppler.nelem()) {
2286  this_f_grid = f_grid;
2287  this_f_grid += doppler[ipr];
2288  }
2289 
2290  {
2291  ostringstream os;
2292  os << " p_grid[" << ipr << "] = " << a_pressure << "\n";
2293  out3 << os.str();
2294  }
2295 
2296  for (Index ila = 0; ila < n_latitudes; ++ila) // Latitude: ila
2297  for (Index ilo = 0; ilo < n_longitudes; ++ilo) // Longitude: ilo
2298  {
2299  Numeric a_temperature = t_field(ipr, ila, ilo);
2300  a_vmr_list = vmr_field(Range(joker), ipr, ila, ilo);
2301  if (!nlte_field.Data().empty())
2302  a_nlte_list = nlte_field(ipr, ila, ilo);
2303 
2304  Vector this_rtp_mag(3, 0.);
2305 
2306  if (mag_u_field.npages() != 0) {
2307  this_rtp_mag[0] = mag_u_field(ipr, ila, ilo);
2308  }
2309  if (mag_v_field.npages() != 0) {
2310  this_rtp_mag[1] = mag_v_field(ipr, ila, ilo);
2311  }
2312  if (mag_w_field.npages() != 0) {
2313  this_rtp_mag[2] = mag_w_field(ipr, ila, ilo);
2314  }
2315 
2316  // Execute agenda to calculate local absorption.
2317  // Agenda input: f_index, a_pressure, a_temperature, a_vmr_list
2318  // Agenda output: abs, nlte
2320  abs,
2321  nlte,
2322  partial_abs,
2323  partial_nlte_source,
2324  nlte_partial_source,
2326  this_f_grid,
2327  this_rtp_mag,
2328  los,
2329  a_pressure,
2330  a_temperature,
2331  a_nlte_list,
2332  a_vmr_list,
2333  l_abs_agenda);
2334 
2335  // Verify, that the number of elements in abs matrix is
2336  // constistent with stokes_dim:
2337  if (abs.nelem() > 0) {
2338  if (stokes_dim != abs[0].StokesDimensions()) {
2339  ostringstream os;
2340  os << "propmat_clearsky_fieldCalc was called with stokes_dim = "
2341  << stokes_dim << ",\n"
2342  << "but the stokes_dim returned by the agenda is "
2343  << abs[0].StokesDimensions() << ".";
2344  throw runtime_error(os.str());
2345  }
2346  }
2347 
2348  // Verify, that the number of species in abs is
2349  // constistent with vmr_field:
2350  if (n_species != abs.nelem()) {
2351  ostringstream os;
2352  os << "The number of gas species in vmr_field is " << n_species
2353  << ",\n"
2354  << "but the number of species returned by the agenda is "
2355  << abs.nelem() << ".";
2356  throw runtime_error(os.str());
2357  }
2358 
2359  // Verify, that the number of frequencies in abs is
2360  // constistent with f_extent:
2361  if (abs.nelem() > 0) {
2362  if (n_frequencies != abs[0].NumberOfFrequencies()) {
2363  ostringstream os;
2364  os << "The number of frequencies desired is " << n_frequencies
2365  << ",\n"
2366  << "but the number of frequencies returned by the agenda is "
2367  << abs[0].NumberOfFrequencies() << ".";
2368  throw runtime_error(os.str());
2369  }
2370  }
2371 
2372  // Store the result in output field.
2373  // We have to transpose abs, because the dimensions of the
2374  // two variables are:
2375  // abs_field: [abs_species, f_grid, stokes, stokes, p_grid, lat_grid, lon_grid]
2376  // abs: [abs_species][f_grid, stokes, stokes]
2377  for (Index i = 0; i < abs.nelem(); i++) {
2378  abs[i].GetTensor3(propmat_clearsky_field(
2379  i, joker, joker, joker, ipr, ila, ilo));
2380 
2381  if (not nlte_field.Data().empty()) {
2382  //If some are NLTE and others not, this might be bad...
2383  nlte_source_field(i, joker, joker, ipr, ila, ilo) =
2384  nlte[i].Data()(0, 0, joker, joker);
2385  }
2386  }
2387  }
2388  } catch (const std::runtime_error& e) {
2389 #pragma omp critical(propmat_clearsky_fieldCalc_fail)
2390  {
2391  fail_msg = e.what();
2392  failed = true;
2393  }
2394  }
2395  }
2396 
2397  if (failed) throw runtime_error(fail_msg);
2398 }
2399 
2400 /* Workspace method: Doxygen documentation will be auto-generated */
2402  const GasAbsLookup& abs_lookup,
2403  const Verbosity&) {
2404  const Vector& lookup_f_grid = abs_lookup.GetFgrid();
2405  f_grid.resize(lookup_f_grid.nelem());
2406  f_grid = lookup_f_grid;
2407 }
2408 
2409 /* Workspace method: Doxygen documentation will be auto-generated */
2411  const GasAbsLookup& abs_lookup,
2412  const Verbosity&) {
2413  const Vector& lookup_p_grid = abs_lookup.GetPgrid();
2414  p_grid.resize(lookup_p_grid.nelem());
2415  p_grid = lookup_p_grid;
2416 }
2417 
2419 
2441 Numeric calc_lookup_error( // Parameters for lookup table:
2442  Workspace& ws,
2443  const GasAbsLookup& al,
2444  const Index& abs_p_interp_order,
2445  const Index& abs_t_interp_order,
2446  const Index& abs_nls_interp_order,
2447  const bool ignore_errors,
2448  // Parameters for LBL:
2449  const Agenda& abs_xsec_agenda,
2450  // Parameters for both:
2451  const Numeric& local_p,
2452  const Numeric& local_t,
2453  const Vector& local_vmrs,
2454  const Verbosity& verbosity) {
2455  // Allocate some matrices. I also tried allocating these (and the
2456  // vectors below) outside, but there was no significant speed
2457  // advantage. (I guess the LBL calculation is expensive enough to
2458  // make the extra time of allocation here insignificant.)
2459  Matrix sga_tab; // Absorption, dimension [n_species,n_f_grid]:
2460  const EnergyLevelMap local_nlte_dummy;
2461 
2462  // Do lookup table first:
2463 
2464  try {
2465  // Absorption, dimension [n_species,n_f_grid]:
2466  // Output variable: sga_tab
2467  al.Extract(sga_tab,
2468  abs_p_interp_order,
2469  abs_t_interp_order,
2470  abs_nls_interp_order,
2471  0, // f_interp_order
2472  local_p,
2473  local_t,
2474  local_vmrs,
2475  al.f_grid,
2476  0.0); // Extpolfac
2477  } catch (const std::runtime_error& x) {
2478  // If ignore_errors is set to true, then we mark this case for
2479  // skipping, and ignore the exceptions.
2480  // Otherwise, we re-throw the exception.
2481  if (ignore_errors)
2482  return -1;
2483  else
2484  throw runtime_error(x.what());
2485  }
2486 
2487  // Get number of frequencies. (We cannot do this earlier, since we
2488  // get it from the output of al.Extract.)
2489  const Index n_f = sga_tab.ncols();
2490 
2491  // Allocate some vectors with this dimension:
2492  Vector abs_tab(n_f);
2493  Vector abs_lbl(n_f, 0.0);
2494  Vector abs_rel_diff(n_f);
2495 
2496  // Sum up for all species, to get total absorption:
2497  for (Index i = 0; i < n_f; ++i) abs_tab[i] = sga_tab(joker, i).sum();
2498 
2499  // Now get absorption line-by-line.
2500 
2501  // Variable to hold result of absorption calculation:
2502  ArrayOfPropagationMatrix propmat_clearsky;
2503  ArrayOfStokesVector nlte_source;
2504  ArrayOfPropagationMatrix dpropmat_clearsky_dx;
2505  ArrayOfStokesVector dnlte_dx_source, nlte_dsource_dx;
2506  ArrayOfMatrix d;
2507  const ArrayOfRetrievalQuantity jacobian_quantities(0);
2508  Index propmat_clearsky_checked = 1,
2509  nlte_do = 0; // FIXME: OLE: Properly pass this through?
2510 
2511  // Initialize propmat_clearsky:
2512  propmat_clearskyInit(propmat_clearsky,
2513  nlte_source,
2514  dpropmat_clearsky_dx,
2515  dnlte_dx_source,
2516  nlte_dsource_dx,
2517  al.species,
2518  jacobian_quantities,
2519  al.f_grid,
2520  1, // Stokes dimension
2521  propmat_clearsky_checked,
2522  nlte_do,
2523  verbosity);
2524 
2525  // Add result of LBL calculation to propmat_clearsky:
2527  propmat_clearsky,
2528  nlte_source,
2529  dpropmat_clearsky_dx,
2530  dnlte_dx_source,
2531  nlte_dsource_dx,
2532  al.f_grid,
2533  al.species,
2534  jacobian_quantities,
2535  local_p,
2536  local_t,
2537  local_nlte_dummy,
2538  local_vmrs,
2539  abs_xsec_agenda,
2540  verbosity);
2541 
2542  // Argument 0 above is the Doppler shift (usually
2543  // rtp_doppler). Should be zero in this case.
2544 
2545  // Sum up for all species, to get total absorption:
2546  for (auto& pm : propmat_clearsky) abs_lbl += pm.Kjj();
2547 
2548  // Ok. What we have to compare is abs_tab and abs_lbl.
2549 
2550  assert(abs_tab.nelem() == n_f);
2551  assert(abs_lbl.nelem() == n_f);
2552  assert(abs_rel_diff.nelem() == n_f);
2553  for (Index i = 0; i < n_f; ++i) {
2554  // Absolute value of relative difference in percent:
2555  abs_rel_diff[i] = fabs((abs_tab[i] - abs_lbl[i]) / abs_lbl[i] * 100);
2556  }
2557 
2558  // Maximum of this:
2559  Numeric max_abs_rel_diff = max(abs_rel_diff);
2560 
2561  // cout << "ma " << max_abs_rel_diff << "\n";
2562 
2563  return max_abs_rel_diff;
2564 }
2565 
2566 /* Workspace method: Doxygen documentation will be auto-generated */
2567 void abs_lookupTestAccuracy( // Workspace reference:
2568  Workspace& ws,
2569  // WS Input:
2570  const GasAbsLookup& abs_lookup,
2571  const Index& abs_lookup_is_adapted,
2572  const Index& abs_p_interp_order,
2573  const Index& abs_t_interp_order,
2574  const Index& abs_nls_interp_order,
2575  const Agenda& abs_xsec_agenda,
2576  // Verbosity object:
2577  const Verbosity& verbosity) {
2578  CREATE_OUT2;
2579 
2580  const GasAbsLookup& al = abs_lookup;
2581 
2582  // Check if the table has been adapted:
2583  if (1 != abs_lookup_is_adapted)
2584  throw runtime_error(
2585  "Gas absorption lookup table must be adapted,\n"
2586  "use method abs_lookupAdapt.");
2587 
2588  // Some important sizes:
2589  const Index n_nls = al.nonlinear_species.nelem();
2590  const Index n_species = al.species.nelem();
2591  // const Index n_f = al.f_grid.nelem();
2592  const Index n_p = al.log_p_grid.nelem();
2593 
2594  if (n_nls <= 0) {
2595  ostringstream os;
2596  os << "This function currently works only with lookup tables\n"
2597  << "containing nonlinear species.";
2598  throw runtime_error(os.str());
2599  }
2600 
2601  // If there are nonlinear species, then at least one species must be
2602  // H2O. We will use that to perturb in the case of nonlinear
2603  // species.
2604  Index h2o_index = -1;
2605  if (n_nls > 0) {
2606  h2o_index = find_first_species_tg(al.species,
2608 
2609  // This is a runtime error, even though it would be more logical
2610  // for it to be an assertion, since it is an internal check on
2611  // the table. The reason is that it is somewhat awkward to check
2612  // for this in other places.
2613  if (h2o_index == -1) {
2614  ostringstream os;
2615  os << "With nonlinear species, at least one species must be a H2O species.";
2616  throw runtime_error(os.str());
2617  }
2618  }
2619 
2620  // Check temperature interpolation
2621 
2622  Vector inbet_t_pert(al.t_pert.nelem() - 1);
2623  for (Index i = 0; i < inbet_t_pert.nelem(); ++i)
2624  inbet_t_pert[i] = (al.t_pert[i] + al.t_pert[i + 1]) / 2.0;
2625 
2626  // To store the temperature error, which we define as the maximum of
2627  // the absolute value of the relative difference between LBL and
2628  // lookup table, in percent.
2629  Numeric err_t = -999;
2630 
2631 #pragma omp parallel for if (!arts_omp_in_parallel())
2632  for (Index pi = 0; pi < n_p; ++pi)
2633  for (Index ti = 0; ti < inbet_t_pert.nelem(); ++ti) {
2634  // Find local conditions:
2635 
2636  // Pressure:
2637  Numeric local_p = al.p_grid[pi];
2638 
2639  // Temperature:
2640  Numeric local_t = al.t_ref[pi] + inbet_t_pert[ti];
2641 
2642  // VMRs:
2643  Vector local_vmrs = al.vmrs_ref(joker, pi);
2644 
2645  // Watch out, the table probably does not have an absorption
2646  // value for exactly the reference H2O profile. We multiply
2647  // with the first perturbation.
2648  local_vmrs[h2o_index] *= al.nls_pert[0];
2649 
2650  Numeric max_abs_rel_diff =
2651  calc_lookup_error(ws,
2652  // Parameters for lookup table:
2653  al,
2654  abs_p_interp_order,
2655  abs_t_interp_order,
2656  abs_nls_interp_order,
2657  true, // ignore errors
2658  // Parameters for LBL:
2659  abs_xsec_agenda,
2660  // Parameters for both:
2661  local_p,
2662  local_t,
2663  local_vmrs,
2664  verbosity);
2665 
2666  // cout << "ma " << max_abs_rel_diff << "\n";
2667 
2668  //Critical directive here is necessary, because all threads
2669  //access the same variable.
2670 #pragma omp critical(abs_lookupTestAccuracy_piti)
2671  {
2672  if (max_abs_rel_diff > err_t) err_t = max_abs_rel_diff;
2673  }
2674 
2675  } // end parallel for loop
2676 
2677  // Check H2O interpolation
2678 
2679  Vector inbet_nls_pert(al.nls_pert.nelem() - 1);
2680  for (Index i = 0; i < inbet_nls_pert.nelem(); ++i)
2681  inbet_nls_pert[i] = (al.nls_pert[i] + al.nls_pert[i + 1]) / 2.0;
2682 
2683  // To store the H2O error, which we define as the maximum of
2684  // the absolute value of the relative difference between LBL and
2685  // lookup table, in percent.
2686  Numeric err_nls = -999;
2687 
2688 #pragma omp parallel for if (!arts_omp_in_parallel())
2689  for (Index pi = 0; pi < n_p; ++pi)
2690  for (Index ni = 0; ni < inbet_nls_pert.nelem(); ++ni) {
2691  // Find local conditions:
2692 
2693  // Pressure:
2694  Numeric local_p = al.p_grid[pi];
2695 
2696  // Temperature:
2697 
2698  // Watch out, the table probably does not have an absorption
2699  // value for exactly the reference temperature. We add
2700  // the first perturbation.
2701 
2702  Numeric local_t = al.t_ref[pi] + al.t_pert[0];
2703 
2704  // VMRs:
2705  Vector local_vmrs = al.vmrs_ref(joker, pi);
2706 
2707  // Now we have to modify the H2O VMR according to nls_pert:
2708  local_vmrs[h2o_index] *= inbet_nls_pert[ni];
2709 
2710  Numeric max_abs_rel_diff =
2711  calc_lookup_error(ws,
2712  // Parameters for lookup table:
2713  al,
2714  abs_p_interp_order,
2715  abs_t_interp_order,
2716  abs_nls_interp_order,
2717  true, // ignore errors
2718  // Parameters for LBL:
2719  abs_xsec_agenda,
2720  // Parameters for both:
2721  local_p,
2722  local_t,
2723  local_vmrs,
2724  verbosity);
2725 
2726  //Critical directive here is necessary, because all threads
2727  //access the same variable.
2728 #pragma omp critical(abs_lookupTestAccuracy_pini)
2729  {
2730  if (max_abs_rel_diff > err_nls) err_nls = max_abs_rel_diff;
2731  }
2732 
2733  } // end parallel for loop
2734 
2735  // Check pressure interpolation
2736 
2737  // IMPORTANT: This does not test the pure pressure interpolation,
2738  // unless we have constant reference profiles for T and
2739  // H2O. Otherwise we have T and H2O interpolation mixed in.
2740 
2741  Vector inbet_p_grid(n_p - 1);
2742  Vector inbet_t_ref(n_p - 1);
2743  Matrix inbet_vmrs_ref(n_species, n_p - 1);
2744  for (Index i = 0; i < inbet_p_grid.nelem(); ++i) {
2745  inbet_p_grid[i] = exp((al.log_p_grid[i] + al.log_p_grid[i + 1]) / 2.0);
2746  inbet_t_ref[i] = (al.t_ref[i] + al.t_ref[i + 1]) / 2.0;
2747  for (Index j = 0; j < n_species; ++j)
2748  inbet_vmrs_ref(j, i) = (al.vmrs_ref(j, i) + al.vmrs_ref(j, i + 1)) / 2.0;
2749  }
2750 
2751  // To store the interpolation error, which we define as the maximum of
2752  // the absolute value of the relative difference between LBL and
2753  // lookup table, in percent.
2754  Numeric err_p = -999;
2755 
2756 #pragma omp parallel for if (!arts_omp_in_parallel())
2757  for (Index pi = 0; pi < n_p - 1; ++pi) {
2758  // Find local conditions:
2759 
2760  // Pressure:
2761  Numeric local_p = inbet_p_grid[pi];
2762 
2763  // Temperature:
2764 
2765  // Watch out, the table probably does not have an absorption
2766  // value for exactly the reference temperature. We add
2767  // the first perturbation.
2768 
2769  Numeric local_t = inbet_t_ref[pi] + al.t_pert[0];
2770 
2771  // VMRs:
2772  Vector local_vmrs = inbet_vmrs_ref(joker, pi);
2773 
2774  // Watch out, the table probably does not have an absorption
2775  // value for exactly the reference H2O profile. We multiply
2776  // with the first perturbation.
2777  local_vmrs[h2o_index] *= al.nls_pert[0];
2778 
2779  Numeric max_abs_rel_diff = calc_lookup_error(ws,
2780  // Parameters for lookup table:
2781  al,
2782  abs_p_interp_order,
2783  abs_t_interp_order,
2784  abs_nls_interp_order,
2785  true, // ignore errors
2786  // Parameters for LBL:
2787  abs_xsec_agenda,
2788  // Parameters for both:
2789  local_p,
2790  local_t,
2791  local_vmrs,
2792  verbosity);
2793 
2794  //Critical directive here is necessary, because all threads
2795  //access the same variable.
2796 #pragma omp critical(abs_lookupTestAccuracy_pi)
2797  {
2798  if (max_abs_rel_diff > err_p) err_p = max_abs_rel_diff;
2799  }
2800  }
2801 
2802  // Check total error
2803 
2804  // To store the interpolation error, which we define as the maximum of
2805  // the absolute value of the relative difference between LBL and
2806  // lookup table, in percent.
2807  Numeric err_tot = -999;
2808 
2809 #pragma omp parallel for if (!arts_omp_in_parallel())
2810  for (Index pi = 0; pi < n_p - 1; ++pi)
2811  for (Index ti = 0; ti < inbet_t_pert.nelem(); ++ti)
2812  for (Index ni = 0; ni < inbet_nls_pert.nelem(); ++ni) {
2813  // Find local conditions:
2814 
2815  // Pressure:
2816  Numeric local_p = inbet_p_grid[pi];
2817 
2818  // Temperature:
2819  Numeric local_t = inbet_t_ref[pi] + inbet_t_pert[ti];
2820 
2821  // VMRs:
2822  Vector local_vmrs = inbet_vmrs_ref(joker, pi);
2823 
2824  // Multiply with perturbation.
2825  local_vmrs[h2o_index] *= inbet_nls_pert[ni];
2826 
2827  Numeric max_abs_rel_diff =
2828  calc_lookup_error(ws,
2829  // Parameters for lookup table:
2830  al,
2831  abs_p_interp_order,
2832  abs_t_interp_order,
2833  abs_nls_interp_order,
2834  true, // ignore errors
2835  // Parameters for LBL:
2836  abs_xsec_agenda,
2837  // Parameters for both:
2838  local_p,
2839  local_t,
2840  local_vmrs,
2841  verbosity);
2842 
2843  //Critical directive here is necessary, because all threads
2844  //access the same variable.
2845 #pragma omp critical(abs_lookupTestAccuracy_pitini)
2846  {
2847  if (max_abs_rel_diff > err_tot) {
2848  err_tot = max_abs_rel_diff;
2849 
2850  // cout << "New max error: pi, ti, ni, err_tot:\n"
2851  // << pi << ", " << ti << ", " << ni << ", " << err_tot << "\n";
2852  }
2853  }
2854  }
2855 
2856  out2 << " Max. of absolute value of relative error in percent:\n"
2857  << " Note: Unless you have constant reference profiles, the\n"
2858  << " pressure interpolation error will have other errors mixed in.\n"
2859  << " Temperature interpolation: " << err_t << "%\n"
2860  << " H2O (NLS) interpolation: " << err_nls << "%\n"
2861  << " Pressure interpolation: " << err_p << "%\n"
2862  << " Total error: " << err_tot << "%\n";
2863 
2864  // Check pressure interpolation
2865 
2866  // assert(p_grid.nelem()==log_p_grid.nelem()); // Make sure that log_p_grid is initialized.
2867  // Vector inbet_log_p_grid(log_p_grid.nelem()-1)
2868  // for (Index i=0; i<log_p_grid.nelem()-1; ++i)
2869  // {
2870  // inbet_log_p_grid[i] = (log_p_grid[i]+log_p_grid[i+1])/2.0;
2871  // }
2872 
2873  // for (Index pi=0; pi<inbet_log_p_grid.nelem(); ++pi)
2874  // {
2875  // for (Index pt=0; pt<)
2876  // }
2877 }
2878 
2879 /* Workspace method: Doxygen documentation will be auto-generated */
2880 void abs_lookupTestAccMC( // Workspace reference:
2881  Workspace& ws,
2882  // WS Input:
2883  const GasAbsLookup& abs_lookup,
2884  const Index& abs_lookup_is_adapted,
2885  const Index& abs_p_interp_order,
2886  const Index& abs_t_interp_order,
2887  const Index& abs_nls_interp_order,
2888  const Index& mc_seed,
2889  const Agenda& abs_xsec_agenda,
2890  // Verbosity object:
2891  const Verbosity& verbosity) {
2892  CREATE_OUT2;
2893  CREATE_OUT3;
2894 
2895  const GasAbsLookup& al = abs_lookup;
2896 
2897  // Check if the table has been adapted:
2898  if (1 != abs_lookup_is_adapted)
2899  throw runtime_error(
2900  "Gas absorption lookup table must be adapted,\n"
2901  "use method abs_lookupAdapt.");
2902 
2903  // Some important sizes:
2904  const Index n_nls = al.nonlinear_species.nelem();
2905  const Index n_species = al.species.nelem();
2906 
2907  if (n_nls <= 0) {
2908  ostringstream os;
2909  os << "This function currently works only with lookup tables\n"
2910  << "containing nonlinear species.";
2911  throw runtime_error(os.str());
2912  }
2913 
2914  // If there are nonlinear species, then at least one species must be
2915  // H2O. We will use that to perturb in the case of nonlinear
2916  // species.
2917  Index h2o_index = -1;
2918  if (n_nls > 0) {
2919  h2o_index = find_first_species_tg(al.species,
2921 
2922  // This is a runtime error, even though it would be more logical
2923  // for it to be an assertion, since it is an internal check on
2924  // the table. The reason is that it is somewhat awkward to check
2925  // for this in other places.
2926  if (h2o_index == -1) {
2927  ostringstream os;
2928  os << "With nonlinear species, at least one species must be a H2O species.";
2929  throw runtime_error(os.str());
2930  }
2931  }
2932 
2933  // How many MC cases to run between each convergence check.
2934  // (It is important for parallelization that this is not too small.)
2935  const Index chunksize = 100;
2936 
2937  //Random Number generator:
2938  Rng rng;
2939  rng.seed(mc_seed, verbosity);
2940  // rng.draw() will draw a double from the uniform distribution [0,1).
2941 
2942  // (Log) Pressure range:
2943  const Numeric lp_max = al.log_p_grid[0];
2944  const Numeric lp_min = al.log_p_grid[al.log_p_grid.nelem() - 1];
2945 
2946  // T perturbation range (additive):
2947  const Numeric dT_min = al.t_pert[0];
2948  const Numeric dT_max = al.t_pert[al.t_pert.nelem() - 1];
2949 
2950  // H2O perturbation range (scaling):
2951  const Numeric dh2o_min = al.nls_pert[0];
2952  const Numeric dh2o_max = al.nls_pert[al.nls_pert.nelem() - 1];
2953 
2954  // We are creating all random numbers for the chunk beforehand, to avoid the
2955  // problem that random number generators in different threads would need
2956  // different seeds to produce independent random numbers.
2957  // (I prefer this solution to the one of having the rng inside the threads,
2958  // because it ensures that the result does not depend on the the number of CPUs.)
2959  Vector rand_lp(chunksize);
2960  Vector rand_dT(chunksize);
2961  Vector rand_dh2o(chunksize);
2962 
2963  // Store the errors for one chunk:
2964  Vector max_abs_rel_diff(chunksize);
2965 
2966  // Flag to break our MC calculation loop eventually
2967  bool keep_looping = true;
2968 
2969  // Total mean and standard deviation. (Is updated after each chunk.)
2970  Numeric total_mean;
2971  Numeric total_std;
2972  Index N_chunk = 0;
2973  while (keep_looping) {
2974  ++N_chunk;
2975 
2976  for (Index i = 0; i < chunksize; ++i) {
2977  // A random pressure, temperature perturbation, and H2O perturbation,
2978  // all with flat PDF between min and max:
2979  rand_lp[i] = rng.draw() * (lp_max - lp_min) + lp_min;
2980  rand_dT[i] = rng.draw() * (dT_max - dT_min) + dT_min;
2981  rand_dh2o[i] = rng.draw() * (dh2o_max - dh2o_min) + dh2o_min;
2982  }
2983 
2984  for (Index i = 0; i < chunksize; ++i) {
2985  // The pressure we work with here:
2986  const Numeric this_lp = rand_lp[i];
2987 
2988  // Now we have to interpolate t_ref and vmrs_ref to this
2989  // pressure, so that we can apply the dT and dh2o perturbations.
2990 
2991  // Pressure grid positions:
2992  ArrayOfGridPosPoly pgp(1);
2993  gridpos_poly(pgp, al.log_p_grid, this_lp, abs_p_interp_order);
2994 
2995  // Pressure interpolation weights:
2996  Vector pitw;
2997  pitw.resize(abs_p_interp_order + 1);
2998  interpweights(pitw, pgp[0]);
2999 
3000  // Interpolated temperature:
3001  const Numeric this_t_ref = interp(pitw, al.t_ref, pgp[0]);
3002 
3003  // Interpolated VMRs:
3004  Vector these_vmrs(n_species);
3005  for (Index j = 0; j < n_species; ++j) {
3006  these_vmrs[j] = interp(pitw, al.vmrs_ref(j, Range(joker)), pgp[0]);
3007  }
3008 
3009  // Now get the actual p, T and H2O values:
3010  const Numeric this_p = exp(this_lp);
3011  const Numeric this_t = this_t_ref + rand_dT[i];
3012  these_vmrs[h2o_index] *= rand_dh2o[i];
3013 
3014  // cout << "p, T, H2O: " << this_p << ", " << this_t << ", " << these_vmrs[h2o_index] << "\n";
3015 
3016  // Get error between table and LBL calculation for these conditions:
3017 
3018  max_abs_rel_diff[i] = calc_lookup_error(ws,
3019  // Parameters for lookup table:
3020  al,
3021  abs_p_interp_order,
3022  abs_t_interp_order,
3023  abs_nls_interp_order,
3024  true, // ignore errors
3025  // Parameters for LBL:
3026  abs_xsec_agenda,
3027  // Parameters for both:
3028  this_p,
3029  this_t,
3030  these_vmrs,
3031  verbosity);
3032  // cout << "max_abs_rel_diff[" << i << "] = " << max_abs_rel_diff[i] << "\n";
3033  }
3034 
3035  // Calculate Mean of the last batch.
3036 
3037  // Total number of valid points in the chunk (not counting negative values,
3038  // which result from failed calculations at the edges of the table.)
3039  Index N = 0;
3040  // Mean (initially sum of all values):
3041  Numeric mean = 0;
3042  for (Index i = 0; i < chunksize; ++i) {
3043  const Numeric x = max_abs_rel_diff[i];
3044  if (x > 0) {
3045  ++N;
3046  mean += x;
3047  }
3048  // else
3049  // {
3050  // cout << "Negative value ignored.\n";
3051  // }
3052  }
3053  // Calculate mean by dividing sum by number of valid points:
3054  mean = mean / (Numeric)N;
3055 
3056  // Now calculate standard deviation:
3057 
3058  // Variance (initially sum of squared differences)
3059  Numeric variance = 0;
3060  for (Index i = 0; i < chunksize; ++i) {
3061  const Numeric x = max_abs_rel_diff[i];
3062  if (x > 0) {
3063  variance += (x - mean) * (x - mean);
3064  }
3065  }
3066  // Divide by N to really calculate variance:
3067  variance = variance / (Numeric)N;
3068 
3069  // cout << "Mean = " << mean << " Std = " << std_dev << "\n";
3070 
3071  if (N_chunk == 1) {
3072  total_mean = mean;
3073  total_std = sqrt(variance);
3074  } else {
3075  const Numeric old_mean = total_mean;
3076 
3077  // This formula assimilates the new chunk mean into the total mean:
3078  total_mean =
3079  (total_mean * ((Numeric)(N_chunk - 1)) + mean) / (Numeric)N_chunk;
3080 
3081  // Do the same for the standard deviation.
3082  // First get rid of the square root:
3083  total_std = total_std * total_std;
3084  // Now multiply with old normalisation:
3085  total_std *= (Numeric)(N_chunk - 1);
3086  // Now add the new sigma
3087  total_std += variance;
3088  // Divide by the new normalisation:
3089  total_std /= (Numeric)N_chunk;
3090  // And finally take the square root:
3091  total_std = sqrt(total_std);
3092 
3093  // Stop the chunk loop if desired accuracy has been reached.
3094  // We take 1% here, no point in trying to be more accurate!
3095  if (abs(total_mean - old_mean) < total_mean / 100) keep_looping = false;
3096  }
3097 
3098  // cout << " Chunk " << N_chunk << ": Mean estimate = " << total_mean
3099  // << " Std estimate = " << total_std << "\n";
3100 
3101  out3 << " Chunk " << N_chunk << ": Mean estimate = " << total_mean
3102  << " Std estimate = " << total_std << "\n";
3103 
3104  } // End of "keep_looping" loop that runs over the chunks
3105 
3106  out2 << " Mean relative error: " << total_mean << "%\n"
3107  << " Standard deviation: " << total_std << "%\n";
3108 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void jacobianAddAbsSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition: m_jacobian.cc:159
ArrayOfIndex equivalent_propmattype_indexes(const ArrayOfRetrievalQuantity &js)
Returns a list of positions for the derivatives in Propagation Matrix calculations.
Definition: jacobian.cc:1099
#define N
Definition: rng.cc:164
void propmat_clearsky_fieldCalc(Workspace &ws, Tensor7 &propmat_clearsky_field, Tensor6 &nlte_source_field, const Index &atmfields_checked, const Vector &f_grid, const Index &stokes_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor4 &vmr_field, const EnergyLevelMap &nlte_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Agenda &abs_agenda, const Vector &doppler, const Vector &los, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearsky_fieldCalc.
void abs_lookupInit(GasAbsLookup &x, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupInit.
Definition: m_abs_lookup.cc:50
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:195
Vector nls_pert
The vector of perturbations for the VMRs of the nonlinear species.
QuantumIdentifier::QType Index LowerQuantumNumbers Species
Vector t_pert
The vector of temperature perturbations [K].
void choose_abs_nls(ArrayOfArrayOfSpeciesTag &abs_nls, const ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &verbosity)
Choose species for abs_nls.
Declarations having to do with the four output streams.
void parse_atmcompact_speciestype(String &species_type, const String &field_name, const String &delim)
Definition: cloudbox.cc:848
void abs_speciesInit(ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &)
WORKSPACE METHOD: abs_speciesInit.
void abs_lookupTestAccuracy(Workspace &ws, const GasAbsLookup &abs_lookup, const Index &abs_lookup_is_adapted, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Agenda &abs_xsec_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupTestAccuracy.
void check_abs_species(const ArrayOfArrayOfSpeciesTag &abs_species)
Check the correctness of abs_species.
void abs_speciesSet(ArrayOfArrayOfSpeciesTag &abs_species, Index &abs_xsec_agenda_checked, Index &propmat_clearsky_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesSet.
The Vector class.
Definition: matpackI.h:860
void f_gridFromGasAbsLookup(Vector &f_grid, const GasAbsLookup &abs_lookup, const Verbosity &)
WORKSPACE METHOD: f_gridFromGasAbsLookup.
#define abs(x)
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
const Tensor4 & Data() const noexcept
Energy level type.
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1312
void abs_lookupSetupBatch(Vector &abs_p, Vector &abs_t, Vector &abs_t_pert, Matrix &abs_vmrs, ArrayOfArrayOfSpeciesTag &abs_nls, Vector &abs_nls_pert, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfGriddedField4 &batch_fields, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Index &atmosphere_dim, const Numeric &p_step10, const Numeric &t_step, const Numeric &h2o_step, const Vector &extremes, const Index &robust, const Index &check_gridnames, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupSetupBatch.
The Tensor4 class.
Definition: matpackIV.h:421
The range class.
Definition: matpackI.h:160
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:402
String get_tag_group_name(const ArrayOfSpeciesTag &tg)
Return the name of a tag group as a string.
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:60
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
bool is_zeeman(const ArrayOfSpeciesTag &tg)
Is this a Zeeman tag group?
void abs_speciesAdd(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, Index &abs_xsec_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesAdd.
The Tensor7 class.
Definition: matpackVII.h:2382
void VectorInsertGridPoints(Vector &og, const Vector &ingrid, const Vector &points, const Verbosity &verbosity)
WORKSPACE METHOD: VectorInsertGridPoints.
ArrayOfIndex nonlinear_species
The species tags with non-linear treatment.
void Adapt(const ArrayOfArrayOfSpeciesTag &current_species, ConstVectorView current_f_grid, const Verbosity &verbosity)
Adapt lookup table to current calculation.
#define min(a, b)
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Index GetSpeciesIndex(const Index &isp) const
Vector log_p_grid
The natural log of the pressure grid.
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:63
void atmfields_checkedCalc(Index &atmfields_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const SpeciesAuxData &partition_functions, const Index &abs_f_interp_order, const Index &negative_vmr_ok, const Index &bad_partition_functions_ok, const Verbosity &)
WORKSPACE METHOD: atmfields_checkedCalc.
Definition: m_checked.cc:125
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
bool is_unique(const ArrayOfIndex &x)
Checks if an ArrayOfIndex is unique, i.e., has no duplicate values.
Definition: logic.cc:274
void choose_abs_nls_pert(Vector &abs_nls_pert, ConstVectorView refprof, ConstVectorView minprof, ConstVectorView maxprof, const Numeric &step, const Index &p_interp_order, const Index &nls_interp_order, const Verbosity &verbosity)
Chose the H2O perturbations abs_nls_pert.
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1120
Tensor4 xsec
Absorption cross sections.
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
The Tensor3 class.
Definition: matpackIII.h:339
The global header file for ARTS.
Matrix vmrs_ref
The reference VMR profiles.
Vector f_grid
The frequency grid [Hz].
void abs_lookupSetupWide(Vector &abs_p, Vector &abs_t, Vector &abs_t_pert, Matrix &abs_vmrs, ArrayOfArrayOfSpeciesTag &abs_nls, Vector &abs_nls_pert, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Numeric &p_min, const Numeric &p_max, const Numeric &p_step10, const Numeric &t_min, const Numeric &t_max, const Numeric &h2o_min, const Numeric &h2o_max, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupSetupWide.
_CS_string_type str() const
Definition: sstream.h:491
void propmat_clearskyAddFromLookup(ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const GasAbsLookup &abs_lookup, const Index &abs_lookup_is_adapted, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Index &abs_f_interp_order, const Vector &f_grid, const Numeric &a_pressure, const Numeric &a_temperature, const Vector &a_vmr_list, const ArrayOfRetrievalQuantity &jacobian_quantities, const Numeric &extpolfac, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddFromLookup.
Declarations for agendas.
void abs_lookupAdapt(GasAbsLookup &abs_lookup, Index &abs_lookup_is_adapted, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &f_grid, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupAdapt.
ArrayOfIndex idx
Numeric calc_lookup_error(Workspace &ws, const GasAbsLookup &al, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const bool ignore_errors, const Agenda &abs_xsec_agenda, const Numeric &local_p, const Numeric &local_t, const Vector &local_vmrs, const Verbosity &verbosity)
Compare lookup and LBL calculation.
void abs_lookupCalc(Workspace &ws, GasAbsLookup &abs_lookup, Index &abs_lookup_is_adapted, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfSpeciesTag &abs_nls, const Vector &f_grid, const Vector &abs_p, const Matrix &abs_vmrs, const Vector &abs_t, const Vector &abs_t_pert, const Vector &abs_nls_pert, const Agenda &abs_xsec_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupCalc.
Definition: m_abs_lookup.cc:60
void p_gridFromGasAbsLookup(Vector &p_grid, const GasAbsLookup &abs_lookup, const Verbosity &)
WORKSPACE METHOD: p_gridFromGasAbsLookup.
void propmat_clearskyAddOnTheFly(Workspace &ws, ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfStokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_dx_source, ArrayOfStokesVector &nlte_dsource_dx, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const EnergyLevelMap &rtp_nlte, const Vector &rtp_vmr, const Agenda &abs_xsec_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddOnTheFly.
Definition: m_abs.cc:1504
#define CREATE_OUT1
Definition: messages.h:205
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
An absorption lookup table.
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1304
const Vector & GetFgrid() const
const Joker joker
#define temp
member functions of the Rng class and gsl_rng code
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
Declarations for the gas absorption lookup table.
void AtmFieldsAndParticleBulkPropFieldFromCompact(Vector &p_grid, Vector &lat_grid, Vector &lon_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, Tensor4 &particle_bulkprop_field, ArrayOfString &particle_bulkprop_names, const ArrayOfArrayOfSpeciesTag &abs_species, const GriddedField4 &atm_fields_compact, const Index &atmosphere_dim, const String &delim, const Numeric &p_min, const Index &check_gridnames, const Verbosity &)
WORKSPACE METHOD: AtmFieldsAndParticleBulkPropFieldFromCompact.
Declarations required for the calculation of absorption coefficients.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
void choose_abs_t_pert(Vector &abs_t_pert, ConstVectorView abs_t, ConstVectorView tmin, ConstVectorView tmax, const Numeric &step, const Index &p_interp_order, const Index &t_interp_order, const Verbosity &verbosity)
Chose the temperature perturbations abs_t_pert.
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1279
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:531
Vector t_ref
The reference temperature profile [K].
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
bool supports_lookup(const ArrayOfRetrievalQuantity &js)
Returns if the array supports lookup table derivatives.
Definition: jacobian.cc:1224
String get_species_name(const ArrayOfSpeciesTag &tg)
Return the species of a tag group as a string.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Vector p_grid
The pressure grid for the table [Pa].
Definition: rng.h:554
Header file for interpolation_poly.cc.
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:466
Index find_next_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec, const Index &start)
Find next occurrence of species in tag groups.
void Extract(Matrix &sga, const Index &p_interp_order, const Index &t_interp_order, const Index &h2o_interp_order, const Index &f_interp_order, const Numeric &p, const Numeric &T, ConstVectorView abs_vmrs, ConstVectorView new_f_grid, const Numeric &extpolfac) const
Extract scalar gas absorption coefficients from the lookup table.
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
bool empty() const
Check if variable is empty.
Definition: matpackIV.cc:52
#define max(a, b)
The Tensor6 class.
Definition: matpackVI.h:1088
void find_nonlinear_continua(ArrayOfIndex &cont, const ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &verbosity)
Find continuum species in abs_species.
A constant view of a Vector.
Definition: matpackI.h:476
void abs_xsec_agendaExecute(Workspace &ws, ArrayOfMatrix &abs_xsec_per_species, ArrayOfMatrix &src_xsec_per_species, ArrayOfArrayOfMatrix &dabs_xsec_per_species_dx, ArrayOfArrayOfMatrix &dsrc_xsec_per_species_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &abs_species_active, const Vector &f_grid, const Vector &abs_p, const Vector &abs_t, const EnergyLevelMap &abs_nlte, const Matrix &abs_vmrs, const Agenda &input_agenda)
Definition: auto_md.cc:23564
Numeric mean(const ConstVectorView &x)
Mean function, vector version.
Definition: matpackI.cc:1589
Index nelem(const Lines &l)
Number of lines.
const Index GFIELD4_P_GRID
const Vector & GetPgrid() const
Workspace class.
Definition: workspace_ng.h:40
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
void resize(Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVI.cc:2175
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1476
ArrayOfGridPosPoly fgp_default
Frequency grid positions.
Structure to store a grid position for higher order interpolation.
#define CREATE_OUT3
Definition: messages.h:207
void seed(unsigned long int n, const Verbosity &verbosity)
Seeds the Rng with the integer argument.
Definition: rng.cc:54
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:23495
Header file for helper functions for OpenMP.
Auxiliary data for isotopologues.
Definition: absorption.h:217
Internal cloudbox functions.
void propmat_clearskyInit(ArrayOfPropagationMatrix &propmat_clearsky, ArrayOfStokesVector &nlte_source, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, ArrayOfStokesVector &dnlte_dx_source, ArrayOfStokesVector &nlte_dsource_dx, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &f_grid, const Index &stokes_dim, const Index &propmat_clearsky_agenda_checked, const Index &nlte_do, const Verbosity &)
WORKSPACE METHOD: propmat_clearskyInit.
Definition: m_abs.cc:1028
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:66
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void parse_atmcompact_speciesname(String &species_name, const String &field_name, const String &delim)
Definition: cloudbox.cc:880
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5484
void abs_lookupTestAccMC(Workspace &ws, const GasAbsLookup &abs_lookup, const Index &abs_lookup_is_adapted, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Index &mc_seed, const Agenda &abs_xsec_agenda, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupTestAccMC.
#define CREATE_OUT2
Definition: messages.h:206
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
void abs_speciesAdd2(Workspace &ws, ArrayOfArrayOfSpeciesTag &abs_species, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, Index &propmat_clearsky_agenda_checked, Index &abs_xsec_agenda_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesAdd2.
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1296
ArrayOfArrayOfSpeciesTag species
The species tags for which the table is valid.
This file contains declerations of functions of physical character.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
void abs_lookupSetup(Vector &abs_p, Vector &abs_t, Vector &abs_t_pert, Matrix &abs_vmrs, ArrayOfArrayOfSpeciesTag &abs_nls, Vector &abs_nls_pert, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &atmfields_checked, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &abs_p_interp_order, const Index &abs_t_interp_order, const Index &abs_nls_interp_order, const Numeric &p_step10, const Numeric &t_step, const Numeric &h2o_step, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupSetup.
const Index GFIELD4_FIELD_NAMES
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1064
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056