ARTS 2.5.9 (git: 825fa5f2)
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 <cmath>
28#include <limits>
29#include <map>
30
31#include "absorption.h"
32#include "agenda_class.h"
33#include "arts.h"
34#include "arts_omp.h"
35#include "auto_md.h"
36#include "check_input.h"
37#include "cloudbox.h"
38#include "debug.h"
39#include "gas_abs_lookup.h"
40#include "global_data.h"
41#include "gridded_fields.h"
43#include "math_funcs.h"
44#include "matpackI.h"
45#include "matpackV.h"
46#include "messages.h"
47#include "physics_funcs.h"
48#include "propagationmatrix.h"
49#include "rng.h"
50
53
54/* Workspace method: Doxygen documentation will be auto-generated */
55void abs_lookupInit(GasAbsLookup& x, const Verbosity& verbosity) {
56 ArtsOut2 out2(verbosity);
57 // Nothing to do here.
58 // That means, we rely on the default constructor.
59
60 x = GasAbsLookup();
61 out2 << " Created an empty gas absorption lookup table.\n";
62}
63
64/* Workspace method: Doxygen documentation will be auto-generated */
65void abs_lookupCalc( // Workspace reference:
66 Workspace& ws,
67 // WS Output:
68 GasAbsLookup& abs_lookup,
69 Index& abs_lookup_is_adapted,
70 // WS Input:
71 const ArrayOfArrayOfSpeciesTag& abs_species,
72 const ArrayOfArrayOfSpeciesTag& abs_nls,
73 const Vector& f_grid,
74 const Vector& abs_p,
75 const Matrix& abs_vmrs,
76 const Vector& abs_t,
77 const Vector& abs_t_pert,
78 const Vector& abs_nls_pert,
79 const Agenda& propmat_clearsky_agenda,
80 // GIN
81 const Numeric& lowest_vmr,
82 // Verbosity object:
83 const Verbosity& verbosity) {
86
88 propmat_clearsky_agenda.name() not_eq "propmat_clearsky_agenda",
89 R"--(
90Not propmat_clearsky_agenda:
91
92)--",
93 propmat_clearsky_agenda)
94
95 ARTS_USER_ERROR_IF(lowest_vmr <= 0, R"--(
96You must give a lowest_vmr value above 0. This is because
97the computations of the absorption lookup table uses true
98absorption calculations but the output is in cross-sections
99
100Your current lowest_vmr value is: )--", lowest_vmr)
101
102 // We need this temporary variable to make a local copy of all VMRs,
103 // where we then perturb the H2O profile as needed
104 Matrix these_all_vmrs = abs_vmrs;
105
106 // We will be calling an absorption agenda one species at a
107 // time. This is better than doing all simultaneously, because is
108 // saves memory and allows for consistent treatment of nonlinear
109 // species. But it means we need local copies of species, line list,
110 // and line shapes for agenda communication.
111
112 // 2. Determine various important sizes:
113 const Index n_species = abs_species.nelem(); // Number of abs species
114 const Index n_nls = abs_nls.nelem(); // Number of nonlinear species
115 const Index n_f_grid = f_grid.nelem(); // Number of frequency grid points
116 const Index n_p_grid = abs_p.nelem(); // Number of presure grid points
117 const Index n_t_pert = abs_t_pert.nelem(); // Number of temp. perturbations
118 const Index n_nls_pert = abs_nls_pert.nelem(); // Number of VMR pert. for NLS
119
120 // 3. Input to absorption calculations:
121
122 // Absorption vmrs and temperature:
123 Matrix this_vmr(1, n_p_grid);
124 Vector abs_h2o(n_p_grid);
125 Vector this_t; // Has same dimension, but is
126 // initialized by assignment later.
127 const EnergyLevelMap this_nlte_dummy;
128
129 // Species list, lines, and line shapes, all with only 1 element:
130 ArrayOfArrayOfSpeciesTag this_species(1);
131
132 // List of active species for agenda call. Will always be filled with only
133 // one species.
134 ArrayOfIndex abs_species_active(1);
135
136 // Local copy of nls_pert and t_pert:
137 Vector these_nls_pert; // Is resized every time when used
138 Vector these_t_pert; // Is resized later on
139
140 // 4. Checks of input parameter correctness:
141 const Index h2o_index = find_first_species(abs_species, Species::fromShortName("H2O"));
142
143 if (h2o_index < 0) {
144 // If there are nonlinear species, then at least one species must be
145 // H2O. We will use that to set h2o_abs, and to perturb in the case
146 // of nonlinear species.
147 ARTS_USER_ERROR_IF(n_nls > 0,
148 "If you have nonlinear species, at least one\n"
149 "species must be a H2O species.")
150
151 out2 << " You have no H2O species. Absorption continua will not work.\n"
152 << " You should get a runtime error if you try them anyway.\n";
153 }
154
155 // abs_species, f_grid, and p_grid should not be empty:
156 if (0 == n_species || 0 == n_f_grid || 0 == n_p_grid) {
157 ostringstream os;
158 os << "One of the required input variables is empty:\n"
159 << "abs_species.nelem() = " << n_species << ",\n"
160 << "f_grid.nelem() = " << n_f_grid << ",\n"
161 << "abs_p.nelem() = " << n_p_grid << ".";
162 throw runtime_error(os.str());
163 }
164
165 // Set up the index array abs_nls from the tag array
166 // abs_nls. Give an error message if these
167 // tags are not included in abs_species.
168 ArrayOfIndex abs_nls_idx;
169 for (Index i = 0; i < n_nls; ++i) {
170 Index s;
171 for (s = 0; s < n_species; ++s) {
172 if (abs_nls[i] == abs_species[s]) {
173 abs_nls_idx.push_back(s);
174 break;
175 }
176 }
177 if (s == n_species) {
178 ostringstream os;
179 os << "Did not find *abs_nls* tag group \""
180 << abs_nls[i] << "\" in *abs_species*.";
181 throw runtime_error(os.str());
182 }
183 }
184
185 // Furthermore abs_nls_idx must not contain duplicate values:
186 if (!is_unique(abs_nls_idx)) {
187 ostringstream os;
188 os << "The variable *abs_nls* must not have duplicate species.\n"
189 << "Value of *abs_nls*: " << abs_nls_idx;
190 throw runtime_error(os.str());
191 }
192
193 // VMR matrix must match species list and pressure levels:
194 chk_size("abs_vmrs", abs_vmrs, n_species, n_p_grid);
195
196 // Temperature vector must match number of pressure levels:
197 chk_size("abs_t", abs_t, n_p_grid);
198
199 // abs_nls_pert should only be not-empty if we have nonlinear species:
200 if (((0 == n_nls) && (0 != n_nls_pert)) ||
201 ((0 != n_nls) && (0 == n_nls_pert))) {
202 ostringstream os;
203 os << "You have to set both abs_nls and abs_nls_pert, or none.";
204 throw runtime_error(os.str());
205 }
206
207 // 4.a Set up a logical array for the nonlinear species.
208 ArrayOfIndex non_linear(n_species, 0);
209 for (Index s = 0; s < n_nls; ++s) {
210 non_linear[abs_nls_idx[s]] = 1;
211 }
212
213 // 5. Set general lookup table properties:
214 abs_lookup.species = abs_species; // Species list
215 abs_lookup.nonlinear_species =
216 abs_nls_idx; // Nonlinear species (e.g., H2O, O2)
217 abs_lookup.f_grid = f_grid; // Frequency grid
218 abs_lookup.p_grid = abs_p; // Pressure grid
219 abs_lookup.vmrs_ref = abs_vmrs;
220 abs_lookup.t_ref = abs_t;
221 abs_lookup.t_pert = abs_t_pert;
222 abs_lookup.nls_pert = abs_nls_pert;
223
224 // 5.a. Set log_p_grid:
225 abs_lookup.log_p_grid.resize(n_p_grid);
226 transform(abs_lookup.log_p_grid, log, abs_lookup.p_grid);
227
228 // 6. Create abs_lookup.xsec with the right dimensions:
229 {
230 Index a, b, c, d;
231
232 if (0 == n_t_pert)
233 a = 1;
234 else
235 a = n_t_pert;
236
237 b = n_species + n_nls * (n_nls_pert - 1);
238
239 c = n_f_grid;
240
241 d = n_p_grid;
242
243 abs_lookup.xsec.resize(a, b, c, d);
244 abs_lookup.xsec = NAN;
245 }
246
247 // 6.a. Set up these_t_pert. This is done so that we can use the
248 // same loop over the perturbations, independent of
249 // whether we have temperature perturbations or not.
250 if (0 != n_t_pert) {
251 out2 << " With temperature perturbations.\n";
252 these_t_pert.resize(n_t_pert);
253 these_t_pert = abs_t_pert;
254 } else {
255 out2 << " No temperature perturbations.\n";
256 these_t_pert.resize(1);
257 these_t_pert = 0;
258 }
259
260 const Index these_t_pert_nelem = these_t_pert.nelem();
261
262 // 7. Now we have to fill abs_lookup.xsec with the right values!
263
264 String fail_msg;
265 bool failed = false;
266
267 // Loop species:
268 for (Index i = 0, spec = 0; i < n_species; ++i) {
269 // Skipping Zeeman and free_electrons species.
270 // (Mixed tag groups between those and other species are not allowed.)
271 if (abs_species[i].Zeeman() || abs_species[i].FreeElectrons() || abs_species[i].Particles()) {
272 spec++;
273 continue;
274 }
275
276 // spec is the index for the second dimension of abs_lookup.xsec.
277
278 // Prepare absorption agenda input for this species:
279 out2 << " Doing species " << i + 1 << " of " << n_species << ": "
280 << abs_species[i] << ".\n";
281
282 // Set active species:
283 abs_species_active[0] = i;
284
285 // Get a dummy list of tag groups with only the current element:
286 this_species[0].resize(abs_species[i].nelem());
287 this_species[0] = abs_species[i];
288
289 // Set up these_nls_pert. This is done so that we can use the
290 // same loop over the perturbations, independent of
291 // whether we have nonlinear species or not.
292 if (non_linear[i]) {
293 out2 << " This is a species with H2O VMR perturbations.\n";
294 these_nls_pert.resize(n_nls_pert);
295 these_nls_pert = abs_nls_pert;
296 } else {
297 these_nls_pert.resize(1);
298 these_nls_pert = 1;
299 }
300
301 // Loop these_nls_pert:
302 for (Index s = 0; s < these_nls_pert.nelem(); ++s, ++spec) {
303 // Remember, spec is the index for the second dimension of
304 // abs_lookup.xsec
305
306 if (non_linear[i]) {
307 out2 << " Doing H2O VMR variant " << s + 1 << " of " << n_nls_pert
308 << ": " << abs_nls_pert[s] << ".\n";
309 }
310
311 // Make a local copy of the VMRs, and manipulate the H2O VMR within it.
312 // Note: We do not need a runtime error check that h2o_index is ok here,
313 // because earlier on we throw an error if there is no H2O species although we
314 // need it. So, if h2o_indes is -1, we here simply assume that there
315 // should not be a perturbation
316 if (h2o_index >= 0) {
317 these_all_vmrs(h2o_index, joker) = abs_vmrs(h2o_index, joker);
318 these_all_vmrs(h2o_index, joker) *=
319 these_nls_pert[s]; // Add perturbation
320 }
321
322 // VMR for this species (still needed by interfact to continua):
323 // FIXME: This variable may go away eventually, when the continuum
324 // part no longer needs it.
325 this_vmr(0, joker) = these_all_vmrs(i, joker);
326
327 // For abs_h2o, we can always add the perturbations (it will
328 // not make a difference if the species itself is also H2O).
329 // Attention, we need to treat here also the case that there
330 // is no H2O species. We will then set abs_h2o to
331 // -1. Absorption routines that do not really need abs_h2o
332 // will still run.
333 //
334 // FIXME: abs_h2o is currently still needed by the continuum part.
335 // Should go away eventually.
336 if (h2o_index == -1) {
337 // The case without H2O species.
338 abs_h2o.resize(1);
339 abs_h2o = -1;
340 } else {
341 // The normal case.
342 abs_h2o = these_all_vmrs(h2o_index, joker);
343 }
344
345 // Loop temperature perturbations
346 // ------------------------------
347
348 // We use a parallel for loop for this.
349
350 // There is something strange here: abs_lookup seems to be
351 // "shared" by default, although I have set default(none). I
352 // suspect that the reason for this behavior is that
353 // abs_lookup is a return by reference parameter of this
354 // function. Anyway, shared is the correct setting for
355 // abs_lookup, so there is no problem.
356
358
359#pragma omp parallel for if ( \
360 !arts_omp_in_parallel() && \
361 these_t_pert_nelem >= \
362 arts_omp_get_max_threads()) private(this_t) \
363 firstprivate(wss)
364 for (Index j = 0; j < these_t_pert_nelem; ++j) {
365 // Skip remaining iterations if an error occurred
366 if (failed) continue;
367
368 // The try block here is necessary to correctly handle
369 // exceptions inside the parallel region.
370 try {
371 if (0 != n_t_pert) {
372 // We first prepare the output in a string here,
373 // so that we can write it to out3 with a single
374 // operation. This avoids messy output from
375 // multiple threads.
376 ostringstream os;
377
378 os << " Doing temperature variant " << j + 1 << " of " << n_t_pert
379 << ": " << these_t_pert[j] << ".\n";
380
381 out3 << os.str();
382 }
383
384 // Create perturbed temperature profile:
385 this_t = abs_lookup.t_ref;
386 this_t += these_t_pert[j];
387
388 // Store in the right place:
389 // Loop through all altitudes
390 for (Index p = 0; p < n_p_grid; ++p) {
392 StokesVector S;
395 Vector rtp_vmr = these_all_vmrs(joker, p);
396 for (auto& x : rtp_vmr) x = std::max(lowest_vmr, x);
397
398 // Perform the propagation matrix computations
400 K,
401 S,
402 dK,
403 dS,
404 {},
405 abs_species[i],
406 f_grid,
407 {},
408 {},
409 abs_p[p],
410 this_t[p],
411 {},
412 rtp_vmr,
413 propmat_clearsky_agenda);
414 K.Kjj() /= rtp_vmr[i] * number_density(abs_p[p], this_t[p]);
415 abs_lookup.xsec(j, spec, Range(joker), p) = K.Kjj();
416
417 // There used to be a division by the number density
418 // n here. This is no longer necessary, since
419 // abs_xsec_per_species now contains true absorption
420 // cross sections.
421
422 // IMPORTANT: There was a bug in my old Matlab
423 // function "create_lookup.m" to generate the lookup
424 // table. (The number density was always the
425 // reference one, and did not change for different
426 // temperatures.) Patricks Atmlab function
427 // "arts_abstable_from_arts1.m" did *not* have this bug.
428
429 // Calculate the number density for the given pressure and
430 // temperature:
431 // n = n0*T0/p0 * p/T or n = p/kB/t, ideal gas law
432 // const Numeric n = number_density( abs_lookup.p_grid[p],
433 // this_t[p] );
434 // abs_lookup.xsec( j, spec, Range(joker), p ) /= n;
435 }
436 } // end of try block
437 catch (const std::runtime_error& e) {
438#pragma omp critical(abs_lookupCalc_fail)
439 {
440 fail_msg = e.what();
441 failed = true;
442 }
443 }
444 } // end of parallel for loop
445
446 if (failed) throw runtime_error(fail_msg);
447 }
448 }
449
450 // 6. Initialize fgp_default.
451 abs_lookup.flag_default = Interpolation::LagrangeVector(abs_lookup.f_grid, abs_lookup.f_grid, 0);
452
453 // Set the abs_lookup_is_adapted flag. After all, the table fits the
454 // current frequency grid and species selection.
455 abs_lookup_is_adapted = 1;
456}
457
459
477 const ArrayOfArrayOfSpeciesTag& abs_species,
478 const Verbosity& verbosity) {
480
481 cont.resize(0);
482
483 // This is quite complicated, unfortunately. The approach here
484 // is copied from abs_xsec_per_speciesAddConts. For explanation,
485 // see there.
486
487 // Loop tag groups:
488 for (Index i = 0; i < abs_species.nelem(); ++i) {
489 // Loop tags in tag group
490 for (Index s = 0; s < abs_species[i].nelem(); ++s) {
491 // Check for continuum tags
492 if (abs_species[i][s].type == Species::TagType::Predefined ||
493 abs_species[i][s].type == Species::TagType::Cia) {
494 const String thisname = abs_species[i][s].Name();
495 // Ok, now we know this is a continuum tag.
496 out3 << " Continuum tag: " << thisname;
497
498 // Check whether we want nonlinear treatment for
499 // this or not. We have three classes of continua:
500 // 1. Those that we know do not require it
501 // 2. Those that we know require h2o_abs
502 // 3. Those for which we do not know
503
504 // The list here is not at all perfect, just a first
505 // guess. If you need particular models, you better
506 // check that the right thing is done with your model.
507
508 // 1. Continua known to not use h2o_abs
509 // We take also H2O itself here, since this is
510 // handled separately
511 if (Species::fromShortName("H2O") == abs_species[i][s].Spec() ||
512 "N2-" == thisname.substr(0, 3) || "CO2-" == thisname.substr(0, 4) ||
513 "O2-CIA" == thisname.substr(0, 6) ||
514 "O2-v0v" == thisname.substr(0, 6) ||
515 "O2-v1v" == thisname.substr(0, 6) ||
516 "H2-CIA" == thisname.substr(0, 6) ||
517 "He-CIA" == thisname.substr(0, 6) ||
518 "CH4-CIA" == thisname.substr(0, 7) ||
519 "liquidcloud-MPM93" == thisname.substr(0, 17) ||
520 "liquidcloud-ELL07" == thisname.substr(0, 17)) {
521 out3 << " --> not added.\n";
522 break;
523 }
524
525 // 2. Continua known to use h2o_abs
526 if ("O2-" == thisname.substr(0, 3)) {
527 cont.push_back(i);
528 out3 << " --> added to abs_nls.\n";
529 break;
530 }
531
532 // 3. abs_species tags that are NOT allowed in LUT
533 // calculations
534 if ("icecloud-" == thisname.substr(0, 9) ||
535 "rain-" == thisname.substr(0, 5)) {
536 ostringstream os;
537 os << "Tag " << thisname << " not allowed in absorption "
538 << "lookup tables.";
539 throw runtime_error(os.str());
540 }
541
542 // If we get here, then the tag was neither in the
543 // posivitive nor in the negative list. We through a
544 // runtime error.
545 out3 << " --> unknown.\n";
546 ostringstream os;
547 os << "Unknown whether tag " << thisname
548 << " is a nonlinear species (i.e. uses h2o_abs) or not.\n"
549 << "Cannot set abs_nls automatically.";
550 throw runtime_error(os.str());
551 }
552 }
553 }
554}
555
557
566 const ArrayOfArrayOfSpeciesTag& abs_species,
567 const Verbosity& verbosity) {
569
570 abs_nls.resize(0);
571
572 // Add all H2O species as non-linear:
573 Index next_h2o = 0;
574 while (-1 !=
575 (next_h2o = find_next_species(
576 abs_species, Species::fromShortName("H2O"), next_h2o))) {
577 abs_nls.push_back(abs_species[next_h2o]);
578 ++next_h2o;
579 }
580
581 // Certain continuum models also depend on abs_h2o. There is a
582 // helper function that contains a list of these.
583 ArrayOfIndex cont;
584 find_nonlinear_continua(cont, abs_species, verbosity);
585
586 // Add these to abs_nls:
587 for (Index i = 0; i < cont.nelem(); ++i) {
588 abs_nls.push_back(abs_species[cont[i]]);
589 }
590
591 out2 << " Species marked for nonlinear treatment:\n";
592 for (Index i = 0; i < abs_nls.nelem(); ++i) {
593 out2 << " ";
594 for (Index j = 0; j < abs_nls[i].nelem(); ++j) {
595 if (j != 0) out2 << ", ";
596 out2 << abs_nls[i][j].Name();
597 }
598 out2 << "\n";
599 }
600}
601
603
616void choose_abs_t_pert(Vector& abs_t_pert,
617 ConstVectorView abs_t,
618 ConstVectorView tmin,
619 ConstVectorView tmax,
620 const Numeric& step,
621 const Index& p_interp_order,
622 const Index& t_interp_order,
623 const Verbosity& verbosity) {
626
627 // The code to find out the range for perturbation is a bit
628 // complicated. The problem is that, since we use higher order
629 // interpolation for p, we may require temperatures well outside the
630 // local min/max range at the reference level. We solve this by
631 // really looking at the min/max range for those levels required by
632 // the p_interp_order.
633
634 Numeric mindev = 1e9;
635 Numeric maxdev = -1e9;
636
637 Vector the_grid(0, abs_t.nelem(), 1);
638 for (Index i = 0; i < the_grid.nelem(); ++i) {
639 const Index idx0 = Interpolation::pos_finder(i, Numeric(i), the_grid, p_interp_order, false, true);
640
641 for (Index j = 0; j < p_interp_order+1; ++j) {
642 // Our pressure grid for the lookup table may be coarser than
643 // the original one for the batch cases. This may lead to max/min
644 // values in the original data that exceed those we assumed
645 // above. We add some extra margin here to account for
646 // that. (The margin is +-10K)
647
648 Numeric delta_min = tmin[i] - abs_t[idx0 + j] - 10;
649 Numeric delta_max = tmax[i] - abs_t[idx0 + j] + 10;
650
651 if (delta_min < mindev) mindev = delta_min;
652 if (delta_max > maxdev) maxdev = delta_max;
653 }
654 }
655
656 out3 << " abs_t_pert: mindev/maxdev : " << mindev << " / " << maxdev << "\n";
657
658 // We divide the interval between mindev and maxdev, so that the
659 // steps are of size *step* or smaller. But we also need at least
660 // *t_interp_order*+1 points.
661 Index div = t_interp_order;
662 Numeric effective_step;
663 do {
664 effective_step = (maxdev - mindev) / (Numeric)div;
665 ++div;
666 } while (effective_step > step);
667
668 abs_t_pert = Vector(mindev, div, effective_step);
669
670 out2 << " abs_t_pert: " << abs_t_pert[0] << " K to "
671 << abs_t_pert[abs_t_pert.nelem() - 1] << " K in steps of "
672 << effective_step << " K (" << abs_t_pert.nelem() << " grid points)\n";
673}
674
676
689void choose_abs_nls_pert(Vector& abs_nls_pert,
690 ConstVectorView refprof,
691 ConstVectorView minprof,
692 ConstVectorView maxprof,
693 const Numeric& step,
694 const Index& p_interp_order,
695 const Index& nls_interp_order,
696 const Verbosity& verbosity) {
699
700 // The code to find out the range for perturbation is a bit
701 // complicated. The problem is that, since we use higher order
702 // interpolation for p, we may require humidities well outside the
703 // local min/max range at the reference level. We solve this by
704 // really looking at the min/max range for those levels required by
705 // the p_interp_order.
706
707 Numeric mindev = 0;
708 Numeric maxdev = -1e9;
709
710 // mindev is set to zero from the start, since we always want to
711 // include 0.
712
713 Vector the_grid(0, refprof.nelem(), 1);
714 for (Index i = 0; i < the_grid.nelem(); ++i) {
715 // cout << "!!!!!! i = " << i << "\n";
716 // cout << " min/ref/max = " << minprof[i] << " / "
717 // << refprof[i] << " / "
718 // << maxprof[i] << "\n";
719
720 const Index idx0 = Interpolation::pos_finder(i, Numeric(i), the_grid, p_interp_order, false, true);
721
722 for (Index j = 0; j < p_interp_order+1; ++j) {
723 // cout << "!!!!!! j = " << j << "\n";
724 // cout << " ref[j] = " << refprof[gp.idx[j]] << " ";
725 // cout << " minfrac[j] = " << minprof[i] / refprof[gp.idx[j]] << " ";
726 // cout << " maxfrac[j] = " << maxprof[i] / refprof[gp.idx[j]] << " \n";
727
728 // Our pressure grid for the lookup table may be coarser than
729 // the original one for the batch cases. This may lead to max/min
730 // values in the original data that exceed those we assumed
731 // above. We add some extra margin to the max value here to account for
732 // that. (The margin is a factor of 2.)
733
734 Numeric delta_min = minprof[i] / refprof[idx0 + j];
735 Numeric delta_max = 2 * maxprof[i] / refprof[idx0 + j];
736
737 if (delta_min < mindev) mindev = delta_min;
738 // do not update maxdev, when delta_max is infinity (this results from
739 // refprof being 0)
740 if (!std::isinf(delta_max) && (delta_max > maxdev)) maxdev = delta_max;
741 }
742 }
743
744 out3 << " abs_nls_pert: mindev/maxdev : " << mindev << " / " << maxdev
745 << "\n";
746
747 bool allownegative = false;
748 if (mindev < 0) {
749 out2
750 << " Warning: I am getting a negative fractional distance to the H2O\n"
751 << " reference profile. Some of your H2O fields may contain negative values.\n"
752 << " Will allow negative values also for abs_nls_pert.\n";
753 allownegative = true;
754 }
755
756 if (!allownegative) {
757 mindev = 0;
758 out3 << " Adjusted mindev : " << mindev << "\n";
759 }
760
761 if (std::isinf(maxdev)) {
762 ostringstream os;
763 os << "Perturbation upper limit is infinity (likely due to the reference\n"
764 << "profile being 0 at at least one pressure level). Can not work\n"
765 << "with that.";
766 throw runtime_error(os.str());
767 }
768
769 // We divide the interval between mindev and maxdev, so that the
770 // steps are of size *step* or smaller. But we also need at least
771 // *nls_interp_order*+1 points.
772 Index div = nls_interp_order;
773 Numeric effective_step;
774 do {
775 effective_step = (maxdev - mindev) / (Numeric)div;
776 ++div;
777 } while (effective_step > step);
778
779 abs_nls_pert = Vector(mindev, div, effective_step);
780
781 // If there are negative values, we also add 0. The reason for this
782 // is that 0 is a turning point.
783 if (allownegative) {
784 VectorInsertGridPoints(abs_nls_pert, abs_nls_pert, {0}, verbosity);
785 out2
786 << " I am including also 0 in the abs_nls_pert, because it is a turning \n"
787 << " point. Consider to use a higher abs_nls_interp_order, for example 4.\n";
788 }
789
790 out2 << " abs_nls_pert: " << abs_nls_pert[0] << " to "
791 << abs_nls_pert[abs_nls_pert.nelem() - 1]
792 << " (fractional units) in steps of "
793 << abs_nls_pert[1] - abs_nls_pert[0] << " (" << abs_nls_pert.nelem()
794 << " grid points)\n";
796
797/* Workspace method: Doxygen documentation will be auto-generated */
798void abs_lookupSetup( // WS Output:
799 Vector& abs_p,
800 Vector& abs_t,
801 Vector& abs_t_pert,
802 Matrix& abs_vmrs,
804 Vector& abs_nls_pert,
805 // WS Input:
806 const Index& atmosphere_dim,
807 const Vector& p_grid,
808 // const Vector& lat_grid,
809 // const Vector& lon_grid,
810 const Tensor3& t_field,
811 const Tensor4& vmr_field,
812 const Index& atmfields_checked,
813 const ArrayOfArrayOfSpeciesTag& abs_species,
814 const Index& abs_p_interp_order,
815 const Index& abs_t_interp_order,
816 const Index& abs_nls_interp_order,
817 // Control Parameters:
818 const Numeric& p_step10,
819 const Numeric& t_step,
820 const Numeric& h2o_step,
821 const Verbosity& verbosity) {
822 // Checks on input parameters:
823
824 if (atmfields_checked != 1)
825 throw runtime_error(
826 "The atmospheric fields must be flagged to have "
827 "passed a consistency check (atmfields_checked=1).");
828
829 // We don't actually need lat_grid and lon_grid, but we have them as
830 // input variables, so that we can use the standard functions to
831 // check atmospheric fields and grids. A bit cheesy, but I don't
832 // want to program all the checks explicitly.
833
834 // Check grids (outcommented the ones that have been done by
835 // atmfields_checkedCalc already):
836 //chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
837
838 if (p_grid.nelem() < 2) {
839 ostringstream os;
840 os << "You need at least two pressure levels.";
841 throw runtime_error(os.str());
842 }
843
844 // Check T field:
845 //chk_atm_field("t_field", t_field, atmosphere_dim,
846 // p_grid, lat_grid, lon_grid);
847
848 // Check VMR field (and abs_species):
849 //chk_atm_field("vmr_field", vmr_field, atmosphere_dim,
850 // abs_species.nelem(), p_grid, lat_grid, lon_grid);
851
852 // Check the keyword arguments:
853 if (p_step10 <= 0 || t_step <= 0 || h2o_step <= 0) {
854 ostringstream os;
855 os << "The keyword arguments p_step, t_step, and h2o_step must be >0.";
856 throw runtime_error(os.str());
857 }
858
859 // Ok, all input parameters seem to be reasonable.
860
861 // For consistency with other code around arts (e.g., correlation
862 // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
863 // convert it here to the natural log:
864 const Numeric p_step = log(pow(10.0, p_step10));
865
866 // We will need the log of the pressure grid:
867 Vector log_p_grid(p_grid.nelem());
868 transform(log_p_grid, log, p_grid);
869
870 // const Numeric epsilon = 0.01 * p_step; // This is the epsilon that
871 // // we use for comparing p grid spacings.
872
873 // Construct abs_p
874 // ---------------
875
876 ArrayOfNumeric log_abs_p_a; // We take log_abs_p_a as an array of
877 // Numeric, so that we can easily
878 // build it up by appending new elements to the end.
879
880 // Check whether there are pressure levels that are further apart
881 // (in log(p)) than p_step, and insert additional levels if
882 // necessary:
883
884 log_abs_p_a.push_back(log_p_grid[0]);
885
886 for (Index i = 1; i < log_p_grid.nelem(); ++i) {
887 const Numeric dp =
888 log_p_grid[i - 1] - log_p_grid[i]; // The grid is descending.
889
890 const Numeric dp_by_p_step = dp / p_step;
891 // cout << "dp_by_p_step: " << dp_by_p_step << "\n";
892
893 // How many times does p_step fit into dp?
894 const Index n = (Index)ceil(dp_by_p_step);
895 // n is the number of intervals that we want to have in the
896 // new grid. The number of additional points to insert is
897 // n-1. But we have to insert the original point as well.
898 // cout << n << "\n";
899
900 const Numeric ddp = dp / (Numeric)n;
901 // cout << "ddp: " << ddp << "\n";
902
903 for (Index j = 1; j <= n; ++j)
904 log_abs_p_a.push_back(log_p_grid[i - 1] - (Numeric)j * ddp);
905 }
906
907 // Copy to a proper vector, we need this also later for
908 // interpolation:
909 Vector log_abs_p(log_abs_p_a.nelem());
910 for (Index i = 0; i < log_abs_p_a.nelem(); ++i) log_abs_p[i] = log_abs_p_a[i];
911
912 // Copy the new grid to abs_p, removing the log:
913 abs_p.resize(log_abs_p.nelem());
914 transform(abs_p, exp, log_abs_p);
915
916 // Check that abs_p has enough points for the interpolation order
917 if (abs_p.nelem() < abs_p_interp_order + 1) {
918 ostringstream os;
919 os << "Your pressure grid does not have enough levels for the desired interpolation order.";
920 throw runtime_error(os.str());
921 }
922
923 // We will also have to interpolate T and VMR profiles to the new
924 // pressure grid. We interpolate in log(p), as usual in ARTS.
925
926 // Grid positions:
927 ArrayOfGridPos gp(log_abs_p.nelem());
928 gridpos(gp, log_p_grid, log_abs_p);
929
930 // Interpolation weights:
931 Matrix itw(gp.nelem(), 2);
932 interpweights(itw, gp);
933
934 // In the 1D case the lookup table is just a lookup table in
935 // pressure. We treat this simple case first.
936 if (1 == atmosphere_dim) {
937 // Reference temperature,
938 // interpolate abs_t from t_field:
939 abs_t.resize(log_abs_p.nelem());
940 interp(abs_t, itw, t_field(joker, 0, 0), gp);
941
942 // Temperature perturbations:
943 abs_t_pert.resize(0);
944
945 // Reference VMR profiles,
946 // interpolate abs_vmrs from vmr_field:
947 abs_vmrs.resize(abs_species.nelem(), log_abs_p.nelem());
948 for (Index i = 0; i < abs_species.nelem(); ++i)
949 interp(abs_vmrs(i, joker), itw, vmr_field(i, joker, 0, 0), gp);
950
951 // Species for which H2O VMR is perturbed:
952 abs_nls.resize(0);
953
954 // H2O VMR perturbations:
955 abs_nls_pert.resize(0);
956 } else {
957 // 2D or 3D case. We have to set up T and nonlinear species variations.
958
959 // Make an intelligent choice for the nonlinear species.
960 choose_abs_nls(abs_nls, abs_species, verbosity);
961
962 // Now comes a part where we analyse the atmospheric fields.
963 // We need to find the max, min, and mean profile for
964 // temperature and VMRs.
965 // We do this on the original p grid, not on the refined p
966 // grid, to be more efficient.
967
968 // Temperature:
969 Vector tmin(p_grid.nelem());
970 Vector tmax(p_grid.nelem());
971 Vector tmean(p_grid.nelem());
972
973 for (Index i = 0; i < p_grid.nelem(); ++i) {
974 tmin[i] = min(t_field(i, joker, joker));
975 tmax[i] = max(t_field(i, joker, joker));
976 tmean[i] = mean(t_field(i, joker, joker));
977 }
978
979 // cout << "Tmin: " << tmin << "\n";
980 // cout << "Tmax: " << tmax << "\n";
981 // cout << "Tmean: " << tmean << "\n";
982
983 // Calculate mean profiles of all species. (We need all for abs_vmrs
984 // not only H2O.)
985 Matrix vmrmean(abs_species.nelem(), p_grid.nelem());
986 for (Index s = 0; s < abs_species.nelem(); ++s)
987 for (Index i = 0; i < p_grid.nelem(); ++i) {
988 vmrmean(s, i) = mean(vmr_field(s, i, joker, joker));
989 }
990
991 // If there are NLS, determine H2O statistics:
992
993 // We have to define these here, outside the if block, because
994 // we need the values later.
995 Vector h2omin(p_grid.nelem());
996 Vector h2omax(p_grid.nelem());
997 const Index h2o_index = find_first_species(abs_species, Species::fromShortName("H2O"));
998 // We need this inside the if clauses for nonlinear species
999 // treatment. The function returns "-1" if there is no H2O
1000 // species. There is a check for that in the next if block, with
1001 // an appropriate runtime error.
1002
1003 if (0 < abs_nls.nelem()) {
1004 // Check if there really is a H2O species.
1005 if (h2o_index < 0) {
1006 ostringstream os;
1007 os << "Some of your species require nonlinear treatment,\n"
1008 << "but you have no H2O species.";
1009 throw runtime_error(os.str());
1010 }
1011
1012 for (Index i = 0; i < p_grid.nelem(); ++i) {
1013 h2omin[i] = min(vmr_field(h2o_index, i, joker, joker));
1014 h2omax[i] = max(vmr_field(h2o_index, i, joker, joker));
1015 }
1016
1017 // cout << "H2Omin: " << h2omin << "\n";
1018 // cout << "H2Omax: " << h2omax << "\n";
1019 // cout << "H2Omean: " << vmrmean(h2o_index,joker) << "\n";
1020 }
1021
1022 // Interpolate in pressure, set abs_t, abs_vmr...
1023
1024 // Reference temperature,
1025 // interpolate abs_t from tmean:
1026 abs_t.resize(log_abs_p.nelem());
1027 interp(abs_t, itw, tmean, gp);
1028
1029 // Temperature perturbations:
1030 choose_abs_t_pert(abs_t_pert,
1031 tmean,
1032 tmin,
1033 tmax,
1034 t_step,
1035 abs_p_interp_order,
1036 abs_t_interp_order,
1037 verbosity);
1038 // cout << "abs_t_pert: " << abs_t_pert << "\n";
1039
1040 // Reference VMR profiles,
1041 // interpolate abs_vmrs from vmrmean:
1042 abs_vmrs.resize(abs_species.nelem(), log_abs_p.nelem());
1043 for (Index i = 0; i < abs_species.nelem(); ++i)
1044 interp(abs_vmrs(i, joker), itw, vmrmean(i, joker), gp);
1045
1046 if (0 < abs_nls.nelem()) {
1047 // Construct abs_nls_pert:
1048 choose_abs_nls_pert(abs_nls_pert,
1049 vmrmean(h2o_index, joker),
1050 h2omin,
1051 h2omax,
1052 h2o_step,
1053 abs_p_interp_order,
1054 abs_nls_interp_order,
1055 verbosity);
1056 } else {
1057 // Empty abs_nls_pert:
1058 abs_nls_pert.resize(0);
1059 }
1060 // cout << "abs_nls_pert: " << abs_nls_pert << "\n";
1061 }
1063
1064/* Workspace method: Doxygen documentation will be auto-generated */
1065void abs_lookupSetupBatch( // WS Output:
1066 Vector& abs_p,
1067 Vector& abs_t,
1068 Vector& abs_t_pert,
1069 Matrix& abs_vmrs,
1070 ArrayOfArrayOfSpeciesTag& abs_nls,
1071 Vector& abs_nls_pert,
1072 // WS Input:
1073 const ArrayOfArrayOfSpeciesTag& abs_species,
1074 const ArrayOfGriddedField4& batch_fields,
1075 const Index& abs_p_interp_order,
1076 const Index& abs_t_interp_order,
1077 const Index& abs_nls_interp_order,
1078 const Index& atmosphere_dim,
1079 // Control Parameters:
1080 const Numeric& p_step10,
1081 const Numeric& t_step,
1082 const Numeric& h2o_step,
1083 const Vector& extremes,
1084 const Index& robust,
1085 const Index& check_gridnames,
1086 const Verbosity& verbosity) {
1090
1091 // For consistency with other code around arts (e.g., correlation
1092 // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
1093 // convert it here to the natural log:
1094 const Numeric p_step = log(pow(10.0, p_step10));
1095
1096 // Derive which abs_species is H2O (required for nonlinear species handling)
1097 // returns -1 if no H2O present
1098 const Index h2o_index = find_first_species(abs_species, Species::fromShortName("H2O"));
1099 // cout << "The " << h2o_index+1 << ". species in abs_species is H2O\n";
1100 // cout << "That is, H2O is expected to be the " << indoff+h2o_index
1101 // << ". column of the atmospheric fields\n";
1102
1103 ArrayOfIndex batch_index(abs_species.nelem());
1104 Index T_index = -1;
1105 Index z_index = -1;
1106
1107 ArrayOfString species_names(abs_species.nelem());
1108 for (Index i = 0; i < abs_species.nelem(); ++i)
1109 species_names[i] = Species::toShortName(abs_species[i].Species());
1110
1111 const ArrayOfString field_names = batch_fields[0].get_string_grid(0);
1112
1113 String species_type;
1114 String species_name;
1115 const String delim = "-";
1116
1117 // Check that the field names in batch_fields are good. (We check
1118 // only the first field in the batch.)
1119 {
1120 ostringstream os;
1121 bool bail = false;
1122
1123 // One of the abs_species has to be H2O
1124 // (ideally that should be checked later, when we know whether there are
1125 // nonlinear species present. however, currently batch mode REQUIRES H2O to
1126 // be present, so we check that immediately)
1127 if (h2o_index < 0) {
1128 os << "One of the atmospheric fields must be H2O.\n";
1129 bail = true;
1130 }
1131
1132 // First simply check dimensions of abs_species and field_names
1133 if (field_names.nelem() < 3) {
1134 os << "Atmospheric states in batch_fields must have at\n"
1135 << "least three fields: T, z, and at least one absorption species.";
1136 throw runtime_error(os.str());
1137 }
1138
1139 if (abs_species.nelem() < 1) {
1140 os << "At least one absorption species needs to be defined "
1141 << "in abs_species.";
1142 throw runtime_error(os.str());
1143 }
1144
1145 // Check that all required fields are present.
1146 const Index nf = field_names.nelem();
1147 bool found;
1148
1149 // Looking for temperature field:
1150 found = false;
1151 for (Index fi = 0; fi < nf; ++fi) {
1152 parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1153 if (species_type == "T") {
1154 if (found) {
1155 os << "Only one temperature ('T') field allowed, "
1156 << "but found at least 2.\n";
1157 bail = true;
1158 } else {
1159 found = true;
1160 T_index = fi;
1161 }
1162 }
1163 }
1164 if (!found) {
1165 os << "One temperature ('T') field required, but none found.\n";
1166 bail = true;
1167 }
1168
1169 // Looking for altitude field:
1170 found = false;
1171 for (Index fi = 0; fi < nf; ++fi) {
1172 parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1173 if (species_type == "z") {
1174 if (found) {
1175 os << "Only one altitude ('z') field allowed, "
1176 << "but found at least 2.\n";
1177 bail = true;
1178 } else {
1179 found = true;
1180 z_index = fi;
1181 }
1182 }
1183 }
1184 if (!found) {
1185 os << "One altitude ('z') field required, but none found.\n";
1186 bail = true;
1187 }
1188
1189 // Now going over all abs_species elements and match the corresponding
1190 // field_name (we always take the first occurence of a matching field). We
1191 // don't care that batch_fields contains further fields, too. And we can't
1192 // expect the fields being in the same order as the abs_species.
1193 // At the same time, we keep the indices of the fields corresponding to the
1194 // abs_species elements, because later we will need the field data.
1195 Index fi;
1196 for (Index j = 0; j < abs_species.nelem(); ++j) {
1197 found = false;
1198 fi = 0;
1199 while (!found && fi < nf) {
1200 parse_atmcompact_speciestype(species_type, field_names[fi], delim);
1201 // do we have an abs_species type field?
1202 if (species_type == "abs_species") {
1203 parse_atmcompact_speciesname(species_name, field_names[fi], delim);
1204 if (species_name == species_names[j]) {
1205 found = true;
1206 batch_index[j] = fi;
1207 }
1208 }
1209 fi++;
1210 }
1211 if (!found) {
1212 os << "No field for absorption species '" << species_names[j]
1213 << "'found.\n";
1214 bail = true;
1215 }
1216 }
1217
1218 os << "Your field names are:\n" << field_names;
1219
1220 if (bail) throw runtime_error(os.str());
1221 }
1222
1223 // FIXME: Adjustment of min/max values for Jacobian perturbations is still missing.
1224
1225 // Make an intelligent choice for the nonlinear species.
1226 choose_abs_nls(abs_nls, abs_species, verbosity);
1227
1228 // Find out maximum and minimum pressure and check that pressure grid is decreasing.
1229 Numeric maxp = batch_fields[0].get_numeric_grid(GFIELD4_P_GRID)[0];
1230 Numeric minp = batch_fields[0].get_numeric_grid(
1231 GFIELD4_P_GRID)[batch_fields[0].get_numeric_grid(GFIELD4_P_GRID).nelem() -
1232 1];
1233
1234 ArrayOfIndex valid_field_indices;
1235 for (Index i = 0; i < batch_fields.nelem(); ++i) {
1236 // Local variables for atmfields_check.
1237 Index atmfields_checked;
1238 Tensor4 t4_dummy;
1239 Tensor3 t3_dummy;
1240
1241 Vector p_grid;
1242 Vector lat_grid;
1243 Vector lon_grid;
1244 Tensor3 t_field;
1245 Tensor3 z_field;
1246 Tensor4 vmr_field;
1247 Tensor4 particle_bulkprop_field;
1248 ArrayOfString particle_bulkprop_names;
1249 GriddedField4 atm_fields_compact;
1250 Index abs_f_interp_order{0};
1251
1252 // Extract fields from atmfield and check their validity.
1253 // This closes the loophole when only calculating lookup tables.
1254 atm_fields_compact = batch_fields[i];
1255
1257 lat_grid,
1258 lon_grid,
1259 t_field,
1260 z_field,
1261 vmr_field,
1262 particle_bulkprop_field,
1263 particle_bulkprop_names,
1264 abs_species,
1265 atm_fields_compact,
1266 atmosphere_dim,
1267 "-",
1268 0,
1269 check_gridnames,
1270 verbosity);
1271
1272 try {
1273 atmfields_checkedCalc(atmfields_checked,
1274 atmosphere_dim,
1275 p_grid,
1276 lat_grid,
1277 lon_grid,
1278 abs_species,
1279 t_field,
1280 vmr_field,
1281 t3_dummy,
1282 t3_dummy,
1283 t3_dummy,
1284 t3_dummy,
1285 t3_dummy,
1286 t3_dummy,
1287 abs_f_interp_order,
1288 0,
1289 verbosity);
1290 } catch (const std::exception& e) {
1291 // If `robust`, skip field and continue, ...
1292 if (robust) {
1293 out1 << " WARNING! Skipped invalid atmfield "
1294 << "at batch_atmfield index " << i << ".\n"
1295 << "The runtime error produced was:\n"
1296 << e.what() << "\n";
1297 continue;
1298 }
1299 // ... else throw an error.
1300 else {
1301 stringstream err;
1302 err << "Invalid atmfield at batch_atmfield index " << i << ".\n"
1303 << "The runtime error produced was:\n"
1304 << e.what() << "\n";
1305 throw std::runtime_error(err.str());
1306 }
1307 };
1308 valid_field_indices.push_back(i); // Append index to list of valid fields.
1309
1310 if (maxp < batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)[0])
1311 maxp = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)[0];
1312 if (minp >
1313 batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)
1314 [batch_fields[i].get_numeric_grid(GFIELD4_P_GRID).nelem() - 1])
1315 minp = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID)
1316 [batch_fields[i].get_numeric_grid(GFIELD4_P_GRID).nelem() - 1];
1317 }
1318 // cout << " minp/maxp: " << minp << " / " << maxp << "\n";
1319
1320 // Information on the number of skipped atmospheres.
1321 if (batch_fields.nelem() > valid_field_indices.nelem()) {
1322 out1 << " " << batch_fields.nelem() - valid_field_indices.nelem()
1323 << " out of " << batch_fields.nelem() << " atmospheres ignored.\n";
1324 }
1325
1326 // Throw error if no atmfield passed the check.
1327 if (valid_field_indices.nelem() < 1) {
1328 stringstream err;
1329 err << "You need at least one valid profile.\n"
1330 << "It seems that no atmfield passed the checks!\n";
1331 throw std::runtime_error(err.str());
1332 }
1333
1334 if (maxp == minp) {
1335 ostringstream os;
1336 os << "You need at least two pressure levels.";
1337 throw runtime_error(os.str());
1338 }
1339
1340 // We construct the pressure grid as follows:
1341 // - Everything is done in log(p).
1342 // - Start with maxp and go down in steps of p_step until we are <= minp.
1343 // - Adjust the final pressure value to be exactly minp, otherwise
1344 // we have problems in getting min, max, and mean values for this
1345 // point later.
1346 Index np = (Index)ceil((log(maxp) - log(minp)) / p_step) + 1;
1347 // The +1 above has to be there because we must include also both
1348 // grid end points.
1349
1350 // If np is too small for the interpolation order, we increase it:
1351 if (np < abs_p_interp_order + 1) np = abs_p_interp_order + 1;
1352
1353 Vector log_abs_p(log(maxp), np, -p_step);
1354 log_abs_p[np - 1] = log(minp);
1355
1356 abs_p.resize(np);
1357 transform(abs_p, exp, log_abs_p);
1358 out2 << " abs_p: " << abs_p[0] << " Pa to " << abs_p[np - 1]
1359 << " Pa in log10 steps of " << p_step10 << " (" << np
1360 << " grid points)\n";
1361
1362 // Now we have to determine the statistics of T and VMRs, we need
1363 // profiles of min, max, and mean of these, on the abs_p grid.
1364
1365 // Index n_variables = batch_fields[0].data.nbooks();
1366 // we will do statistics for data fields excluding the scat_species fields
1367 Index n_variables = 2 + abs_species.nelem();
1368
1369 // The first dimension of datamin, datamax, and datamean is the
1370 // variable (T,Z,H2O,O3,...). The second dimension is pressure. We
1371 // assume all elements of the batch have the same variables.
1372
1373 Matrix datamin(n_variables, np, numeric_limits<Numeric>::max());
1374 Matrix datamax(n_variables, np, numeric_limits<Numeric>::min());
1375 // The limits here are from the header file <limits>
1376 Matrix datamean(n_variables, np, 0);
1377 Vector mean_norm(np, 0); // Divide by this to get mean.
1378
1379 // We will loop over all batch cases to analyze the statistics and
1380 // calculate reference profiles. As a little side project, we will
1381 // also calculate the absolute min and max of temperature and
1382 // humidity. This is handy as input to abs_lookupSetupWide.
1383 Numeric mint = +1e99;
1384 Numeric maxt = -1e99;
1385 Numeric minh2o = +1e99;
1386 Numeric maxh2o = -1e99;
1387
1388 // Loop over valid atmfields.
1389 for (Index vi = 0; vi < valid_field_indices.nelem(); ++vi) {
1390 // Get internal field index.
1391 Index i = valid_field_indices[vi];
1392
1393 // Check that really each case has the same variables (see
1394 // comment above.)
1395 if (batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES) !=
1396 batch_fields[0].get_string_grid(GFIELD4_FIELD_NAMES))
1397 throw runtime_error(
1398 "All batch atmospheres must contain the same field names.");
1399
1400 for (Index j = 0;
1401 j < batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES).nelem();
1402 ++j)
1403 if (batch_fields[i].get_string_grid(GFIELD4_FIELD_NAMES)[j] !=
1404 batch_fields[0].get_string_grid(GFIELD4_FIELD_NAMES)[j])
1405 throw runtime_error(
1406 "All batch atmospheres must contain the same individual field names.");
1407
1408 // Create convenient handles:
1409 const Vector& p_grid = batch_fields[i].get_numeric_grid(GFIELD4_P_GRID);
1410 const Tensor4& data = batch_fields[i].data;
1411
1412 // Update our global max/min values for T and H2O:
1413
1414 // We have to loop over pressure, latitudes, and longitudes
1415 // here. The dimensions of data are:
1416 // data[N_fields][N_p][N_lat][N_lon]
1417
1418 for (Index ip = 0; ip < data.npages(); ++ip)
1419 for (Index ilat = 0; ilat < data.nrows(); ++ilat)
1420 for (Index ilon = 0; ilon < data.ncols(); ++ilon) {
1421 // Field T_index is temperature:
1422 if (data(T_index, ip, ilat, ilon) < mint)
1423 mint = data(T_index, ip, ilat, ilon);
1424 if (data(T_index, ip, ilat, ilon) > maxt)
1425 maxt = data(T_index, ip, ilat, ilon);
1426 // Field batch_index[h2o_index] is H2O:
1427 if (data(batch_index[h2o_index], ip, ilat, ilon) < minh2o) {
1428 minh2o = data(batch_index[h2o_index], ip, ilat, ilon);
1429 }
1430 if (data(batch_index[h2o_index], ip, ilat, ilon) > maxh2o) {
1431 maxh2o = data(batch_index[h2o_index], ip, ilat, ilon);
1432 }
1433 }
1434
1435 // Interpolate the current batch fields to the abs_p grid. We
1436 // have to do this for each batch case, since the grids could
1437 // all be different.
1438
1439 Vector log_p_grid(p_grid.nelem());
1440 transform(log_p_grid, log, p_grid);
1441
1442 // There is a catch here: We can only use the part of abs_p that
1443 // is inside the current pressure grid p_grid, otherwise we
1444 // would have to extrapolate.
1445 // The eps_bottom and eps_top account for the little bit of
1446 // extrapolation that is allowed.
1447
1448 const Numeric eps_bottom = (log_p_grid[0] - log_p_grid[1]) / 2.1;
1449 Index first_p = 0;
1450 while (log_abs_p[first_p] > log_p_grid[0] + eps_bottom) ++first_p;
1451
1452 const Numeric eps_top = (log_p_grid[log_p_grid.nelem() - 2] -
1453 log_p_grid[log_p_grid.nelem() - 1]) /
1454 2.1;
1455 Index last_p = log_abs_p.nelem() - 1;
1456 while (log_abs_p[last_p] < log_p_grid[log_p_grid.nelem() - 1] - eps_top)
1457 --last_p;
1458
1459 const Index extent_p = last_p - first_p + 1;
1460
1461 // This was too complicated to get right:
1462 // const Index first_p = (Index) round ( (log_abs_p[0] - log_p_grid[0]) / p_step);
1463 // const Index extent_p = (Index) round ( (log_abs_p[first_p] - log_p_grid[log_p_grid.nelem()-1]) / p_step) + 1;
1464
1465 ConstVectorView active_log_abs_p = log_abs_p[Range(first_p, extent_p)];
1466
1467 // cout << "first_p / last_p / extent_p : " << first_p << " / " << last_p << " / " << extent_p << "\n";
1468 // cout << "log_p_grid: " << log_p_grid << "\n";
1469 // cout << "log_abs_p: " << log_abs_p << "\n";
1470 // cout << "active_log_abs_p: " << active_log_abs_p << "\n";
1471 // cout << "=============================================================\n";
1472 // arts_exit();
1473
1474 // Grid positions:
1475 ArrayOfGridPos gp(active_log_abs_p.nelem());
1476 gridpos(gp, log_p_grid, active_log_abs_p);
1477 // gridpos(gp, log_p_grid, active_log_abs_p, 100);
1478 // We allow much more extrapolation here than normal (0.5 is
1479 // normal). If we do not do this, then we get problems for
1480 // p_grids that are much finer than abs_p.
1481
1482 // Interpolation weights:
1483 Matrix itw(gp.nelem(), 2);
1484 interpweights(itw, gp);
1485
1486 // We have to loop over fields, latitudes, and longitudes here, doing the
1487 // pressure interpolation for all. The dimensions of data are:
1488 // data[N_fields][N_p][N_lat][N_lon]
1489 // For data_interp we reduce data by the particle related fields, i.e.,
1490 // the ones in between the first two and the last N=abs_species.nelem().
1491 // Tensor4 data_interp(data.nbooks(),
1492 Tensor4 data_interp(
1493 n_variables, active_log_abs_p.nelem(), data.nrows(), data.ncols());
1494
1495 for (Index lo = 0; lo < data.ncols(); ++lo)
1496 for (Index la = 0; la < data.nrows(); ++la) {
1497 // for (Index fi=0; fi<data.nbooks(); ++fi)
1498 // we have to handle T/z and abs_species parts separately since they
1499 // can have scat_species in between, which we want to ignore
1500 interp(data_interp(0, joker, la, lo),
1501 itw,
1502 data(T_index, joker, la, lo),
1503 gp);
1504 interp(data_interp(1, joker, la, lo),
1505 itw,
1506 data(z_index, joker, la, lo),
1507 gp);
1508 for (Index fi = 0; fi < abs_species.nelem(); ++fi)
1509 interp(data_interp(fi + 2, joker, la, lo),
1510 itw,
1511 data(batch_index[fi], joker, la, lo),
1512 gp);
1513 }
1514
1515 // Now update our datamin, datamax, and datamean variables.
1516 // We need the min and max only for the T and H2O profile,
1517 // not for others. But we need the mean for all. We are just
1518 // hopping over the case that we do not need below. This is not
1519 // very clean, but efficient. And it avoids handling all the
1520 // different cases explicitly.
1521 for (Index lo = 0; lo < data_interp.ncols(); ++lo)
1522 for (Index la = 0; la < data_interp.nrows(); ++la) {
1523 for (Index fi = 0; fi < data_interp.nbooks(); ++fi) {
1524 if (1 != fi) // We skip the z field, which we do not need
1525 for (Index pr = 0; pr < data_interp.npages(); ++pr) {
1526 // Min and max only needed for T and H2o
1527 if (0 == fi || (h2o_index + 2) == fi) {
1528 if (data_interp(fi, pr, la, lo) < datamin(fi, first_p + pr))
1529 datamin(fi, first_p + pr) = data_interp(fi, pr, la, lo);
1530 if (data_interp(fi, pr, la, lo) > datamax(fi, first_p + pr))
1531 datamax(fi, first_p + pr) = data_interp(fi, pr, la, lo);
1532 }
1533
1534 datamean(fi, first_p + pr) += data_interp(fi, pr, la, lo);
1535 }
1536 }
1537
1538 // The mean_norm is actually a bit tricky. It depends on
1539 // pressure, since different numbers of cases contribute
1540 // to the mean for different pressures. At the very eges
1541 // of the grid, typically only a single case contributes.
1542
1543 mean_norm[Range(first_p, extent_p)] += 1;
1544 }
1545 }
1546
1547 out2 << " Global statistics:\n"
1548 << " min(p) / max(p) [Pa]: " << minp << " / " << maxp << "\n"
1549 << " min(T) / max(T) [K]: " << mint << " / " << maxt << "\n"
1550 << " min(H2O) / max(H2O) [VMR]: " << minh2o << " / " << maxh2o << "\n";
1551
1552 // Divide mean by mean_norm to get the mean:
1553 ARTS_ASSERT(np == mean_norm.nelem());
1554 for (Index fi = 0; fi < datamean.nrows(); ++fi)
1555 if (1 != fi) // We skip the z field, which we do not need
1556 for (Index pi = 0; pi < np; ++pi) {
1557 // We do this in an explicit loop here, since we have to
1558 // check whether there really were data points to calculate
1559 // the mean at each level.
1560 if (0 < mean_norm[pi])
1561 datamean(fi, pi) /= mean_norm[pi];
1562 else {
1563 ostringstream os;
1564 os << "No data at pressure level " << pi + 1 << " of " << np << " ("
1565 << abs_p[pi] << " hPa).";
1566 throw runtime_error(os.str());
1567 }
1568 }
1569
1570 // If we do higher order pressure interpolation, then we should
1571 // smooth the reference profiles with a boxcar filter of width
1572 // p_interp_order+1. Otherwise we get numerical problems if there
1573 // are any sharp spikes in the reference profiles.
1574 ARTS_ASSERT(log_abs_p.nelem() == np);
1575 Matrix smooth_datamean(datamean.nrows(), datamean.ncols(), 0);
1576 for (Index i = 0; i < np; ++i) {
1577 const Index idx0 = Interpolation::pos_finder(i, log_abs_p[i], log_abs_p, abs_p_interp_order, false, false);
1578
1579 for (Index fi = 0; fi < datamean.nrows(); ++fi)
1580 if (1 != fi) // We skip the z field, which we do not need
1581 {
1582 for (Index j = 0; j < abs_p_interp_order+1; ++j) {
1583 smooth_datamean(fi, i) += datamean(fi, idx0 + j);
1584 }
1585 smooth_datamean(fi, i) /= Numeric(abs_p_interp_order + 1);
1586 }
1587 // cout << "H2O-raw / H2O-smooth: "
1588 // << datamean(h2o_index+2,i) << " / "
1589 // << smooth_datamean(h2o_index+2,i) << "\n";
1590 }
1591
1592 // There is another complication: If the (smoothed) mean for the H2O
1593 // reference profile is 0, then we have to adjust both mean and max
1594 // value to a non-zero number, otherwise the reference profile will
1595 // be zero, and we will get numerical problems later on when we
1596 // divide by the reference value. So, we set it here to 1e-9.
1597 for (Index i = 0; i < np; ++i) {
1598 // Assert that really H2O has index batch_index[h2o_index] VMR field list
1599 // and h2o_index in abs_species
1601 species_type, field_names[batch_index[h2o_index]], delim);
1602 ARTS_ASSERT(species_type == "abs_species");
1604 species_name, field_names[batch_index[h2o_index]], delim);
1605 ARTS_ASSERT("H2O" == species_name);
1606 ARTS_ASSERT("H2O" == species_names[h2o_index]);
1607
1608 // Find mean and max H2O for this level:
1609 Numeric& mean_h2o = smooth_datamean(h2o_index + 2, i);
1610 Numeric& max_h2o = datamax(h2o_index + 2, i);
1611 if (mean_h2o <= 0) {
1612 mean_h2o = 1e-9;
1613 max_h2o = 1e-9;
1614 out3 << " H2O profile contained zero values, adjusted to 1e-9.\n";
1615 }
1616 }
1617
1618 // Set abs_t:
1619 abs_t.resize(np);
1620 abs_t = smooth_datamean(0, joker);
1621 // cout << "abs_t: " << abs_t << "\n\n";
1622 // cout << "tmin: " << datamin(0,joker) << "\n\n";
1623 // cout << "tmax: " << datamax(0,joker) << "\n";
1624
1625 // Set abs_vmrs:
1626 ARTS_ASSERT(abs_species.nelem() == smooth_datamean.nrows() - 2);
1627 abs_vmrs.resize(abs_species.nelem(), np);
1628 abs_vmrs = smooth_datamean(Range(2, abs_species.nelem()), joker);
1629 // cout << "\n\nabs_vmrs: " << abs_vmrs << "\n\n";
1630
1631 // Construct abs_t_pert:
1632 ConstVectorView tmin = datamin(0, joker);
1633 ConstVectorView tmax = datamax(0, joker);
1634 choose_abs_t_pert(abs_t_pert,
1635 abs_t,
1636 tmin,
1637 tmax,
1638 t_step,
1639 abs_p_interp_order,
1640 abs_t_interp_order,
1641 verbosity);
1642 // cout << "abs_t_pert: " << abs_t_pert << "\n";
1643
1644 // Construct abs_nls_pert:
1645 ConstVectorView h2omin = datamin(h2o_index + 2, joker);
1646 ConstVectorView h2omax = datamax(h2o_index + 2, joker);
1647 choose_abs_nls_pert(abs_nls_pert,
1648 abs_vmrs(h2o_index, joker),
1649 h2omin,
1650 h2omax,
1651 h2o_step,
1652 abs_p_interp_order,
1653 abs_nls_interp_order,
1654 verbosity);
1655 // cout << "abs_nls_pert: " << abs_nls_pert << "\n";
1656
1657 // Append the explicitly given user extreme values, if necessary:
1658 if (0 != extremes.nelem()) {
1659 // There must be 4 values in this case: t_min, t_max, h2o_min, h2o_max
1660 if (4 != extremes.nelem()) {
1661 ostringstream os;
1662 os << "There must be exactly 4 elements in extremes:\n"
1663 << "min(abs_t_pert), max(abs_t_pert), min(abs_nls_pert), max(abs_nls_pert)";
1664 throw runtime_error(os.str());
1665 }
1666
1667 // t_min:
1668 if (extremes[0] < abs_t_pert[0]) {
1669 Vector dummy = abs_t_pert;
1670 abs_t_pert.resize(abs_t_pert.nelem() + 1);
1671 abs_t_pert[0] = extremes[0];
1672 abs_t_pert[Range(1, dummy.nelem())] = dummy;
1673 out2 << " Added min extreme value for abs_t_pert: " << abs_t_pert[0]
1674 << "\n";
1675 }
1676
1677 // t_max:
1678 if (extremes[1] > abs_t_pert[abs_t_pert.nelem() - 1]) {
1679 Vector dummy = abs_t_pert;
1680 abs_t_pert.resize(abs_t_pert.nelem() + 1);
1681 abs_t_pert[Range(0, dummy.nelem())] = dummy;
1682 abs_t_pert[abs_t_pert.nelem() - 1] = extremes[1];
1683 out2 << " Added max extreme value for abs_t_pert: "
1684 << abs_t_pert[abs_t_pert.nelem() - 1] << "\n";
1685 }
1686
1687 // nls_min:
1688 if (extremes[2] < abs_nls_pert[0]) {
1689 Vector dummy = abs_nls_pert;
1690 abs_nls_pert.resize(abs_nls_pert.nelem() + 1);
1691 abs_nls_pert[0] = extremes[2];
1692 abs_nls_pert[Range(1, dummy.nelem())] = dummy;
1693 out2 << " Added min extreme value for abs_nls_pert: " << abs_nls_pert[0]
1694 << "\n";
1695 }
1696
1697 // nls_max:
1698 if (extremes[3] > abs_nls_pert[abs_nls_pert.nelem() - 1]) {
1699 Vector dummy = abs_nls_pert;
1700 abs_nls_pert.resize(abs_nls_pert.nelem() + 1);
1701 abs_nls_pert[Range(0, dummy.nelem())] = dummy;
1702 abs_nls_pert[abs_nls_pert.nelem() - 1] = extremes[3];
1703 out2 << " Added max extreme value for abs_nls_pert: "
1704 << abs_nls_pert[abs_nls_pert.nelem() - 1] << "\n";
1705 }
1707}
1708
1709void abs_lookupSetupWide( // WS Output:
1710 Vector& abs_p,
1711 Vector& abs_t,
1712 Vector& abs_t_pert,
1713 Matrix& abs_vmrs,
1714 ArrayOfArrayOfSpeciesTag& abs_nls,
1715 Vector& abs_nls_pert,
1716 // WS Input:
1717 const ArrayOfArrayOfSpeciesTag& abs_species,
1718 const Index& abs_p_interp_order,
1719 const Index& abs_t_interp_order,
1720 const Index& abs_nls_interp_order,
1721 // Control Parameters:
1722 const Numeric& p_min,
1723 const Numeric& p_max,
1724 const Numeric& p_step10,
1725 const Numeric& t_min,
1726 const Numeric& t_max,
1727 const Numeric& h2o_min,
1728 const Numeric& h2o_max,
1729 const Verbosity& verbosity) {
1731
1732 // For consistency with other code around arts (e.g., correlation
1733 // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
1734 // convert it here to the natural log:
1735 const Numeric p_step = log(pow(10.0, p_step10));
1736
1737 // Make an intelligent choice for the nonlinear species.
1738 choose_abs_nls(abs_nls, abs_species, verbosity);
1739
1740 // 1. Fix pressure grid abs_p
1741 // --------------------------
1742
1743 // We construct the pressure grid as follows:
1744 // - Everything is done in log(p).
1745 // - Start with p_max and go down in steps of p_step until we are <= p_min.
1746
1747 Index np = (Index)ceil((log(p_max) - log(p_min)) / p_step) + 1;
1748 // The +1 above has to be there because we must include also both
1749 // grid end points.
1750
1751 // If np is too small for the interpolation order, we increase it:
1752 if (np < abs_p_interp_order + 1) np = abs_p_interp_order + 1;
1753
1754 Vector log_abs_p(log(p_max), np, -p_step);
1755
1756 abs_p.resize(np);
1757 transform(abs_p, exp, log_abs_p);
1758 out2 << " abs_p: " << abs_p[0] << " Pa to " << abs_p[np - 1]
1759 << " Pa in log10 steps of " << p_step10 << " (" << np
1760 << " grid points)\n";
1761
1762 // 2. Fix reference temperature profile abs_t and temperature perturbations
1763 // ------------------------------------------------------------------------
1764
1765 // We simply take a constant reference profile.
1766
1767 Numeric t_ref = (t_min + t_max) / 2;
1768
1769 abs_t.resize(np);
1770 abs_t = t_ref; // Assign same value to all elements.
1771
1772 // We have to make vectors out of t_min and t_max, so we can use
1773 // them in the choose_abs_t_pert function call:
1774 Vector min_prof(np), max_prof(np);
1775 min_prof = t_min;
1776 max_prof = t_max;
1777
1778 // Chose temperature perturbations:
1779 choose_abs_t_pert(abs_t_pert,
1780 abs_t,
1781 min_prof,
1782 max_prof,
1783 20,
1784 abs_p_interp_order,
1785 abs_t_interp_order,
1786 verbosity);
1787
1788 // 3. Fix reference H2O profile and abs_nls_pert
1789 // ---------------------------------------------
1790
1791 // We take a constant reference profile of 1000ppm (=1e-3) for H2O
1792 Numeric const h2o_ref = 1e-3;
1793
1794 // And 1 ppt (1e-9) as default for all VMRs
1795 Numeric const other_ref = 1e-9;
1796
1797 // We have to assign this value to all pressures of the H2O profile,
1798 // and 0 to all other profiles.
1799
1800 // abs_vmrs has dimension [n_species, np].
1801 abs_vmrs.resize(abs_species.nelem(), np);
1802 abs_vmrs = other_ref;
1803
1804 // We look for O2 and N2, and assign constant values to them.
1805 // The values are from Wallace&Hobbs, 2nd edition.
1806
1807 const Index o2_index =
1808 find_first_species(abs_species, Species::fromShortName("O2"));
1809 if (o2_index >= 0) {
1810 abs_vmrs(o2_index, joker) = 0.2095;
1811 }
1812
1813 const Index n2_index =
1814 find_first_species(abs_species, Species::fromShortName("N2"));
1815 if (n2_index >= 0) {
1816 abs_vmrs(n2_index, joker) = 0.7808;
1817 }
1818
1819 // Which species is H2O?
1820 const Index h2o_index = find_first_species(
1821 abs_species, Species::fromShortName("H2O"));
1822
1823 // The function returns "-1" if there is no H2O
1824 // species.
1825 if (0 < abs_nls.nelem()) {
1826 if (h2o_index < 0) {
1827 ostringstream os;
1828 os << "Some of your species require nonlinear treatment,\n"
1829 << "but you have no H2O species.";
1830 throw runtime_error(os.str());
1831 }
1832
1833 // Assign constant reference value to all H2O levels:
1834 abs_vmrs(h2o_index, joker) = h2o_ref;
1835
1836 // We have to make vectors out of h2o_min and h2o_max, so we can use
1837 // them in the choose_abs_nls_pert function call.
1838 // We re-use the vectors min_prof and max_prof that we have
1839 // defined above.
1840 min_prof = h2o_min;
1841 max_prof = h2o_max;
1842
1843 // Construct abs_nls_pert:
1844 choose_abs_nls_pert(abs_nls_pert,
1845 abs_vmrs(h2o_index, joker),
1846 min_prof,
1847 max_prof,
1848 1e99,
1849 abs_p_interp_order,
1850 abs_nls_interp_order,
1851 verbosity);
1852 } else {
1854 out1 << " WARNING:\n"
1855 << " You have no species that require H2O variations.\n"
1856 << " This case might work, but it has never been tested.\n"
1857 << " Please test it, then remove this warning.\n";
1858 }
1860
1861/* Workspace method: Doxygen documentation will be auto-generated */
1862void abs_speciesAdd( // WS Output:
1863 ArrayOfArrayOfSpeciesTag& abs_species,
1864 Index& propmat_clearsky_agenda_checked,
1865 // Control Parameters:
1866 const ArrayOfString& names,
1867 const Verbosity& verbosity) {
1869
1870 // Invalidate agenda check flags
1871 propmat_clearsky_agenda_checked = false;
1872
1873 // Size of initial array
1874 Index n_gs = abs_species.nelem();
1875
1876 // Temporary ArrayOfSpeciesTag
1877 ArrayOfSpeciesTag temp;
1878
1879 // Each element of the array of Strings names defines one tag
1880 // group. Let's work through them one by one.
1881 for (Index i = 0; i < names.nelem(); ++i) {
1882 abs_species.emplace_back(names[i]);
1883 }
1884
1885 check_abs_species(abs_species);
1886
1887 // Print list of tag groups to the most verbose output stream:
1888 out3 << " Added tag groups:";
1889 for (Index i = n_gs; i < abs_species.nelem(); ++i) {
1890 out3 << "\n " << i << ":";
1891 for (Index s = 0; s < abs_species[i].nelem(); ++s) {
1892 out3 << " " << abs_species[i][s].Name();
1893 }
1894 }
1895 out3 << '\n';
1897
1898/* Workspace method: Doxygen documentation will be auto-generated */
1899void abs_speciesAdd2( // WS Output:
1900 Workspace& ws,
1901 ArrayOfArrayOfSpeciesTag& abs_species,
1903 Agenda& jacobian_agenda,
1904 Index& propmat_clearsky_agenda_checked,
1905 // WS Input:
1906 const Index& atmosphere_dim,
1907 const Vector& p_grid,
1908 const Vector& lat_grid,
1909 const Vector& lon_grid,
1910 // WS Generic Input:
1911 const Vector& rq_p_grid,
1912 const Vector& rq_lat_grid,
1913 const Vector& rq_lon_grid,
1914 // Control Parameters:
1915 const String& species,
1916 const String& mode,
1917 const Verbosity& verbosity) {
1919
1920 // Invalidate agenda check flags
1921 propmat_clearsky_agenda_checked = false;
1922
1923 // Add species to *abs_species*
1924 abs_species.emplace_back(species);
1925
1926 check_abs_species(abs_species);
1927
1928 // Print list of added tag group to the most verbose output stream:
1929 out3 << " Appended tag group:";
1930 out3 << "\n " << abs_species.nelem() - 1 << ":";
1931 for (Index s = 0; s < abs_species.back().nelem(); ++s) {
1932 out3 << " " << abs_species.back()[s].Name();
1933 }
1934 out3 << '\n';
1935
1936 // Do retrieval part
1938 jq,
1939 jacobian_agenda,
1940 atmosphere_dim,
1941 p_grid,
1942 lat_grid,
1943 lon_grid,
1944 rq_p_grid,
1945 rq_lat_grid,
1946 rq_lon_grid,
1947 species,
1948 mode,
1949 1,
1950 verbosity);
1952
1953/* Workspace method: Doxygen documentation will be auto-generated */
1954void abs_speciesInit(ArrayOfArrayOfSpeciesTag& abs_species, const Verbosity&) {
1955 abs_species.resize(0);
1957
1958/* Workspace method: Doxygen documentation will be auto-generated */
1959void abs_speciesSet( // WS Output:
1960 ArrayOfArrayOfSpeciesTag& abs_species,
1961 Index& propmat_clearsky_agenda_checked,
1962 // Control Parameters:
1963 const ArrayOfString& names,
1964 const Verbosity& verbosity) {
1966
1967 // Invalidate agenda check flags
1968 propmat_clearsky_agenda_checked = false;
1969
1970 abs_species.resize(names.nelem());
1971
1972 //cout << "Names: " << names << "\n";
1973
1974 // Each element of the array of Strings names defines one tag
1975 // group. Let's work through them one by one.
1976 for (Index i = 0; i < names.nelem(); ++i) {
1977 // This part has now been moved to array_species_tag_from_string.
1978 // Call this function.
1979 abs_species[i] = ArrayOfSpeciesTag(names[i]);
1980 }
1981
1982 check_abs_species(abs_species);
1983
1984 // Print list of tag groups to the most verbose output stream:
1985 out3 << " Defined tag groups: ";
1986 for (Index i = 0; i < abs_species.nelem(); ++i) {
1987 out3 << "\n " << i << ":";
1988 for (Index s = 0; s < abs_species[i].nelem(); ++s)
1989 out3 << " " << abs_species[i][s].Name();
1990 }
1991 out3 << '\n';
1993
1994/* Workspace method: Doxygen documentation will be auto-generated */
1995void abs_lookupAdapt(GasAbsLookup& abs_lookup,
1996 Index& abs_lookup_is_adapted,
1997 const ArrayOfArrayOfSpeciesTag& abs_species,
1998 const Vector& f_grid,
1999 const Verbosity& verbosity) {
2000 abs_lookup.Adapt(abs_species, f_grid, verbosity);
2001 abs_lookup_is_adapted = 1;
2003
2004/* Workspace method: Doxygen documentation will be auto-generated */
2006 PropagationMatrix& propmat_clearsky,
2007 ArrayOfPropagationMatrix& dpropmat_clearsky_dx,
2008 const GasAbsLookup& abs_lookup,
2009 const Index& abs_lookup_is_adapted,
2010 const Index& abs_p_interp_order,
2011 const Index& abs_t_interp_order,
2012 const Index& abs_nls_interp_order,
2013 const Index& abs_f_interp_order,
2014 const Vector& f_grid,
2015 const Numeric& a_pressure,
2016 const Numeric& a_temperature,
2017 const Vector& a_vmr_list,
2018 const ArrayOfRetrievalQuantity& jacobian_quantities,
2019 const ArrayOfArrayOfSpeciesTag& abs_species,
2020 const ArrayOfSpeciesTag& select_abs_species,
2021 const Numeric& extpolfac,
2022 const Index& no_negatives,
2023 const Verbosity& verbosity) {
2025
2026 // Variables needed by abs_lookup.Extract:
2027 Matrix abs_scalar_gas, dabs_scalar_gas_df, dabs_scalar_gas_dt;
2028
2029 // Check if the table has been adapted:
2030 if (1 != abs_lookup_is_adapted)
2031 throw runtime_error(
2032 "Gas absorption lookup table must be adapted,\n"
2033 "use method abs_lookupAdapt.");
2034
2035 const bool do_jac = supports_lookup(jacobian_quantities);
2036 const bool do_freq_jac = do_frequency_jacobian(jacobian_quantities);
2037 const bool do_temp_jac = do_temperature_jacobian(jacobian_quantities);
2038 const Numeric df = frequency_perturbation(jacobian_quantities);
2039 const Numeric dt = temperature_perturbation(jacobian_quantities);
2040
2041 // The combination of doing frequency jacobian together with an
2042 // absorption lookup table is quite dangerous. If the frequency
2043 // interpolation order for the table is zero, the Jacobian will be
2044 // zero, and the cause for this may be difficult for a user to
2045 // find. So we do not allow this combination.
2046 if (do_freq_jac and (1 > abs_f_interp_order))
2047 throw std::runtime_error("Wind/frequency Jacobian is not possible without at least first\n"
2048 "order frequency interpolation in the lookup table. Please use\n"
2049 "abs_f_interp_order>0 or remove wind/frequency Jacobian.");
2050
2051 // The function we are going to call here is one of the few helper
2052 // functions that adjust the size of their output argument
2053 // automatically.
2054 abs_lookup.Extract(abs_scalar_gas,
2055 select_abs_species,
2056 abs_p_interp_order,
2057 abs_t_interp_order,
2058 abs_nls_interp_order,
2059 abs_f_interp_order,
2060 a_pressure,
2061 a_temperature,
2062 a_vmr_list,
2063 f_grid,
2064 extpolfac);
2065 if (do_freq_jac) {
2066 Vector dfreq = f_grid;
2067 dfreq += df;
2068 abs_lookup.Extract(dabs_scalar_gas_df,
2069 select_abs_species,
2070 abs_p_interp_order,
2071 abs_t_interp_order,
2072 abs_nls_interp_order,
2073 abs_f_interp_order,
2074 a_pressure,
2075 a_temperature,
2076 a_vmr_list,
2077 dfreq,
2078 extpolfac);
2079 }
2080 if (do_temp_jac) {
2081 const Numeric dtemp = a_temperature + dt;
2082 abs_lookup.Extract(dabs_scalar_gas_dt,
2083 select_abs_species,
2084 abs_p_interp_order,
2085 abs_t_interp_order,
2086 abs_nls_interp_order,
2087 abs_f_interp_order,
2088 a_pressure,
2089 dtemp,
2090 a_vmr_list,
2091 f_grid,
2092 extpolfac);
2093 }
2094
2095 if (no_negatives){
2096 //Check for negative values due to interpolation and set them to zero
2097 for (Index ir = 0; ir < abs_scalar_gas.nrows(); ir++){
2098 for (Index ic = 0; ic < abs_scalar_gas.ncols(); ic++){
2099 if (abs_scalar_gas(ir, ic)<0) abs_scalar_gas(ir, ic)=0;
2100 }
2101 }
2102 }
2103
2104 // Now add to the right place in the absorption matrix.
2105
2106 if (not do_jac) {
2107 for (Index ii = 0; ii < abs_scalar_gas.nrows(); ii++) {
2108 propmat_clearsky.Kjj() += abs_scalar_gas(ii, joker);
2109 }
2110 } else {
2111 for (Index isp = 0; isp < abs_scalar_gas.nrows(); isp++) {
2112 propmat_clearsky.Kjj() += abs_scalar_gas(isp, joker);
2113
2114 for (Index iv = 0; iv < propmat_clearsky.NumberOfFrequencies();
2115 iv++) {
2116 for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {
2117 const auto& deriv = jacobian_quantities[iq];
2118
2119 if (not deriv.propmattype()) continue;
2120
2121 if (deriv == Jacobian::Atm::Temperature) {
2122 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
2123 (dabs_scalar_gas_dt(isp, iv) - abs_scalar_gas(isp, iv)) / dt;
2124 } else if (is_frequency_parameter(deriv)) {
2125 dpropmat_clearsky_dx[iq].Kjj()[iv] +=
2126 (dabs_scalar_gas_df(isp, iv) - abs_scalar_gas(isp, iv)) / df;
2127 } else if (deriv == abs_species[isp]) {
2128 // WARNING: If CIA in list, this scales wrong by factor 2
2129 dpropmat_clearsky_dx[iq].Kjj()[iv] += (std::isnormal(a_vmr_list[isp])) ? abs_scalar_gas(isp, iv) / a_vmr_list[isp] : 0;
2130 }
2131 }
2132 }
2133 }
2134 }
2136
2137/* Workspace method: Doxygen documentation will be auto-generated */
2139 // WS Output:
2140 Tensor7& propmat_clearsky_field,
2141 // WS Input:
2142 const Index& atmfields_checked,
2143 const Vector& f_grid,
2144 const Index& stokes_dim,
2145 const Vector& p_grid,
2146 const Vector& lat_grid,
2147 const Vector& lon_grid,
2148 const Tensor3& t_field,
2149 const Tensor4& vmr_field,
2150 const EnergyLevelMap& nlte_field,
2151 const Tensor3& mag_u_field,
2152 const Tensor3& mag_v_field,
2153 const Tensor3& mag_w_field,
2154 const ArrayOfArrayOfSpeciesTag& abs_species,
2155 const Agenda& abs_agenda,
2156 // WS Generic Input:
2157 const Vector& doppler,
2158 const Vector& los,
2159 const Verbosity& verbosity) {
2162
2163 chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
2164 if (atmfields_checked != 1)
2165 throw runtime_error(
2166 "The atmospheric fields must be flagged to have "
2167 "passed a consistency check (atmfields_checked=1).");
2168
2169 ArrayOfPropagationMatrix partial_abs;
2170 ArrayOfStokesVector partial_nlte; // FIXME: This is not stored!
2172 StokesVector nlte;
2173 Vector a_vmr_list;
2174 EnergyLevelMap a_nlte_list;
2175
2176 // Get the number of species from the leading dimension of vmr_field:
2177 const Index n_species = vmr_field.nbooks();
2178
2179 // Number of frequencies:
2180 const Index n_frequencies = f_grid.nelem();
2181
2182 // Number of pressure levels:
2183 const Index n_pressures = p_grid.nelem();
2184
2185 // Number of latitude grid points (must be at least one):
2186 const Index n_latitudes = max(Index(1), lat_grid.nelem());
2187
2188 // Number of longitude grid points (must be at least one):
2189 const Index n_longitudes = max(Index(1), lon_grid.nelem());
2190
2191 // Check that doppler is empty or matches p_grid
2192 if (0 != doppler.nelem() && p_grid.nelem() != doppler.nelem()) {
2193 ostringstream os;
2194 os << "Variable doppler must either be empty, or match the dimension of "
2195 << "p_grid.";
2196 throw runtime_error(os.str());
2197 }
2198
2199 // Resize output field.
2200 // The dimension in lat and lon must be at least one, even if these
2201 // grids are empty.
2202 out2 << " Creating propmat field with dimensions:\n"
2203 << " " << n_species << " gas species,\n"
2204 << " " << n_frequencies << " frequencies,\n"
2205 << " " << stokes_dim << " stokes dimension,\n"
2206 << " " << stokes_dim << " stokes dimension,\n"
2207 << " " << n_pressures << " pressures,\n"
2208 << " " << n_latitudes << " latitudes,\n"
2209 << " " << n_longitudes << " longitudes.\n";
2210
2211 propmat_clearsky_field.resize(n_species,
2212 n_frequencies,
2213 stokes_dim,
2214 stokes_dim,
2215 n_pressures,
2216 n_latitudes,
2217 n_longitudes);
2218
2219 // Fake jacobian_quantities to deal with partial absorption
2220 ArrayOfRetrievalQuantity jacobian_quantities;
2221 for (auto& species_list: abs_species) {
2222 jacobian_quantities.emplace_back();
2223 auto& rq = jacobian_quantities.back();
2224 rq.Subtag(species_list.Name());
2225 rq.Target(Jacobian::Target(Jacobian::Special::ArrayOfSpeciesTagVMR, species_list));
2226 rq.Target().perturbation = 0.001;
2227 }
2228
2229 String fail_msg;
2230 bool failed = false;
2231
2232 // Make local copy of f_grid, so that we can apply Dopler if we want.
2233 Vector this_f_grid = f_grid;
2234
2236
2237 // Now we have to loop all points in the atmosphere:
2238 if (n_pressures)
2239#pragma omp parallel for if (!arts_omp_in_parallel() && \
2240 n_pressures >= arts_omp_get_max_threads()) \
2241 firstprivate(wss, this_f_grid) private( \
2242 abs, nlte, partial_abs, partial_nlte, a_vmr_list)
2243 for (Index ipr = 0; ipr < n_pressures; ++ipr) // Pressure: ipr
2244 {
2245 // Skip remaining iterations if an error occurred
2246 if (failed) continue;
2247
2248 // The try block here is necessary to correctly handle
2249 // exceptions inside the parallel region.
2250 try {
2251 Numeric a_pressure = p_grid[ipr];
2252
2253 if (0 != doppler.nelem()) {
2254 this_f_grid = f_grid;
2255 this_f_grid += doppler[ipr];
2256 }
2257
2258 {
2259 ostringstream os;
2260 os << " p_grid[" << ipr << "] = " << a_pressure << "\n";
2261 out3 << os.str();
2262 }
2263
2264 for (Index ila = 0; ila < n_latitudes; ++ila) // Latitude: ila
2265 for (Index ilo = 0; ilo < n_longitudes; ++ilo) // Longitude: ilo
2266 {
2267 Numeric a_temperature = t_field(ipr, ila, ilo);
2268 a_vmr_list = vmr_field(Range(joker), ipr, ila, ilo);
2269 if (!nlte_field.value.empty())
2270 a_nlte_list = nlte_field(ipr, ila, ilo);
2271
2272 Vector this_rtp_mag(3, 0.);
2273
2274 if (mag_u_field.npages() != 0) {
2275 this_rtp_mag[0] = mag_u_field(ipr, ila, ilo);
2276 }
2277 if (mag_v_field.npages() != 0) {
2278 this_rtp_mag[1] = mag_v_field(ipr, ila, ilo);
2279 }
2280 if (mag_w_field.npages() != 0) {
2281 this_rtp_mag[2] = mag_w_field(ipr, ila, ilo);
2282 }
2283
2284 // Execute agenda to calculate local absorption.
2285 // Agenda input: f_index, a_pressure, a_temperature, a_vmr_list
2286 // Agenda output: abs, nlte
2288 abs,
2289 nlte,
2290 partial_abs,
2291 partial_nlte,
2292 jacobian_quantities,
2293 {},
2294 this_f_grid,
2295 this_rtp_mag,
2296 los,
2297 a_pressure,
2298 a_temperature,
2299 a_nlte_list,
2300 a_vmr_list,
2301 abs_agenda);
2302
2303 // Convert from derivative to absorption
2304 for (Index ispec=0; ispec<partial_abs.nelem(); ispec++) {
2305 partial_abs[ispec] *= a_vmr_list[ispec];
2306 }
2307
2308 // Verify, that the number of elements in abs matrix is
2309 // constistent with stokes_dim:
2310 ARTS_USER_ERROR_IF (stokes_dim != abs.StokesDimensions(),
2311 "propmat_clearsky_fieldCalc was called with stokes_dim = ",
2312 stokes_dim, ",\n"
2313 "but the stokes_dim returned by the agenda is ",
2314 abs.StokesDimensions(), ".")
2315
2316 // Verify, that the number of species in abs is
2317 // constistent with vmr_field:
2318 ARTS_USER_ERROR_IF (n_species != partial_abs.nelem(),
2319 "The number of gas species in vmr_field is ", n_species,
2320 ",\n"
2321 "but the number of species returned by the agenda is ",
2322 partial_abs.nelem(), ".")
2323
2324 // Verify, that the number of frequencies in abs is
2325 // constistent with f_extent:
2326 ARTS_USER_ERROR_IF (n_frequencies != abs.NumberOfFrequencies(),
2327 "The number of frequencies desired is ", n_frequencies,
2328 ",\n"
2329 "but the number of frequencies returned by the agenda is ",
2330 abs.NumberOfFrequencies(), ".")
2331
2332 // Store the result in output field.
2333 // We have to transpose abs, because the dimensions of the
2334 // two variables are:
2335 // abs_field: [abs_species, f_grid, stokes, stokes, p_grid, lat_grid, lon_grid]
2336 // abs: [abs_species][f_grid, stokes, stokes]
2337 for (Index i = 0; i < partial_abs.nelem(); i++) {
2338 partial_abs[i].GetTensor3(propmat_clearsky_field(
2339 i, joker, joker, joker, ipr, ila, ilo));
2340 }
2341 }
2342 } catch (const std::runtime_error& e) {
2343#pragma omp critical(propmat_clearsky_fieldCalc_fail)
2344 {
2345 fail_msg = e.what();
2346 failed = true;
2347 }
2348 }
2349 }
2350
2351 if (failed) throw runtime_error(fail_msg);
2353
2354/* Workspace method: Doxygen documentation will be auto-generated */
2355void f_gridFromGasAbsLookup(Vector& f_grid,
2356 const GasAbsLookup& abs_lookup,
2357 const Verbosity&) {
2358 const Vector& lookup_f_grid = abs_lookup.GetFgrid();
2359 f_grid.resize(lookup_f_grid.nelem());
2360 f_grid = lookup_f_grid;
2362
2363/* Workspace method: Doxygen documentation will be auto-generated */
2364void p_gridFromGasAbsLookup(Vector& p_grid,
2365 const GasAbsLookup& abs_lookup,
2366 const Verbosity&) {
2367 const Vector& lookup_p_grid = abs_lookup.GetPgrid();
2368 p_grid.resize(lookup_p_grid.nelem());
2369 p_grid = lookup_p_grid;
2370}
Declarations required for the calculation of absorption coefficients.
Declarations for agendas.
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
The global header file for ARTS.
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 ArrayOfSpeciesTag &select_abs_species, 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:25665
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:69
String name() const
Agenda name.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index ncols() const noexcept
Definition: matpackIV.h:146
bool empty() const noexcept
Definition: matpackIV.h:152
Index nrows() const noexcept
Definition: matpackIV.h:145
Index nbooks() const noexcept
Definition: matpackIV.h:143
Index npages() const noexcept
Definition: matpackIV.h:144
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
An absorption lookup table.
const Vector & GetPgrid() const
const Vector & GetFgrid() const
void Adapt(const ArrayOfArrayOfSpeciesTag &current_species, ConstVectorView current_f_grid, const Verbosity &verbosity)
Adapt lookup table to current calculation.
void Extract(Matrix &sga, const ArrayOfSpeciesTag &select_abs_species, 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.
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
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:160
Stokes vector is as Propagation matrix but only has 4 possible values.
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4 class.
Definition: matpackIV.h:435
The Tensor7 class.
Definition: matpackVII.h:2407
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5488
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Workspace class.
Definition: workspace_ng.h:53
void parse_atmcompact_speciesname(String &species_name, const String &field_name, const String &delim)
Definition: cloudbox.cc:830
void parse_atmcompact_speciestype(String &species_type, const String &field_name, const String &delim)
Definition: cloudbox.cc:799
Internal cloudbox functions.
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
Declarations for the gas absorption lookup table.
Implementation of gridded fields.
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:1175
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
bool supports_lookup(const ArrayOfRetrievalQuantity &js)
Returns if the array supports lookup table derivatives.
Definition: jacobian.cc:1109
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1006
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1183
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1191
bool is_unique(const ArrayOfIndex &x)
Checks if an ArrayOfIndex is unique, i.e., has no duplicate values.
Definition: logic.cc:277
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_speciesAdd(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesAdd.
void f_gridFromGasAbsLookup(Vector &f_grid, const GasAbsLookup &abs_lookup, const Verbosity &)
WORKSPACE METHOD: f_gridFromGasAbsLookup.
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 &propmat_clearsky_agenda, const Numeric &lowest_vmr, const Verbosity &verbosity)
WORKSPACE METHOD: abs_lookupCalc.
Definition: m_abs_lookup.cc:65
void abs_speciesSet(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesSet.
void abs_speciesAdd2(Workspace &ws, ArrayOfArrayOfSpeciesTag &abs_species, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, Index &propmat_clearsky_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.
void find_nonlinear_continua(ArrayOfIndex &cont, const ArrayOfArrayOfSpeciesTag &abs_species, const Verbosity &verbosity)
Find continuum species in abs_species.
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 ArrayOfSpeciesTag &select_abs_species, const Numeric &extpolfac, const Index &no_negatives, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddFromLookup.
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_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 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:55
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 AtmFieldsAndParticleBulkPropFieldFromCompact(Vector &p_grid, Vector &lat_grid, Vector &lon_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, Tensor4 &particle_bulkprop_field, ArrayOfString &particle_bulkprop_names, const ArrayOfArrayOfSpeciesTag &abs_species, const GriddedField4 &atm_fields_compact, const Index &atmosphere_dim, const String &delim, const Numeric &p_min, const Index &check_gridnames, const Verbosity &)
WORKSPACE METHOD: AtmFieldsAndParticleBulkPropFieldFromCompact.
void VectorInsertGridPoints(Vector &og, const Vector &ingrid, 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 &)
WORKSPACE METHOD: atmfields_checkedCalc.
Definition: m_checked.cc:97
void jacobianAddAbsSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition: m_jacobian.cc:106
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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:1433
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT1
Definition: messages.h:204
#define CREATE_OUT3
Definition: messages.h:206
#define CREATE_OUT2
Definition: messages.h:205
Namespace to contain things required for absorption calculations.
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric e
Elementary charge convenience name [C].
constexpr Index GFIELD4_P_GRID
Global constant, Index of the pressure grid in GriddedField4.
constexpr Index GFIELD4_FIELD_NAMES
Global constant, Index of the field names in GriddedField4.
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)
Implements Zeeman modeling.
Definition: zeemandata.cc:333
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
Stuff related to the propagation matrix.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
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)
Holds all information required for individual partial derivatives.
Definition: jacobian.h:108
#define d
#define a
#define c
#define b