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