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