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