ARTS  2.3.1285(git:92a29ea9-dirty)
m_atmosphere.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012
2  Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3  Stefan Buehler <sbuehler@ltu.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
36 /*===========================================================================
37  === External declarations
38  ===========================================================================*/
39 
40 #include <cfloat>
41 #include <cmath>
42 #include "abs_species_tags.h"
43 #include "absorption.h"
44 #include "agenda_class.h"
45 #include "arts.h"
46 #include "auto_md.h"
47 #include "check_input.h"
48 #include "cloudbox.h"
49 #include "geodetic.h"
50 #include "global_data.h"
51 #include "gridded_fields.h"
52 #include "interpolation.h"
53 #include "interpolation_poly.h"
54 #include "linescaling.h"
55 #include "matpackIII.h"
56 #include "messages.h"
57 #include "rte.h"
58 #include "special_interp.h"
59 #include "xml_io.h"
60 
61 extern const Index GFIELD3_P_GRID;
62 extern const Index GFIELD3_LAT_GRID;
63 extern const Index GFIELD3_LON_GRID;
64 extern const Index GFIELD4_FIELD_NAMES;
65 extern const Index GFIELD4_P_GRID;
66 extern const Index GFIELD4_LAT_GRID;
67 extern const Index GFIELD4_LON_GRID;
68 
69 extern const Numeric DEG2RAD;
70 extern const Numeric PI;
71 extern const Numeric GAS_CONSTANT;
72 
74 
76 extern const Numeric EPSILON_LON_CYCLIC = 2 * DBL_EPSILON;
77 
78 /*===========================================================================
79  *=== Helper functions
80  *===========================================================================*/
81 
83 
100  Index& nf,
101  const String& name,
102  const Index& prepend,
103  const Verbosity&) {
104  // Number of fields already present:
106 
107  if (prepend) {
108  // Add name of new field to field name list:
110  as.insert(as.begin(), name);
111  } else {
112  // Add name of new field to field name list:
113  af.get_string_grid(GFIELD4_FIELD_NAMES).push_back(name);
114  }
115 
116  // Save original fields:
117  const Tensor4 dummy = af.data;
118 
119  // Adjust size:
120  af.resize(nf + 1, dummy.npages(), dummy.nrows(), dummy.ncols());
121  nf++; // set to new number of fields
122 
123  // Copy back original field:
124  af.data(Range((prepend && nf > 1) ? 1 : 0, nf - 1), joker, joker, joker) =
125  dummy;
126 }
127 
128 /*===========================================================================
129  === The functions (in alphabetical order)
130  ===========================================================================*/
131 
133 /*
134  This helper function is used by all AtmFieldPRegrid WSMs to do the common
135  calculation of the grid positions and interpolation weights for pressures.
136  It is an adaptation of GriddedFieldPRegridHelper for Tensors instead of
137  GriddedFields.
138 
139  \param[out] ing_min Index in the new grid with first value covered
140  by the old grid.
141  \param[out] ing_max Index in the new grid with last value covered
142  by the old grid.
143  \param[out] gp_p Pressure grid positions
144  \param[out] itw Interpolation weights
145  \param[in] p_grid_out New pressure grid
146  \param[in] p_grid_in Old pressure grid
147  \param[in] interp_order Interpolation order
148  \param[in] verbosity Verbosity levels
149  */
151  Index& ing_max,
152  ArrayOfGridPosPoly& gp_p,
153  Matrix& itw,
154  ConstVectorView p_grid_out,
155  ConstVectorView p_grid_in,
156  const Index& interp_order,
157  const Verbosity& verbosity) {
158  CREATE_OUT2;
159 
160  out2 << " Interpolation order: " << interp_order << "\n";
161 
162  ing_min = 0;
163  ing_max = p_grid_out.nelem() - 1;
165  "Atmospheric field to p_grid_out", p_grid_in, p_grid_out, interp_order);
166 
167  Index nelem_in_range = ing_max - ing_min + 1;
168 
169  // Calculate grid positions:
170  if (nelem_in_range > 0) {
171  gp_p.resize(nelem_in_range);
172  p2gridpos_poly(gp_p,
173  p_grid_in,
174  p_grid_out[Range(ing_min, nelem_in_range)],
175  interp_order);
176 
177  // Interpolation weights:
178  itw.resize(nelem_in_range, interp_order + 1);
179  interpweights(itw, gp_p);
180  }
181 }
182 
183 /* Workspace method: Doxygen documentation will be auto-generated */
184 void AtmFieldPRegrid( // WS Generic Output:
185  Tensor3& atmtensor_out,
186  // WS Generic Input:
187  const Tensor3& atmtensor_in_orig,
188  const Vector& p_grid_new,
189  const Vector& p_grid_old,
190  const Index& interp_order,
191  const Verbosity& verbosity) {
192  // check that p_grid_old is the p_grid associated with atmtensor_in_orig. we can
193  // only check for consistent size with p_grid_old.
194  if (atmtensor_in_orig.npages() != p_grid_old.nelem()) {
195  ostringstream os;
196  os << "p_grid_old is supposed to be the p_grid associated with the "
197  << "atmospheric field.\n"
198  << "However, it is not as their sizes are inconsistent.\n";
199  throw runtime_error(os.str());
200  }
201 
202  const Tensor3* atmtensor_in_pnt;
203  Tensor3 atmtensor_in_copy;
204 
205  if (&atmtensor_in_orig == &atmtensor_out) {
206  atmtensor_in_copy = atmtensor_in_orig;
207  atmtensor_in_pnt = &atmtensor_in_copy;
208  } else
209  atmtensor_in_pnt = &atmtensor_in_orig;
210 
211  const Tensor3& atmtensor_in = *atmtensor_in_pnt;
212 
213  // Resize output tensor
214  atmtensor_out.resize(
215  p_grid_new.nelem(), atmtensor_in.nrows(), atmtensor_in.ncols());
216 
217  ArrayOfGridPosPoly gp_p;
218  Matrix itw;
219 
220  Index ing_min, ing_max;
221 
222  AtmFieldPRegridHelper(ing_min,
223  ing_max,
224  gp_p,
225  itw,
226  p_grid_new,
227  p_grid_old,
228  interp_order,
229  verbosity);
230 
231  // Interpolate:
232  if ((ing_max - ing_min < 0) ||
233  (ing_max - ing_min + 1 != p_grid_new.nelem())) {
234  ostringstream os;
235  os << "New grid seems not to be sufficiently covered by old grid.\n";
236  throw runtime_error(os.str());
237  }
238 
239  for (Index i = 0; i < atmtensor_in.nrows(); i++)
240  for (Index j = 0; j < atmtensor_in.ncols(); j++)
241  interp(atmtensor_out(joker, i, j), itw, atmtensor_in(joker, i, j), gp_p);
242 }
243 
244 /* Workspace method: Doxygen documentation will be auto-generated */
245 void AtmFieldPRegrid( // WS Generic Output:
246  Tensor4& atmtensor_out,
247  // WS Generic Input:
248  const Tensor4& atmtensor_in_orig,
249  const Vector& p_grid_new,
250  const Vector& p_grid_old,
251  const Index& interp_order,
252  const Verbosity& verbosity) {
253  const Tensor4* atmtensor_in_pnt;
254  Tensor4 atmtensor_in_copy;
255 
256  if (&atmtensor_in_orig == &atmtensor_out) {
257  atmtensor_in_copy = atmtensor_in_orig;
258  atmtensor_in_pnt = &atmtensor_in_copy;
259  } else
260  atmtensor_in_pnt = &atmtensor_in_orig;
261 
262  const Tensor4& atmtensor_in = *atmtensor_in_pnt;
263 
264  // Resize output tensor
265  atmtensor_out.resize(atmtensor_in.nbooks(),
266  p_grid_new.nelem(),
267  atmtensor_in.nrows(),
268  atmtensor_in.ncols());
269 
270  ArrayOfGridPosPoly gp_p;
271  Matrix itw;
272 
273  Index ing_min, ing_max;
274 
275  AtmFieldPRegridHelper(ing_min,
276  ing_max,
277  gp_p,
278  itw,
279  p_grid_new,
280  p_grid_old,
281  interp_order,
282  verbosity);
283 
284  // Interpolate:
285  if ((ing_max - ing_min < 0) ||
286  (ing_max - ing_min + 1 != p_grid_new.nelem())) {
287  ostringstream os;
288  os << "New grid seems not to be sufficiently covered by old grid.\n";
289  throw runtime_error(os.str());
290  }
291 
292  for (Index b = 0; b < atmtensor_in.nbooks(); b++)
293  for (Index i = 0; i < atmtensor_in.nrows(); i++)
294  for (Index j = 0; j < atmtensor_in.ncols(); j++)
295  interp(atmtensor_out(b, joker, i, j),
296  itw,
297  atmtensor_in(b, joker, i, j),
298  gp_p);
299 }
300 
302 /*
303  Helper function for FieldFromGriddedField functions to ensure correct
304  dimension and values of latitude/longitude grids
305 
306  \param[in] lat_grid Latitude grid
307  \param[in] lon_grid Longitude grid
308  \param[in] ilat Latitude grid index in gfield
309  \param[in] ilon Longitude grid index in gfield
310  \param[in] gfield GriddedField
311  */
313  const Vector& lon_grid,
314  const Index ilat,
315  const Index ilon,
316  const GriddedField& gfield) {
317  chk_griddedfield_gridname(gfield, ilat, "Latitude");
318  chk_griddedfield_gridname(gfield, ilon, "Longitude");
319 
320  if (lon_grid.empty()) {
321  chk_size("gfield.lon_grid", gfield.get_numeric_grid(ilon), 1);
322 
323  if (lat_grid.empty())
324  chk_size("gfield.lat_grid", gfield.get_numeric_grid(ilat), 1);
325  else
326  chk_if_equal("lat_grid",
327  "gfield.lat_grid",
328  lat_grid,
329  gfield.get_numeric_grid(ilat));
330  } else {
331  chk_if_equal(
332  "lat_grid", "gfield.lat_grid", lat_grid, gfield.get_numeric_grid(ilat));
333  chk_if_equal(
334  "lon_grid", "gfield.lon_grid", lon_grid, gfield.get_numeric_grid(ilon));
335  }
336 }
337 
338 /* Workspace method: Doxygen documentation will be auto-generated */
339 void FieldFromGriddedField( // WS Generic Output:
340  Matrix& field_out,
341  // WS Input:
342  const Vector& p_grid _U_,
343  const Vector& lat_grid,
344  const Vector& lon_grid,
345  // WS Generic Input:
346  const GriddedField2& gfraw_in,
347  const Verbosity&) {
348  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 0, 1, gfraw_in);
349 
350  field_out = gfraw_in.data;
351 }
352 
353 /* Workspace method: Doxygen documentation will be auto-generated */
354 void FieldFromGriddedField( // WS Generic Output:
355  Tensor3& field_out,
356  // WS Input:
357  const Vector& p_grid,
358  const Vector& lat_grid,
359  const Vector& lon_grid,
360  // WS Generic Input:
361  const GriddedField3& gfraw_in,
362  const Verbosity&) {
363  chk_griddedfield_gridname(gfraw_in, 0, "Pressure");
364  chk_if_equal("p_grid", "gfield.p_grid", p_grid, gfraw_in.get_numeric_grid(0));
365 
366  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 1, 2, gfraw_in);
367 
368  field_out = gfraw_in.data;
369 }
370 
371 /* Workspace method: Doxygen documentation will be auto-generated */
372 void FieldFromGriddedField( // WS Generic Output:
373  Tensor4& field_out,
374  // WS Input:
375  const Vector& p_grid,
376  const Vector& lat_grid,
377  const Vector& lon_grid,
378  // WS Generic Input:
379  const GriddedField4& gfraw_in,
380  const Verbosity&) {
381  chk_griddedfield_gridname(gfraw_in, 1, "Pressure");
382  chk_if_equal("p_grid", "gfield.p_grid", p_grid, gfraw_in.get_numeric_grid(0));
383 
384  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 2, 3, gfraw_in);
385 
386  field_out = gfraw_in.data;
387 }
388 
389 /* Workspace method: Doxygen documentation will be auto-generated */
390 void FieldFromGriddedField( // WS Generic Output:
391  Tensor4& field_out,
392  // WS Input:
393  const Vector& p_grid,
394  const Vector& lat_grid,
395  const Vector& lon_grid,
396  // WS Generic Input:
397  const ArrayOfGriddedField3& gfraw_in,
398  const Verbosity& verbosity) {
399  if (!gfraw_in.nelem()) {
400  CREATE_OUT1;
401  out1 << " Warning: gfraw_in is empty, proceeding anyway\n";
402  field_out.resize(0, 0, 0, 0);
403  } else {
404  field_out.resize(gfraw_in.nelem(),
405  p_grid.nelem(),
406  gfraw_in[0].data.nrows(),
407  gfraw_in[0].data.ncols());
408  }
409 
410  for (Index i = 0; i < gfraw_in.nelem(); i++) {
411  chk_griddedfield_gridname(gfraw_in[i], 0, "Pressure");
412  chk_if_equal(
413  "p_grid", "gfield.p_grid", p_grid, gfraw_in[i].get_numeric_grid(0));
414 
416  lat_grid, lon_grid, 1, 2, gfraw_in[i]);
417 
418  field_out(i, joker, joker, joker) = gfraw_in[i].data;
419  }
420 }
421 
422 /* Workspace method: Doxygen documentation will be auto-generated */
423 void GriddedFieldLatLonExpand( // WS Generic Output:
424  GriddedField2& gfraw_out,
425  // WS Generic Input:
426  const GriddedField2& gfraw_in_orig,
427  const Verbosity&) {
428  const GriddedField2* gfraw_in_pnt;
429  GriddedField2 gfraw_in_copy;
430 
431  if (&gfraw_in_orig == &gfraw_out) {
432  gfraw_in_copy = gfraw_in_orig;
433  gfraw_in_pnt = &gfraw_in_copy;
434  } else
435  gfraw_in_pnt = &gfraw_in_orig;
436 
437  const GriddedField2& gfraw_in = *gfraw_in_pnt;
438 
439  chk_griddedfield_gridname(gfraw_in, 0, "Latitude");
440  chk_griddedfield_gridname(gfraw_in, 1, "Longitude");
441 
442  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
443  throw runtime_error(
444  "Can't expand data because number of Latitudes and Longitudes is greater than 1");
445 
446  gfraw_out.set_grid_name(0, "Latitude");
447  gfraw_out.set_grid_name(1, "Longitude");
448 
449  Vector v(2);
450  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1) {
451  v[0] = -90;
452  v[1] = 90;
453  gfraw_out.set_grid(0, v);
454  gfraw_out.resize(2, gfraw_in.data.ncols());
455 
456  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
457  gfraw_out.data(joker, j) = gfraw_in.data(0, j);
458  } else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1) {
459  v[0] = 0;
460  v[1] = 360;
461  gfraw_out.set_grid(1, v);
462  gfraw_out.resize(gfraw_in.data.nrows(), 2);
463 
464  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
465  gfraw_out.data(j, joker) = gfraw_in.data(j, 0);
466  } else {
467  v[0] = -90;
468  v[1] = 90;
469  gfraw_out.set_grid(0, v);
470  v[0] = 0;
471  v[1] = 360;
472  gfraw_out.set_grid(1, v);
473  gfraw_out.resize(2, 2);
474 
475  gfraw_out.data = gfraw_in.data(0, 0);
476  }
477 }
478 
479 /* Workspace method: Doxygen documentation will be auto-generated */
480 void GriddedFieldLatLonExpand( // WS Generic Output:
481  GriddedField3& gfraw_out,
482  // WS Generic Input:
483  const GriddedField3& gfraw_in_orig,
484  const Verbosity&) {
485  const GriddedField3* gfraw_in_pnt;
486  GriddedField3 gfraw_in_copy;
487 
488  if (&gfraw_in_orig == &gfraw_out) {
489  gfraw_in_copy = gfraw_in_orig;
490  gfraw_in_pnt = &gfraw_in_copy;
491  } else
492  gfraw_in_pnt = &gfraw_in_orig;
493 
494  const GriddedField3& gfraw_in = *gfraw_in_pnt;
495 
496  chk_griddedfield_gridname(gfraw_in, 1, "Latitude");
497  chk_griddedfield_gridname(gfraw_in, 2, "Longitude");
498 
499  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
500  throw runtime_error(
501  "Can't expand data because number of Latitudes and Longitudes is greater than 1");
502 
503  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
504  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
505  gfraw_out.set_grid_name(1, "Latitude");
506  gfraw_out.set_grid_name(2, "Longitude");
507 
508  Vector v(2);
509  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1) {
510  v[0] = -90;
511  v[1] = 90;
512  gfraw_out.set_grid(1, v);
513  gfraw_out.resize(gfraw_in.data.npages(), 2, gfraw_in.data.ncols());
514 
515  for (Index i = 0; i < gfraw_in.data.npages(); i++)
516  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
517  gfraw_out.data(i, joker, j) = gfraw_in.data(i, 0, j);
518  } else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1) {
519  v[0] = 0;
520  v[1] = 360;
521  gfraw_out.set_grid(2, v);
522  gfraw_out.resize(gfraw_in.data.npages(), gfraw_in.data.nrows(), 2);
523 
524  for (Index i = 0; i < gfraw_in.data.npages(); i++)
525  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
526  gfraw_out.data(i, j, joker) = gfraw_in.data(i, j, 0);
527  } else {
528  v[0] = -90;
529  v[1] = 90;
530  gfraw_out.set_grid(1, v);
531  v[0] = 0;
532  v[1] = 360;
533  gfraw_out.set_grid(2, v);
534  gfraw_out.resize(gfraw_in.data.npages(), 2, 2);
535 
536  for (Index i = 0; i < gfraw_in.data.npages(); i++)
537  gfraw_out.data(i, joker, joker) = gfraw_in.data(i, 0, 0);
538  }
539 }
540 
541 /* Workspace method: Doxygen documentation will be auto-generated */
542 void GriddedFieldLatLonExpand( // WS Generic Output:
543  GriddedField4& gfraw_out,
544  // WS Generic Input:
545  const GriddedField4& gfraw_in_orig,
546  const Verbosity&) {
547  const GriddedField4* gfraw_in_pnt;
548  GriddedField4 gfraw_in_copy;
549 
550  if (&gfraw_in_orig == &gfraw_out) {
551  gfraw_in_copy = gfraw_in_orig;
552  gfraw_in_pnt = &gfraw_in_copy;
553  } else
554  gfraw_in_pnt = &gfraw_in_orig;
555 
556  const GriddedField4& gfraw_in = *gfraw_in_pnt;
557 
558  chk_griddedfield_gridname(gfraw_in, 2, "Latitude");
559  chk_griddedfield_gridname(gfraw_in, 3, "Longitude");
560 
561  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
562  throw runtime_error(
563  "Can't expand data because number of Latitudes and Longitudes is greater than 1");
564 
565  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
566  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
567 
568  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
569  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
570 
571  gfraw_out.set_grid_name(2, "Latitude");
572  gfraw_out.set_grid_name(3, "Longitude");
573 
574  Vector v(2);
575  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1) {
576  v[0] = -90;
577  v[1] = 90;
578  gfraw_out.set_grid(2, v);
579  gfraw_out.resize(gfraw_in.data.nbooks(),
580  gfraw_in.data.npages(),
581  2,
582  gfraw_in.data.ncols());
583 
584  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
585  for (Index i = 0; i < gfraw_in.data.npages(); i++)
586  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
587  gfraw_out.data(k, i, joker, j) = gfraw_in.data(k, i, 0, j);
588  } else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1) {
589  v[0] = 0;
590  v[1] = 360;
591  gfraw_out.set_grid(3, v);
592  gfraw_out.resize(gfraw_in.data.nbooks(),
593  gfraw_in.data.npages(),
594  gfraw_in.data.nrows(),
595  2);
596 
597  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
598  for (Index i = 0; i < gfraw_in.data.npages(); i++)
599  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
600  gfraw_out.data(k, i, j, joker) = gfraw_in.data(k, i, j, 0);
601  } else {
602  v[0] = -90;
603  v[1] = 90;
604  gfraw_out.set_grid(2, v);
605  v[0] = 0;
606  v[1] = 360;
607  gfraw_out.set_grid(3, v);
608  gfraw_out.resize(gfraw_in.data.nbooks(), gfraw_in.data.npages(), 2, 2);
609 
610  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
611  for (Index i = 0; i < gfraw_in.data.npages(); i++)
612  gfraw_out.data(k, i, joker, joker) = gfraw_in.data(k, i, 0, 0);
613  }
614 }
615 
616 /* Workspace method: Doxygen documentation will be auto-generated */
617 void GriddedFieldLatLonExpand( // WS Generic Output:
618  ArrayOfGriddedField3& gfraw_out,
619  // WS Generic Input:
620  const ArrayOfGriddedField3& gfraw_in,
621  const Verbosity& verbosity) {
622  gfraw_out.resize(gfraw_in.nelem());
623 
624  for (Index i = 0; i < gfraw_in.nelem(); i++)
625  GriddedFieldLatLonExpand(gfraw_out[i], gfraw_in[i], verbosity);
626 }
627 
629 /*
630  This helper function is used by all GriddedFieldPRegrid WSMs to do the common
631  calculation of the grid positions and interpolation weights for pressures.
632 
633  \param[out] ing_min Index in the new grid with first value covered
634  by the old grid.
635  \param[out] ing_max Index in the new grid with last value covered
636  by the old grid.
637  \param[out] gp_p Pressure grid positions
638  \param[out] itw Interpolation weights
639  \param[in,out] gfraw_out Output GriddedField
640  \param[in] gfraw_in Input GriddedField
641  \param[in] p_grid_index Index of pressure grid
642  \param[in] p_grid New pressure grid
643  \param[in] interp_order Interpolation order
644  \param[in] zeropadding Allow zero padding
645  \param[in] verbosity Verbosity levels
646  */
648  Index& ing_max,
649  ArrayOfGridPosPoly& gp_p,
650  Matrix& itw,
651  GriddedField& gfraw_out,
652  const GriddedField& gfraw_in,
653  const Index p_grid_index,
654  ConstVectorView p_grid,
655  const Index& interp_order,
656  const Index& zeropadding,
657  const Verbosity& verbosity) {
658  CREATE_OUT2;
659 
660  chk_griddedfield_gridname(gfraw_in, p_grid_index, "Pressure");
661 
662  out2 << " Interpolation order: " << interp_order << "\n";
663 
664  const Vector& in_p_grid = gfraw_in.get_numeric_grid(p_grid_index);
665 
666  // Initialize output field. Set grids and copy grid names
667  gfraw_out.set_grid(p_grid_index, p_grid);
668  gfraw_out.set_grid_name(p_grid_index, gfraw_in.get_grid_name(p_grid_index));
669 
670  if (zeropadding) {
671  if (in_p_grid[0] < p_grid[p_grid.nelem() - 1] ||
672  in_p_grid[in_p_grid.nelem() - 1] > p_grid[0]) {
673  ing_min = 0;
674  ing_max = ing_min - 1;
675  } else
677  ing_max,
678  "Raw field to p_grid",
679  in_p_grid,
680  p_grid,
681  interp_order);
682  } else {
683  ing_min = 0;
684  ing_max = p_grid.nelem() - 1;
686  "Raw field to p_grid", in_p_grid, p_grid, interp_order);
687  }
688  Index nelem_in_range = ing_max - ing_min + 1;
689 
690  // Calculate grid positions:
691  if (nelem_in_range > 0) {
692  gp_p.resize(nelem_in_range);
694  gp_p, in_p_grid, p_grid[Range(ing_min, nelem_in_range)], interp_order);
695 
696  // Interpolation weights:
697  itw.resize(nelem_in_range, interp_order + 1);
698  interpweights(itw, gp_p);
699  }
700 }
701 
702 /* Workspace method: Doxygen documentation will be auto-generated */
703 void GriddedFieldPRegrid( // WS Generic Output:
704  GriddedField3& gfraw_out,
705  // WS Input:
706  const Vector& p_grid,
707  // WS Generic Input:
708  const GriddedField3& gfraw_in_orig,
709  const Index& interp_order,
710  const Index& zeropadding,
711  const Verbosity& verbosity) {
712  const GriddedField3* gfraw_in_pnt;
713  GriddedField3 gfraw_in_copy;
714 
715  if (&gfraw_in_orig == &gfraw_out) {
716  gfraw_in_copy = gfraw_in_orig;
717  gfraw_in_pnt = &gfraw_in_copy;
718  } else
719  gfraw_in_pnt = &gfraw_in_orig;
720 
721  const GriddedField3& gfraw_in = *gfraw_in_pnt;
722 
723  const Index p_grid_index = 0;
724 
725  // Resize output GriddedField and copy all non-latitude/longitude grids
726  gfraw_out.resize(
727  p_grid.nelem(), gfraw_in.data.nrows(), gfraw_in.data.ncols());
728  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
729  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
730  gfraw_out.set_grid(2, gfraw_in.get_numeric_grid(2));
731  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
732 
733  ArrayOfGridPosPoly gp_p;
734  Matrix itw;
735 
736  Index ing_min, ing_max;
737 
739  ing_max,
740  gp_p,
741  itw,
742  gfraw_out,
743  gfraw_in,
744  p_grid_index,
745  p_grid,
746  interp_order,
747  zeropadding,
748  verbosity);
749 
750  // Interpolate:
751  if (ing_max - ing_min < 0)
752  gfraw_out.data = 0.;
753  else if (ing_max - ing_min + 1 != p_grid.nelem()) {
754  gfraw_out.data = 0.;
755  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
756  for (Index j = 0; j < gfraw_in.data.ncols(); j++) {
757  // chk_interpolation_grids_loose_check_data(ing_min, ing_max,
758  // "Raw field to p_grid",
759  // gfraw_in.get_numeric_grid(p_grid_index),
760  // p_grid, gfraw_in.data(joker, i, j));
761  interp(gfraw_out.data(Range(ing_min, ing_max - ing_min + 1), i, j),
762  itw,
763  gfraw_in.data(joker, i, j),
764  gp_p);
765  }
766  } else
767  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
768  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
769  interp(
770  gfraw_out.data(joker, i, j), itw, gfraw_in.data(joker, i, j), gp_p);
771 }
772 
773 /* Workspace method: Doxygen documentation will be auto-generated */
774 void GriddedFieldPRegrid( // WS Generic Output:
775  GriddedField4& gfraw_out,
776  // WS Input:
777  const Vector& p_grid,
778  // WS Generic Input:
779  const GriddedField4& gfraw_in_orig,
780  const Index& interp_order,
781  const Index& zeropadding,
782  const Verbosity& verbosity) {
783  const GriddedField4* gfraw_in_pnt;
784  GriddedField4 gfraw_in_copy;
785 
786  if (&gfraw_in_orig == &gfraw_out) {
787  gfraw_in_copy = gfraw_in_orig;
788  gfraw_in_pnt = &gfraw_in_copy;
789  } else
790  gfraw_in_pnt = &gfraw_in_orig;
791 
792  const GriddedField4& gfraw_in = *gfraw_in_pnt;
793 
794  const Index p_grid_index = 1;
795 
796  // Resize output GriddedField and copy all non-latitude/longitude grids
797  gfraw_out.resize(gfraw_in.data.nbooks(),
798  p_grid.nelem(),
799  gfraw_in.data.nrows(),
800  gfraw_in.data.ncols());
801  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
802  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
803  gfraw_out.set_grid(2, gfraw_in.get_numeric_grid(2));
804  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
805  gfraw_out.set_grid(3, gfraw_in.get_numeric_grid(3));
806  gfraw_out.set_grid_name(3, gfraw_in.get_grid_name(3));
807 
808  ArrayOfGridPosPoly gp_p;
809  Matrix itw;
810 
811  Index ing_min, ing_max;
812 
814  ing_max,
815  gp_p,
816  itw,
817  gfraw_out,
818  gfraw_in,
819  p_grid_index,
820  p_grid,
821  interp_order,
822  zeropadding,
823  verbosity);
824 
825  // Interpolate:
826  if (ing_max - ing_min < 0)
827  gfraw_out.data = 0.;
828  else if (ing_max - ing_min + 1 != p_grid.nelem()) {
829  gfraw_out.data = 0.;
830  for (Index b = 0; b < gfraw_in.data.nbooks(); b++)
831  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
832  for (Index j = 0; j < gfraw_in.data.ncols(); j++) {
833  // chk_interpolation_grids_loose_check_data(ing_min, ing_max,
834  // "Raw field to p_grid",
835  // gfraw_in.get_numeric_grid(p_grid_index),
836  // p_grid, gfraw_in.data(b, joker, i, j));
837  interp(gfraw_out.data(b, Range(ing_min, ing_max - ing_min), i, j),
838  itw,
839  gfraw_in.data(b, joker, i, j),
840  gp_p);
841  }
842  } else
843  for (Index b = 0; b < gfraw_in.data.nbooks(); b++)
844  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
845  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
846  interp(gfraw_out.data(b, joker, i, j),
847  itw,
848  gfraw_in.data(b, joker, i, j),
849  gp_p);
850 }
851 
852 /* Workspace method: Doxygen documentation will be auto-generated */
853 void GriddedFieldPRegrid( // WS Generic Output:
854  ArrayOfGriddedField3& agfraw_out,
855  // WS Input:
856  const Vector& p_grid,
857  // WS Generic Input:
858  const ArrayOfGriddedField3& agfraw_in,
859  const Index& interp_order,
860  const Index& zeropadding,
861  const Verbosity& verbosity) {
862  agfraw_out.resize(agfraw_in.nelem());
863 
864  for (Index i = 0; i < agfraw_in.nelem(); i++) {
865  GriddedFieldPRegrid(agfraw_out[i],
866  p_grid,
867  agfraw_in[i],
868  interp_order,
869  zeropadding,
870  verbosity);
871  }
872 }
873 
875 /*
876  This helper function is used by all GriddedFieldLatLonRegrid WSMs to do the common
877  calculation of the grid positions and interpolation weights for latitudes and longitudes.
878 
879  \param[out] gp_lat Latitude grid positions
880  \param[out] gp_lon Longitude grid positions
881  \param[out] itw Interpolation weights
882  \param[in,out] gfraw_out Output GriddedField
883  \param[in] gfraw_in Input GriddedField
884  \param[in] lat_grid_index Index of latitude grid
885  \param[in] lon_grid_index Index of longitude grid
886  \param[in] lat_true New latitude grid
887  \param[in] lon_true New longitude grid
888  \param[in] interp_order Interpolation order
889  \param[in] verbosity Verbosity levels
890  */
892  ArrayOfGridPosPoly& gp_lon,
893  Tensor3& itw,
894  GriddedField& gfraw_out,
895  const GriddedField& gfraw_in,
896  const Index lat_grid_index,
897  const Index lon_grid_index,
898  ConstVectorView lat_true,
899  ConstVectorView lon_true,
900  const Index& interp_order,
901  const Verbosity& verbosity) {
902  CREATE_OUT2;
903 
904  if (!lat_true.nelem())
905  throw runtime_error("The new latitude grid is not allowed to be empty.");
906  if (!lon_true.nelem())
907  throw runtime_error("The new longitude grid is not allowed to be empty.");
908 
909  chk_griddedfield_gridname(gfraw_in, lat_grid_index, "Latitude");
910  chk_griddedfield_gridname(gfraw_in, lon_grid_index, "Longitude");
911  if (gfraw_in.get_grid_size(lat_grid_index) == 1 ||
912  gfraw_in.get_grid_size(lon_grid_index) == 1)
913  throw runtime_error("Raw data has to be true 3D data (nlat>1 and nlon>1).");
914 
915  out2 << " Interpolation order: " << interp_order << "\n";
916 
917  const Vector& in_lat_grid = gfraw_in.get_numeric_grid(lat_grid_index);
918  const Vector& in_lon_grid = gfraw_in.get_numeric_grid(lon_grid_index);
919 
920  // Initialize output field. Set grids and copy grid names
921  gfraw_out.set_grid(lat_grid_index, lat_true);
922  gfraw_out.set_grid_name(lat_grid_index,
923  gfraw_in.get_grid_name(lat_grid_index));
924  gfraw_out.set_grid(lon_grid_index, lon_true);
925  gfraw_out.set_grid_name(lon_grid_index,
926  gfraw_in.get_grid_name(lon_grid_index));
927 
929  "Raw field to lat_grid, 3D case", in_lat_grid, lat_true, interp_order);
930 
931  // Calculate grid positions:
932  gp_lat.resize(lat_true.nelem());
933  gp_lon.resize(lon_true.nelem());
934  gridpos_poly(gp_lat, in_lat_grid, lat_true, interp_order);
935 
936  gridpos_poly_longitudinal("Raw field to lon_grid, 3D case",
937  gp_lon,
938  in_lon_grid,
939  lon_true,
940  interp_order);
941 
942  // Interpolation weights:
943  itw.resize(lat_true.nelem(),
944  lon_true.nelem(),
945  (interp_order + 1) * (interp_order + 1));
946  interpweights(itw, gp_lat, gp_lon);
947 }
948 
949 /* Workspace method: Doxygen documentation will be auto-generated */
950 void GriddedFieldLatLonRegrid( // WS Generic Output:
951  GriddedField2& gfraw_out,
952  // WS Input:
953  const Vector& lat_true,
954  const Vector& lon_true,
955  // WS Generic Input:
956  const GriddedField2& gfraw_in_orig,
957  const Index& interp_order,
958  const Verbosity& verbosity) {
959  if (!lat_true.nelem())
960  throw runtime_error("The new latitude grid is not allowed to be empty.");
961  if (!lon_true.nelem())
962  throw runtime_error("The new longitude grid is not allowed to be empty.");
963 
964  const GriddedField2* gfraw_in_pnt;
965  GriddedField2 gfraw_in_copy;
966 
967  if (&gfraw_in_orig == &gfraw_out) {
968  gfraw_in_copy = gfraw_in_orig;
969  gfraw_in_pnt = &gfraw_in_copy;
970  } else
971  gfraw_in_pnt = &gfraw_in_orig;
972 
973  const GriddedField2& gfraw_in = *gfraw_in_pnt;
974 
975  const Index lat_grid_index = 0;
976  const Index lon_grid_index = 1;
977 
978  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
979  gfraw_in.get_grid_size(lon_grid_index) < 2) {
980  ostringstream os;
981  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
982  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
983  throw runtime_error(os.str());
984  }
985 
986  // Resize output GriddedField
987  gfraw_out.resize(lat_true.nelem(), lon_true.nelem());
988 
989  ArrayOfGridPosPoly gp_lat;
990  ArrayOfGridPosPoly gp_lon;
991  Tensor3 itw;
992 
993  // If lon grid is cyclic, the data values at 0 and 360 must match
994  const Vector& in_lat_grid =
995  gfraw_in.get_numeric_grid(lat_grid_index);
996  const Vector& in_lon_grid =
997  gfraw_in.get_numeric_grid(lon_grid_index);
998 
999  if (is_lon_cyclic(in_lon_grid)) {
1000  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++) {
1001  if (!is_same_within_epsilon(gfraw_in.data(lat, 0),
1002  gfraw_in.data(lat, in_lon_grid.nelem() - 1),
1003  EPSILON_LON_CYCLIC)) {
1004  ostringstream os;
1005  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
1006  << "Mismatch at latitude index : " << lat << " ("
1007  << in_lat_grid[lat] << " degrees)\n"
1008  << "Value at 0 degrees longitude : " << gfraw_in.data(lat, 0)
1009  << "\n"
1010  << "Value at 360 degrees longitude: "
1011  << gfraw_in.data(lat, in_lon_grid.nelem() - 1) << "\n"
1012  << "Difference : "
1013  << gfraw_in.data(lat, in_lon_grid.nelem() - 1) -
1014  gfraw_in.data(lat, 0)
1015  << "\n"
1016  << "Allowed difference : " << EPSILON_LON_CYCLIC;
1017  throw runtime_error(os.str());
1018  }
1019  }
1020  }
1021 
1023  gp_lon,
1024  itw,
1025  gfraw_out,
1026  gfraw_in,
1027  lat_grid_index,
1028  lon_grid_index,
1029  lat_true,
1030  lon_true,
1031  interp_order,
1032  verbosity);
1033 
1034  // Interpolate:
1035  interp(gfraw_out.data, itw, gfraw_in.data, gp_lat, gp_lon);
1036 }
1037 
1038 /* Workspace method: Doxygen documentation will be auto-generated */
1039 void GriddedFieldLatLonRegrid( // WS Generic Output:
1040  GriddedField3& gfraw_out,
1041  // WS Input:
1042  const Vector& lat_true,
1043  const Vector& lon_true,
1044  // WS Generic Input:
1045  const GriddedField3& gfraw_in_orig,
1046  const Index& interp_order,
1047  const Verbosity& verbosity) {
1048  if (!lat_true.nelem())
1049  throw runtime_error("The new latitude grid is not allowed to be empty.");
1050  if (!lon_true.nelem())
1051  throw runtime_error("The new longitude grid is not allowed to be empty.");
1052 
1053  const GriddedField3* gfraw_in_pnt;
1054  GriddedField3 gfraw_in_copy;
1055 
1056  if (&gfraw_in_orig == &gfraw_out) {
1057  gfraw_in_copy = gfraw_in_orig;
1058  gfraw_in_pnt = &gfraw_in_copy;
1059  } else
1060  gfraw_in_pnt = &gfraw_in_orig;
1061 
1062  const GriddedField3& gfraw_in = *gfraw_in_pnt;
1063 
1064  const Index lat_grid_index = 1;
1065  const Index lon_grid_index = 2;
1066 
1067  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
1068  gfraw_in.get_grid_size(lon_grid_index) < 2) {
1069  ostringstream os;
1070  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
1071  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
1072  throw runtime_error(os.str());
1073  }
1074 
1075  // Resize output GriddedField and copy all non-latitude/longitude grids
1076  gfraw_out.resize(gfraw_in.data.npages(), lat_true.nelem(), lon_true.nelem());
1077  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
1078  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
1079 
1080  ArrayOfGridPosPoly gp_lat;
1081  ArrayOfGridPosPoly gp_lon;
1082  Tensor3 itw;
1083 
1084  // If lon grid is cyclic, the data values at 0 and 360 must match
1085  const Vector& in_grid0 = gfraw_in.get_numeric_grid(0);
1086  const Vector& in_lat_grid =
1087  gfraw_in.get_numeric_grid(lat_grid_index);
1088  const Vector& in_lon_grid =
1089  gfraw_in.get_numeric_grid(lon_grid_index);
1090 
1091  if (is_lon_cyclic(in_lon_grid)) {
1092  for (Index g0 = 0; g0 < in_grid0.nelem(); g0++)
1093  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++) {
1095  gfraw_in.data(g0, lat, 0),
1096  gfraw_in.data(g0, lat, in_lon_grid.nelem() - 1),
1097  EPSILON_LON_CYCLIC)) {
1098  ostringstream os;
1099  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
1100  << "Mismatch at 1st grid index : " << g0 << " (" << in_grid0[g0]
1101  << ")\n"
1102  << " at latitude index : " << lat << " ("
1103  << in_lat_grid[lat] << " degrees)\n"
1104  << "Value at 0 degrees longitude : " << gfraw_in.data(g0, lat, 0)
1105  << "\n"
1106  << "Value at 360 degrees longitude: "
1107  << gfraw_in.data(g0, lat, in_lon_grid.nelem() - 1) << "\n"
1108  << "Difference : "
1109  << gfraw_in.data(g0, lat, in_lon_grid.nelem() - 1) -
1110  gfraw_in.data(g0, lat, 0)
1111  << "\n"
1112  << "Allowed difference : " << EPSILON_LON_CYCLIC;
1113  throw runtime_error(os.str());
1114  }
1115  }
1116  }
1117 
1119  gp_lon,
1120  itw,
1121  gfraw_out,
1122  gfraw_in,
1123  lat_grid_index,
1124  lon_grid_index,
1125  lat_true,
1126  lon_true,
1127  interp_order,
1128  verbosity);
1129 
1130  // Interpolate:
1131  for (Index i = 0; i < gfraw_in.data.npages(); i++)
1132  interp(gfraw_out.data(i, joker, joker),
1133  itw,
1134  gfraw_in.data(i, joker, joker),
1135  gp_lat,
1136  gp_lon);
1137 }
1138 
1139 /* Workspace method: Doxygen documentation will be auto-generated */
1140 void GriddedFieldLatLonRegrid( // WS Generic Output:
1141  GriddedField4& gfraw_out,
1142  // WS Input:
1143  const Vector& lat_true,
1144  const Vector& lon_true,
1145  // WS Generic Input:
1146  const GriddedField4& gfraw_in_orig,
1147  const Index& interp_order,
1148  const Verbosity& verbosity) {
1149  if (!lat_true.nelem())
1150  throw runtime_error("The new latitude grid is not allowed to be empty.");
1151  if (!lon_true.nelem())
1152  throw runtime_error("The new longitude grid is not allowed to be empty.");
1153 
1154  const GriddedField4* gfraw_in_pnt;
1155  GriddedField4 gfraw_in_copy;
1156 
1157  if (&gfraw_in_orig == &gfraw_out) {
1158  gfraw_in_copy = gfraw_in_orig;
1159  gfraw_in_pnt = &gfraw_in_copy;
1160  } else
1161  gfraw_in_pnt = &gfraw_in_orig;
1162 
1163  const GriddedField4& gfraw_in = *gfraw_in_pnt;
1164 
1165  const Index lat_grid_index = 2;
1166  const Index lon_grid_index = 3;
1167 
1168  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
1169  gfraw_in.get_grid_size(lon_grid_index) < 2) {
1170  ostringstream os;
1171  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
1172  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
1173  throw runtime_error(os.str());
1174  }
1175 
1176  // Resize output GriddedField and copy all non-latitude/longitude grids
1177  gfraw_out.resize(gfraw_in.data.nbooks(),
1178  gfraw_in.data.npages(),
1179  lat_true.nelem(),
1180  lon_true.nelem());
1181  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
1182  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
1183  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
1184  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
1185 
1186  ArrayOfGridPosPoly gp_lat;
1187  ArrayOfGridPosPoly gp_lon;
1188  Tensor3 itw;
1189 
1191  gp_lon,
1192  itw,
1193  gfraw_out,
1194  gfraw_in,
1195  lat_grid_index,
1196  lon_grid_index,
1197  lat_true,
1198  lon_true,
1199  interp_order,
1200  verbosity);
1201 
1202  // If lon grid is cyclic, the data values at 0 and 360 must match
1203  const Vector& in_grid0 = gfraw_in.get_numeric_grid(0);
1204  const Vector& in_grid1 = gfraw_in.get_numeric_grid(1);
1205  const Vector& in_lat_grid =
1206  gfraw_in.get_numeric_grid(lat_grid_index);
1207  const Vector& in_lon_grid =
1208  gfraw_in.get_numeric_grid(lon_grid_index);
1209 
1210  if (is_lon_cyclic(in_lon_grid)) {
1211  for (Index g0 = 0; g0 < in_grid0.nelem(); g0++)
1212  for (Index g1 = 0; g1 < in_grid1.nelem(); g1++)
1213  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++) {
1215  gfraw_in.data(g0, g1, lat, 0),
1216  gfraw_in.data(g0, g1, lat, in_lon_grid.nelem() - 1),
1217  EPSILON_LON_CYCLIC)) {
1218  ostringstream os;
1219  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
1220  << "Mismatch at 1st grid index : " << g0 << " ("
1221  << in_grid0[g0] << ")\n"
1222  << " at 2nd grid index : " << g1 << " ("
1223  << in_grid1[g1] << ")\n"
1224  << " at latitude index : " << lat << " ("
1225  << in_lat_grid[lat] << " degrees)\n"
1226  << "Value at 0 degrees longitude : "
1227  << gfraw_in.data(g0, g1, lat, 0) << "\n"
1228  << "Value at 360 degrees longitude: "
1229  << gfraw_in.data(g0, g1, lat, in_lon_grid.nelem() - 1) << "\n"
1230  << "Difference : "
1231  << gfraw_in.data(g0, g1, lat, in_lon_grid.nelem() - 1) -
1232  gfraw_in.data(g0, g1, lat, 0)
1233  << "\n"
1234  << "Allowed difference : " << EPSILON_LON_CYCLIC;
1235  throw runtime_error(os.str());
1236  }
1237  }
1238  }
1239 
1240  // Interpolate:
1241  for (Index i = 0; i < gfraw_in.data.nbooks(); i++)
1242  for (Index j = 0; j < gfraw_in.data.npages(); j++)
1243  interp(gfraw_out.data(i, j, joker, joker),
1244  itw,
1245  gfraw_in.data(i, j, joker, joker),
1246  gp_lat,
1247  gp_lon);
1248 }
1249 
1250 /* Workspace method: Doxygen documentation will be auto-generated */
1251 void GriddedFieldLatLonRegrid( // WS Generic Output:
1252  ArrayOfGriddedField3& agfraw_out,
1253  // WS Input:
1254  const Vector& lat_true,
1255  const Vector& lon_true,
1256  // WS Generic Input:
1257  const ArrayOfGriddedField3& agfraw_in,
1258  const Index& interp_order,
1259  const Verbosity& verbosity) {
1260  agfraw_out.resize(agfraw_in.nelem());
1261 
1262  for (Index i = 0; i < agfraw_in.nelem(); i++) {
1263  GriddedFieldLatLonRegrid(agfraw_out[i],
1264  lat_true,
1265  lon_true,
1266  agfraw_in[i],
1267  interp_order,
1268  verbosity);
1269  }
1270 }
1271 
1273 /*
1274  This helper function is used by GriddedFieldZToPRegrid WSM to do the common
1275  calculation of the grid positions and interpolation weights for latitudes and longitudes.
1276 
1277  \param[out] ing_min Index in the new grid with first value covered
1278  by the old grid.
1279  \param[out] ing_max Index in the new grid with last value covered
1280  by the old grid.
1281  \param[out] gp_p Altitude grid positions
1282  \param[out] itw Interpolation weights
1283  \param[in] gfraw_in Input GriddedField
1284  \param[in] z_grid_index Index of altitude grid
1285  \param[in] z_grid New z_grid grid
1286  \param[in] interp_order Interpolation order
1287  \param[in] zeropadding Allow zero padding
1288  \param[in] verbosity Verbosity levels
1289  */
1291  Index& ing_max,
1292  ArrayOfGridPosPoly& gp_p,
1293  Matrix& itw,
1294  const GriddedField& gfraw_in,
1295  const Index z_grid_index,
1296  ConstVectorView z_grid,
1297  const Index& interp_order,
1298  const Index& zeropadding,
1299  const Verbosity& verbosity) {
1300  CREATE_OUT2;
1301 
1302  chk_griddedfield_gridname(gfraw_in, z_grid_index, "Altitude");
1303 
1304  out2 << " Interpolation order: " << interp_order << "\n";
1305 
1306  const Vector& in_z_grid = gfraw_in.get_numeric_grid(z_grid_index);
1307 
1308  if (zeropadding) {
1309  if (in_z_grid[0] > z_grid[z_grid.nelem() - 1] ||
1310  in_z_grid[in_z_grid.nelem() - 1] < z_grid[0]) {
1311  ing_min = 0;
1312  ing_max = ing_min - 1;
1313  } else
1315  ing_max,
1316  "Raw field to z_grid",
1317  in_z_grid,
1318  z_grid,
1319  interp_order);
1320  } else {
1321  ing_min = 0;
1322  ing_max = z_grid.nelem() - 1;
1324  "Raw field to p_grid", in_z_grid, z_grid, interp_order);
1325  }
1326 
1327  Index nelem_in_range = ing_max - ing_min + 1;
1328 
1329  // Calculate grid positions:
1330  if (nelem_in_range > 0) {
1331  gp_p.resize(nelem_in_range);
1332  gridpos_poly(
1333  gp_p, in_z_grid, z_grid[Range(ing_min, nelem_in_range)], interp_order);
1334 
1335  // Interpolation weights:
1336  itw.resize(nelem_in_range, interp_order + 1);
1337  interpweights(itw, gp_p);
1338  }
1339 }
1340 
1341 /* Workspace method: Doxygen documentation will be auto-generated */
1342 void GriddedFieldZToPRegrid( // WS Generic Output:
1343  GriddedField3& gfraw_out, //grid in P
1344  // WS Input:
1345  const Vector& p_grid,
1346  const Vector& lat_grid,
1347  const Vector& lon_grid,
1348  const Tensor3& z_field,
1349  // WS Generic Input:
1350  const GriddedField3& gfraw_in_orig, //grid in Z
1351  const Index& interp_order, // Only linear interpolation allowed
1352  const Index& zeropadding,
1353  const Verbosity& verbosity) {
1354  // z_field must be of the same size as its grids
1355  if (!((z_field.npages() == p_grid.nelem() &&
1356  z_field.nrows() == lat_grid.nelem()) &&
1357  z_field.ncols() == lon_grid.nelem()))
1358  throw std::runtime_error(
1359  "*z_field* must be of the same size as *p_grid*, *lat_grid*, and *lon_grid* in *GriddedFieldZToPRegrid*.");
1360 
1361  // Must name the dimension "Altitude" to ensure user is aware of what they are doing.
1362  chk_griddedfield_gridname(gfraw_in_orig, 0, "Altitude");
1363 
1364  // Each lat and lon grid must be identical between z_field and in field
1365  const Vector& lat_in = gfraw_in_orig.get_numeric_grid(1);
1366  const Vector& lon_in = gfraw_in_orig.get_numeric_grid(2);
1367 
1368  if (lat_grid.nelem() != lat_in.nelem() || lon_grid.nelem() != lon_in.nelem())
1369  throw std::runtime_error(
1370  "Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude to be on the same grid as *z_field*.");
1371  else {
1372  for (Index ii = 0; ii < lat_grid.nelem(); ii++)
1373  if (lat_grid[ii] != lat_in[ii])
1374  throw std::runtime_error(
1375  "Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude of the gridded field to be the same as for *z_field*.");
1376  for (Index ii = 0; ii < lon_grid.nelem(); ii++)
1377  if (lon_grid[ii] != lon_in[ii])
1378  throw std::runtime_error(
1379  "Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude of the gridded field to be the same as for *z_field*.");
1380  }
1381 
1382  // Pointer in case output is input variable (same memory allocated)
1383  const GriddedField3* gfraw_in_pnt;
1384  GriddedField3 gfraw_in_copy;
1385 
1386  if (&gfraw_in_orig == &gfraw_out) {
1387  gfraw_in_copy = gfraw_in_orig;
1388  gfraw_in_pnt = &gfraw_in_copy;
1389  } else
1390  gfraw_in_pnt = &gfraw_in_orig;
1391 
1392  // Now output and input are separate variables (not allocating the same memory)
1393  const GriddedField3& gfraw_in = *gfraw_in_pnt;
1394 
1395  // Right size and order
1396  gfraw_out.resize(p_grid.nelem(), lat_grid.nelem(), lon_grid.nelem());
1397  gfraw_out.set_grid(0, p_grid);
1398  gfraw_out.set_grid_name(0, "Pressure");
1399  gfraw_out.set_grid(1, lat_grid);
1400  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
1401  gfraw_out.set_grid(2, lon_grid);
1402  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
1403  gfraw_out.data = 0.;
1404 
1405  ArrayOfGridPosPoly gp_p;
1406  Matrix itw;
1407 
1408  Index ing_min, ing_max;
1409 
1410  for (Index lat_index = 0; lat_index < lat_grid.nelem(); lat_index++) {
1411  for (Index lon_index = 0; lon_index < lon_grid.nelem(); lon_index++) {
1412  const Vector z_out = z_field(joker, lat_index, lon_index);
1413 
1415  ing_max,
1416  gp_p,
1417  itw,
1418  gfraw_in,
1419  0,
1420  z_out,
1421  interp_order,
1422  zeropadding,
1423  verbosity);
1424 
1425  if (ing_max - ing_min >= 0) {
1426  Range r = joker;
1427 
1428  if (ing_max - ing_min + 1 != z_out.nelem()) {
1429  r = Range(ing_min, ing_max - ing_min + 1);
1430  }
1431 
1432  interp(gfraw_out.data(r, lat_index, lon_index),
1433  itw,
1434  gfraw_in.data(joker, lat_index, lon_index),
1435  gp_p);
1436  }
1437  }
1438  }
1439 }
1440 
1441 // Workspace method, doxygen header will be auto-generated.
1442 // 2007-07-25 Stefan Buehler
1443 void atm_fields_compactFromMatrix( // WS Output:
1444  GriddedField4& af, // atm_fields_compact
1445  // WS Input:
1446  const Index& atmosphere_dim,
1447  // WS Generic Input:
1448  const Matrix& im,
1449  // Control Parameters:
1450  const ArrayOfString& field_names,
1451  const Verbosity&) {
1452  if (1 != atmosphere_dim) {
1453  ostringstream os;
1454  os << "Atmospheric dimension must be 1.";
1455  throw runtime_error(os.str());
1456  }
1457 
1458  const Index np = im.nrows(); // Number of pressure levels.
1459  const Index nf = im.ncols() - 1; // Total number of fields.
1460  ArrayOfIndex f_1; // indices of non-ignored fields
1461  // All fields called "ignore" will be ignored.
1462  String fn_upper; // Temporary variable to hold upper case field_names.
1463 
1464  if (field_names.nelem() != nf) {
1465  ostringstream os;
1466  os << "Cannot extract fields from Matrix.\n"
1467  << "*field_names* must have one element less than there are\n"
1468  << "matrix columns.";
1469  throw runtime_error(os.str());
1470  }
1471 
1472  // Remove additional fields from the field_names. All fields that
1473  // are flagged by 'ignore' in the field names, small or large letters,
1474  // are removed.
1475  for (Index f = 0; f < field_names.nelem(); f++) {
1476  fn_upper = field_names[f];
1477  fn_upper.toupper();
1478  //cout << "fieldname[" << f << "]: " << fn_upper;
1479  if (fn_upper != "IGNORE") {
1480  f_1.push_back(f);
1481  //cout << " put as element " << f_1.size()-1 << " into selection\n";
1482  }
1483  }
1484 
1485  // Copy required field_names to a new variable called field_names_1
1486  Index nf_1 = f_1.size(); // Number of required fields.
1487  ArrayOfString field_names_1(nf_1);
1488  for (Index f = 0; f < nf_1; f++) {
1489  field_names_1[f] = field_names[f_1[f]];
1490  //cout << "new fieldname[" << f << "] (formerly [" << f_1[f] << "]): "
1491  // << field_names_1[f] << "\n";
1492  }
1493 
1494  // out3 << "Copying *" << im_name << "* to *atm_fields_compact*.\n";
1495 
1496  af.set_grid(GFIELD4_FIELD_NAMES, field_names_1);
1497 
1498  af.set_grid(GFIELD4_P_GRID, im(Range(joker), 0));
1499 
1502 
1503  af.resize(nf_1, np, 1, 1); // Resize it according to the required fields
1504  for (Index f = 0; f < nf_1; f++)
1505  af.data(f, Range(joker), 0, 0) = im(Range(joker), f_1[f] + 1);
1506 }
1507 
1508 // Workspace method, doxygen header is auto-generated.
1509 // 2007-07-31 Stefan Buehler
1510 // 2011-05-04 Adapted by Gerrit Holl
1512  GriddedField4& af,
1513  // Control Parameters:
1514  const String& name,
1515  const Numeric& value,
1516  const Index& prepend,
1517  const ArrayOfString& condensibles,
1518  const Verbosity& verbosity) {
1519  Index nf; // Will hold new size
1520 
1521  // Add book
1522  atm_fields_compactExpand(af, nf, name, prepend, verbosity);
1523 
1524  if (condensibles.nelem()) {
1525  const Tensor4& vmrs = af.data;
1526  const ArrayOfString& species = af.get_string_grid(GFIELD4_FIELD_NAMES);
1527  Tensor3 condensible_sum(vmrs.npages(), vmrs.nrows(), vmrs.ncols(), 1.);
1528  for (Index c = 0; c < condensibles.nelem(); c++) {
1529  bool species_found = false;
1530  for (Index i = 0; !species_found && i < species.nelem(); i++) {
1531  if (species[i] == condensibles[c]) {
1532  condensible_sum -= vmrs(i, joker, joker, joker);
1533  species_found = true;
1534  }
1535  }
1536  if (!species_found) {
1537  std::ostringstream os;
1538  os << "Condensible species \"" << condensibles[c] << "\" not found ";
1539  os << "in input data.";
1540  throw runtime_error(os.str());
1541  }
1542  }
1543  condensible_sum *= value;
1544  // Add the constant value:
1545  if (prepend)
1546  af.data(0, joker, joker, joker) = condensible_sum;
1547  else
1548  af.data(nf - 1, joker, joker, joker) = condensible_sum;
1549  } else {
1550  // Add the constant value:
1551  if (prepend)
1552  af.data(0, joker, joker, joker) = value;
1553  else
1554  af.data(nf - 1, joker, joker, joker) = value;
1555  }
1556 }
1557 
1558 // Workspace method, doxygen header is auto-generated
1559 // 2011-05-02 Gerrit Holl
1560 void atm_fields_compactAddSpecies( // WS Output:
1561  GriddedField4& atm_fields_compact,
1562  // WS Generic Input:
1563  const String& name,
1564  const GriddedField3& species,
1565  const Index& prepend,
1566  const Verbosity& verbosity) {
1567  assert(atm_fields_compact.checksize());
1568  assert(species.checksize());
1569 
1570  ConstVectorView af_p_grid =
1571  atm_fields_compact.get_numeric_grid(GFIELD4_P_GRID);
1572  ConstVectorView af_lat_grid =
1573  atm_fields_compact.get_numeric_grid(GFIELD4_LAT_GRID);
1574  ConstVectorView af_lon_grid =
1575  atm_fields_compact.get_numeric_grid(GFIELD4_LON_GRID);
1576  ConstVectorView sp_p_grid = species.get_numeric_grid(GFIELD3_P_GRID);
1577  ConstVectorView sp_lat_grid = species.get_numeric_grid(GFIELD3_LAT_GRID);
1578  ConstVectorView sp_lon_grid = species.get_numeric_grid(GFIELD3_LON_GRID);
1579 
1580  Index new_n_fields; // To be set in next line
1582  atm_fields_compact, new_n_fields, name, prepend, verbosity);
1583 
1584  const Index insert_pos = (prepend) ? 0 : new_n_fields - 1;
1585 
1586  // Interpolate species to atm_fields_compact
1587 
1588  // Common for all dim
1590  "species p_grid to atm_fields_compact p_grid", sp_p_grid, af_p_grid);
1591  ArrayOfGridPos p_gridpos(af_p_grid.nelem());
1592  // gridpos(p_gridpos, sp_p_grid, af_p_grid);
1593  p2gridpos(p_gridpos, sp_p_grid, af_p_grid);
1594 
1595  if (sp_lat_grid.nelem() > 1) {
1596  // Common for all dim>=2
1597  chk_interpolation_grids("species lat_grid to atm_fields_compact lat_grid",
1598  sp_lat_grid,
1599  af_lat_grid);
1600  ArrayOfGridPos lat_gridpos(af_lat_grid.nelem());
1601  gridpos(lat_gridpos, sp_lat_grid, af_lat_grid);
1602 
1603  if (sp_lon_grid.nelem() > 1) { // 3D-case
1604  chk_interpolation_grids("species lon_grid to atm_fields_compact lon_grid",
1605  sp_lon_grid,
1606  af_lon_grid);
1607  ArrayOfGridPos lon_gridpos(af_lon_grid.nelem());
1608  gridpos(lon_gridpos, sp_lon_grid, af_lon_grid);
1609 
1610  Tensor4 itw(
1611  p_gridpos.nelem(), lat_gridpos.nelem(), lon_gridpos.nelem(), 8);
1612  interpweights(itw, p_gridpos, lat_gridpos, lon_gridpos);
1613 
1614  Tensor3 newfield(
1615  af_p_grid.nelem(), af_lat_grid.nelem(), af_lon_grid.nelem());
1616  interp(newfield, itw, species.data, p_gridpos, lat_gridpos, lon_gridpos);
1617 
1618  atm_fields_compact.data(insert_pos, joker, joker, joker) = newfield;
1619  } else { // 2D-case
1620 
1621  Tensor3 itw(p_gridpos.nelem(), lat_gridpos.nelem(), 4);
1622  interpweights(itw, p_gridpos, lat_gridpos);
1623 
1624  Matrix newfield(af_p_grid.nelem(), af_lat_grid.nelem());
1625  interp(
1626  newfield, itw, species.data(joker, joker, 0), p_gridpos, lat_gridpos);
1627 
1628  atm_fields_compact.data(insert_pos, joker, joker, 0) = newfield;
1629  }
1630  } else { // 1D-case
1631  Matrix itw(p_gridpos.nelem(), 2);
1632  interpweights(itw, p_gridpos);
1633 
1634  Vector newfield(af_p_grid.nelem());
1635  interp(newfield, itw, species.data(joker, 0, 0), p_gridpos);
1636 
1637  atm_fields_compact.data(insert_pos, joker, 0, 0) = newfield;
1638  }
1639 }
1640 
1641 // Workspace method, doxygen header is auto-generated
1642 // 2015-06-30 Jana Mendrok
1643 void atm_fields_compactCleanup( //WS Output:
1644  GriddedField4& atm_fields_compact,
1645  //WS Input:
1646  const Numeric& threshold,
1647  const Verbosity&) {
1648  assert(atm_fields_compact.checksize());
1649  Tensor4View afd = atm_fields_compact.data;
1650 
1651  // Check that the data tensor does not contain realistically low (e.g.
1652  // negative) values. Values smaller than threshold will be set to 0).
1653  // Ignore T and z, though.
1654  for (Index i = 0; i < afd.nbooks(); i++)
1655  if (atm_fields_compact.get_string_grid(GFIELD4_FIELD_NAMES)[i] != "T" &&
1656  atm_fields_compact.get_string_grid(GFIELD4_FIELD_NAMES)[i] != "z")
1657  for (Index j = 0; j < afd.npages(); j++)
1658  for (Index k = 0; k < afd.nrows(); k++)
1659  for (Index l = 0; l < afd.ncols(); l++)
1660  if (afd(i, j, k, l) < threshold) afd(i, j, k, l) = 0.0;
1661 }
1662 
1663 // Workspace method, doxygen header is auto-generated
1665  GriddedField4& atm_fields_compact,
1666  // WS Generic Input:
1667  const String& name,
1668  const GriddedField3& field,
1669  const Verbosity&) {
1670  assert(field.checksize());
1671 
1672  ConstVectorView sp_p_grid = field.get_numeric_grid(GFIELD3_P_GRID);
1673  ConstVectorView sp_lat_grid = field.get_numeric_grid(GFIELD3_LAT_GRID);
1674  ConstVectorView sp_lon_grid = field.get_numeric_grid(GFIELD3_LON_GRID);
1675  ArrayOfString sp_name_grid(1);
1676  sp_name_grid[0] = name;
1677 
1678  atm_fields_compact.set_grid(0, sp_name_grid);
1679  atm_fields_compact.set_grid(1, sp_p_grid);
1680  atm_fields_compact.set_grid(2, sp_lat_grid);
1681  atm_fields_compact.set_grid(3, sp_lon_grid);
1682 
1683  atm_fields_compact.data.resize(
1684  1, sp_p_grid.nelem(), sp_lat_grid.nelem(), sp_lon_grid.nelem());
1685 
1686  atm_fields_compact.data(0, joker, joker, joker) = field.data;
1687 }
1688 
1689 // Workspace method, doxygen header is auto-generated
1690 // 2011-05-11 Gerrit Holl
1692  ArrayOfGriddedField4& batch_atm_fields_compact,
1693  // WS Generic Input:
1694  const String& name,
1695  const Numeric& value,
1696  const Index& prepend,
1697  const ArrayOfString& condensibles,
1698  const Verbosity& verbosity) {
1699  for (Index i = 0; i < batch_atm_fields_compact.nelem(); i++) {
1700  atm_fields_compactAddConstant(batch_atm_fields_compact[i],
1701  name,
1702  value,
1703  prepend,
1704  condensibles,
1705  verbosity);
1706  }
1707 }
1708 
1709 // Workspace method, doxygen header is auto-generated
1710 // 2011-05-09 Gerrit Holl
1712  ArrayOfGriddedField4& batch_atm_fields_compact,
1713  // WS Generic Input:
1714  const String& name,
1715  const GriddedField3& species,
1716  const Index& prepend,
1717  const Verbosity& verbosity) {
1718  const Index nelem = batch_atm_fields_compact.nelem();
1719 
1720  String fail_msg;
1721  bool failed = false;
1722 
1723  // Parallelise this for-loop (some interpolation is being done, so it may
1724  // be beneficial)
1725 #pragma omp parallel for if (!arts_omp_in_parallel() && \
1726  nelem >= arts_omp_get_max_threads())
1727  for (Index i = 0; i < nelem; i++) {
1728  try {
1730  batch_atm_fields_compact[i], name, species, prepend, verbosity);
1731  } catch (const std::exception& e) {
1732 #pragma omp critical(batch_atm_fields_compactAddSpecies_fail)
1733  {
1734  fail_msg = e.what();
1735  failed = true;
1736  }
1737  }
1738  }
1739 
1740  if (failed) throw runtime_error(fail_msg);
1741 }
1742 
1743 // Workspace method, doxygen header is auto-generated
1744 // 2015-06-30 Jana Mendrok
1746  ArrayOfGriddedField4& batch_atm_fields_compact,
1747  //WS Input:
1748  const Numeric& threshold,
1749  const Verbosity& verbosity) {
1750  for (Index i = 0; i < batch_atm_fields_compact.nelem(); i++) {
1752  batch_atm_fields_compact[i], threshold, verbosity);
1753  }
1754 }
1755 
1756 // Workspace method, doxygen header is auto-generated.
1758  ArrayOfGriddedField4& batch_atm_fields_compact,
1759  // WS Input:
1760  const Index& atmosphere_dim,
1761  // WS Generic Input:
1762  const ArrayOfMatrix& am,
1763  // Control Parameters:
1764  const ArrayOfString& field_names,
1765  const Verbosity& verbosity) {
1766  const Index amnelem = am.nelem();
1767 
1768  if (amnelem == 0) {
1769  ostringstream os;
1770  os << "No elements in atmospheric scenario batch.\n"
1771  << "Check, whether any batch atmosphere file has been read!";
1772  throw runtime_error(os.str());
1773  }
1774 
1775  // We use the existing WSMs atm_fields_compactFromMatrix and
1776  // atm_fields_compactAddConstant to do most of the work.
1777 
1778  // Make output variable the proper size:
1779  batch_atm_fields_compact.resize(amnelem);
1780 
1781  String fail_msg;
1782  bool failed = false;
1783 
1784  // Loop the batch cases:
1785 #pragma omp parallel for if (!arts_omp_in_parallel() && \
1786  amnelem >= arts_omp_get_max_threads())
1787  for (Index i = 0; i < amnelem; ++i) {
1788  // Skip remaining iterations if an error occurred
1789  if (failed) continue;
1790 
1791  // All the input variables are visible here, despite the
1792  // "default(none)". The reason is that they are return by
1793  // reference arguments of this function, which are shared
1794  // automatically.
1795 
1796  // The try block here is necessary to correctly handle
1797  // exceptions inside the parallel region.
1798  try {
1799  atm_fields_compactFromMatrix(batch_atm_fields_compact[i],
1800  atmosphere_dim,
1801  am[i],
1802  field_names,
1803  verbosity);
1804  } catch (const std::exception& e) {
1805 #pragma omp critical(batch_atm_fields_compactFromArrayOfMatrix_fail)
1806  {
1807  fail_msg = e.what();
1808  failed = true;
1809  }
1810  }
1811  }
1812 
1813  if (failed) throw runtime_error(fail_msg);
1814 }
1815 
1816 /* Workspace method: Doxygen documentation will be auto-generated */
1818  Vector& p_grid,
1819  Vector& lat_grid,
1820  Vector& lon_grid,
1821  Tensor3& t_field,
1822  Tensor3& z_field,
1823  Tensor4& vmr_field,
1824  Tensor4& particle_bulkprop_field,
1825  ArrayOfString& particle_bulkprop_names,
1826 
1827  // WS Input:
1828  const ArrayOfArrayOfSpeciesTag& abs_species,
1829  const GriddedField4& atm_fields_compact,
1830  const Index& atmosphere_dim,
1831  const String& delim,
1832  const Numeric& p_min,
1833  // Control parameters:
1834  const Index& check_gridnames,
1835  const Verbosity&) {
1836  // Make a handle on atm_fields_compact to save typing:
1837  const GriddedField4& c = atm_fields_compact;
1838 
1839  // Check if the grids in our data match atmosphere_dim
1840  // (throws an error if the dimensionality is not correct):
1841  chk_atm_grids(atmosphere_dim,
1845 
1846  // Optional check for gridnames.
1847  if (check_gridnames == 1) {
1848  chk_griddedfield_gridname(c, 1, "Pressure");
1849  chk_griddedfield_gridname(c, 2, "Latitude");
1850  chk_griddedfield_gridname(c, 3, "Longitude");
1851  }
1852 
1853  const Index nf = c.get_grid_size(GFIELD4_FIELD_NAMES);
1854  const Index np = c.get_grid_size(GFIELD4_P_GRID);
1855  // const Index nlat = c.get_grid_size(GFIELD4_LAT_GRID);
1856  // const Index nlon = c.get_grid_size(GFIELD4_LON_GRID);
1859  if (nlat == 0) nlat++;
1860  if (nlon == 0) nlon++;
1861 
1862  // Grids:
1863  p_grid = c.get_numeric_grid(GFIELD4_P_GRID);
1864  lat_grid = c.get_numeric_grid(GFIELD4_LAT_GRID);
1865  lon_grid = c.get_numeric_grid(GFIELD4_LON_GRID);
1866 
1867  // Reduce p_grid to region below p_min
1868  Index l = np - 1;
1869  bool search_toa = true;
1870  while (search_toa && l > 0) {
1871  if (p_grid[l - 1] < p_min)
1872  l--;
1873  else
1874  search_toa = false;
1875  }
1876  if (search_toa) {
1877  ostringstream os;
1878  os << "At least one atmospheric level with pressure larger p_min (="
1879  << p_min << ")\n"
1880  << "is needed, but none is found.";
1881  throw runtime_error(os.str());
1882  }
1883  const Index npn = l + 1;
1884  p_grid = p_grid[Range(0, npn)];
1885 
1886  const Index nsa = abs_species.nelem();
1887 
1888  // Check that there is at least one VMR species:
1889  if (nsa < 1) {
1890  ostringstream os;
1891  os << "There must be at least one absorption species.";
1892  throw runtime_error(os.str());
1893  }
1894 
1895  // In contrast to older versions, we allow the data entries to be in arbitrary
1896  // order. that is, we have to look for all fields individually. we require the
1897  // abs_species related fields as well as the scat_species related fields to
1898  // have leading identifiers, namely 'abs_species' and 'scat_species'. it is
1899  // not mandatory that all fields available in atm_field_compact are applied
1900  // through abs_species and scat_species; left over fields are silently
1901  // ignored.
1902  // For temperature and altitude, occurence of exactly one field entry is
1903  // ensured. For other fields, the first match is used.
1904  // FIXME: or should a match be removed such that a second species occurence
1905  // would match a later-in-line field?
1906 
1907  bool found;
1908  const String as_type = "abs_species";
1909  const String ss_type = "scat_species";
1910 
1911  // Find temperature field:
1912  found = false;
1913  t_field.resize(npn, nlat, nlon);
1914  for (Index i = 0; i < nf; ++i) {
1915  if (c.get_string_grid(GFIELD4_FIELD_NAMES)[i] == "T") {
1916  if (found) {
1917  ostringstream os;
1918  os << "Only one temperature ('T') field allowed,\n"
1919  << "but found at least 2.";
1920  throw runtime_error(os.str());
1921  } else {
1922  found = true;
1923  t_field = c.data(i, Range(0, npn), Range(joker), Range(joker));
1924  }
1925  }
1926  }
1927  if (!found) {
1928  ostringstream os;
1929  os << "One temperature ('T') field required, but none found";
1930  throw runtime_error(os.str());
1931  }
1932 
1933  // Find Altitude field:
1934  found = false;
1935  z_field.resize(npn, nlat, nlon);
1936  for (Index i = 0; i < nf; ++i) {
1937  if (c.get_string_grid(GFIELD4_FIELD_NAMES)[i] == "z") {
1938  if (found) {
1939  ostringstream os;
1940  os << "Only one altitude ('z') field allowed,\n"
1941  << "but found at least 2.";
1942  throw runtime_error(os.str());
1943  } else {
1944  found = true;
1945  z_field = c.data(i, Range(0, npn), Range(joker), Range(joker));
1946  }
1947  }
1948  }
1949  if (!found) {
1950  ostringstream os;
1951  os << "One altitude ('z') field required, but none found";
1952  throw runtime_error(os.str());
1953  }
1954 
1955  // Extracting the required abs_species fields:
1956 
1957  vmr_field.resize(nsa, npn, nlat, nlon);
1958  for (Index j = 0; j < nsa; ++j) {
1959  using global_data::species_data; // The species lookup data:
1960  const String as_name = species_data[abs_species[j][0].Species()].Name();
1961  found = false;
1962  Index i = 0;
1963  String species_type;
1964  String species_name;
1965  while (!found && i < nf) {
1967  species_type, c.get_string_grid(GFIELD4_FIELD_NAMES)[i], delim);
1968  // do we have an abs_species type field?
1969  if (species_type == as_type) {
1971  species_name, c.get_string_grid(GFIELD4_FIELD_NAMES)[i], delim);
1972  if (species_name == as_name) {
1973  found = true;
1974  vmr_field(j, Range(joker), Range(joker), Range(joker)) =
1975  c.data(i, Range(0, npn), Range(joker), Range(joker));
1976  }
1977  }
1978  i++;
1979  }
1980  if (!found) {
1981  ostringstream os;
1982  os << "No field for absorption species '" << as_name << "' found.";
1983  throw runtime_error(os.str());
1984  }
1985  }
1986 
1987  //get the number of scat_species entries
1988  std::vector<Index> Idx;
1989 
1990  String species_type;
1991  for (Index i = 0; i < nf; ++i) {
1993  species_type, c.get_string_grid(GFIELD4_FIELD_NAMES)[i], delim);
1994 
1995  if (species_type == ss_type) {
1996  Idx.push_back(i);
1997  }
1998  }
1999 
2000  const Index nsp = Idx.size();
2001 
2002  // Extracting the required scattering species fields:
2003  particle_bulkprop_field.resize(nsp, npn, nlat, nlon);
2004  particle_bulkprop_field = NAN;
2005  particle_bulkprop_names.resize(nsp);
2006 
2007  // put scat_species entries in particle_bulkprop_field
2008 
2009  for (Index j = 0; j < nsp; ++j) {
2010  String species_name;
2011  String scat_type;
2012 
2014  scat_type, c.get_string_grid(GFIELD4_FIELD_NAMES)[Idx[j]], delim);
2015 
2017  species_name, c.get_string_grid(GFIELD4_FIELD_NAMES)[Idx[j]], delim);
2018 
2019  particle_bulkprop_field(j, Range(joker), Range(joker), Range(joker)) =
2020  c.data(Idx[j], Range(0, npn), Range(joker), Range(joker));
2021 
2022  particle_bulkprop_names[j] = species_name + delim + scat_type;
2023  }
2024 }
2025 
2026 /* Workspace method: Doxygen documentation will be auto-generated */
2027 void AtmosphereSet1D( // WS Output:
2028  Index& atmosphere_dim,
2029  Vector& lat_grid,
2030  Vector& lon_grid,
2031  const Verbosity& verbosity) {
2032  CREATE_OUT2;
2033  CREATE_OUT3;
2034 
2035  out2 << " Sets the atmospheric dimensionality to 1.\n";
2036  out3 << " atmosphere_dim = 1\n";
2037  out3 << " lat_grid is set to be an empty vector\n";
2038  out3 << " lon_grid is set to be an empty vector\n";
2039 
2040  atmosphere_dim = 1;
2041  lat_grid.resize(0);
2042  lon_grid.resize(0);
2043 }
2044 
2045 /* Workspace method: Doxygen documentation will be auto-generated */
2046 void AtmosphereSet2D( // WS Output:
2047  Index& atmosphere_dim,
2048  Vector& lon_grid,
2049  const Verbosity& verbosity) {
2050  CREATE_OUT2;
2051  CREATE_OUT3;
2052 
2053  out2 << " Sets the atmospheric dimensionality to 2.\n";
2054  out3 << " atmosphere_dim = 2\n";
2055  out3 << " lon_grid is set to be an empty vector\n";
2056 
2057  atmosphere_dim = 2;
2058  lon_grid.resize(0);
2059 }
2060 
2061 /* Workspace method: Doxygen documentation will be auto-generated */
2062 void AtmosphereSet3D( // WS Output:
2063  Index& atmosphere_dim,
2064  Vector& lat_true,
2065  Vector& lon_true,
2066  const Verbosity& verbosity) {
2067  CREATE_OUT2;
2068  CREATE_OUT3;
2069 
2070  out2 << " Sets the atmospheric dimensionality to 3.\n";
2071  out3 << " atmosphere_dim = 3\n";
2072 
2073  atmosphere_dim = 3;
2074  lat_true.resize(0);
2075  lon_true.resize(0);
2076 }
2077 
2078 /* Workspace method: Doxygen documentation will be auto-generated */
2079 void AtmFieldsCalc( //WS Output:
2080  Tensor3& t_field,
2081  Tensor3& z_field,
2082  Tensor4& vmr_field,
2083  EnergyLevelMap& nlte_field,
2084  //WS Input
2085  const Vector& p_grid,
2086  const Vector& lat_grid,
2087  const Vector& lon_grid,
2088  const GriddedField3& t_field_raw,
2089  const GriddedField3& z_field_raw,
2090  const ArrayOfGriddedField3& vmr_field_raw,
2091  const ArrayOfGriddedField3& nlte_field_raw,
2092  const ArrayOfQuantumIdentifier& nlte_ids,
2093  const Vector& nlte_energies,
2094  const Index& atmosphere_dim,
2095  // WS Generic Input:
2096  const Index& interp_order,
2097  const Index& vmr_zeropadding,
2098  const Index& vmr_nonegative,
2099  const Index& nlte_when_negative,
2100  const Verbosity& verbosity) {
2101  CREATE_OUT2;
2102 
2103  const Vector& tfr_p_grid =
2104  t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2105  const Vector& tfr_lat_grid =
2106  t_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2107  const Vector& tfr_lon_grid =
2108  t_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2109  const Vector& zfr_p_grid =
2110  z_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2111  const Vector& zfr_lat_grid =
2112  z_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2113  const Vector& zfr_lon_grid =
2114  z_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2115 
2116  out2 << " Interpolation order: " << interp_order << "\n";
2117 
2118  // Basic checks of input variables
2119  //
2120  // Atmosphere
2121  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2122  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
2123 
2124  // NLTE basics
2125  nlte_field.Type() = nlte_ids.nelem() == nlte_field_raw.nelem() ? EnergyLevelMapType::Tensor3_t : EnergyLevelMapType::None_t;
2126  nlte_field.Levels() = nlte_ids.nelem() == nlte_field_raw.nelem() ? nlte_ids : ArrayOfQuantumIdentifier(0);
2127  nlte_field.Energies() = nlte_ids.nelem() == nlte_field_raw.nelem() ? nlte_energies : Vector(0);
2128 
2129  //==========================================================================
2130  if (atmosphere_dim == 1) {
2131  if (!(tfr_lat_grid.nelem() == 1 && tfr_lon_grid.nelem() == 1))
2132  throw runtime_error(
2133  "Temperature data (T_field) has wrong dimension "
2134  "(2D or 3D).\n");
2135 
2136  if (!(zfr_lat_grid.nelem() == 1 && zfr_lon_grid.nelem() == 1))
2137  throw runtime_error(
2138  "Altitude data (z_field) has wrong dimension "
2139  "(2D or 3D).\n");
2140 
2141  GriddedField3 temp_gfield3;
2142 
2144  temp_gfield3, p_grid, t_field_raw, interp_order, 0, verbosity);
2145  t_field = temp_gfield3.data;
2146 
2148  temp_gfield3, p_grid, z_field_raw, interp_order, 0, verbosity);
2149  z_field = temp_gfield3.data;
2150 
2151  ArrayOfGriddedField3 temp_agfield3;
2152  try {
2153  GriddedFieldPRegrid(temp_agfield3,
2154  p_grid,
2155  vmr_field_raw,
2156  interp_order,
2157  vmr_zeropadding,
2158  verbosity);
2159  } catch (const std::runtime_error& e) {
2160  ostringstream os;
2161  os << e.what() << "\n"
2162  << "Note that you can explicitly set vmr_zeropadding "
2163  << "to 1 in the method call.";
2164  throw runtime_error(os.str());
2165  }
2167  vmr_field, p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
2168 
2169  // Non-LTE interpolation
2170  if (nlte_ids.nelem() == nlte_field_raw.nelem()) {
2172  temp_agfield3, p_grid, nlte_field_raw, interp_order, 0, verbosity);
2174  nlte_field.Data(), p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
2175  }
2176  else
2177  nlte_field.Data().resize(0, 0, 0, 0);
2178 
2179  }
2180 
2181  //=========================================================================
2182  else if (atmosphere_dim == 2) {
2183  if (tfr_lat_grid.nelem() == 1 && tfr_lon_grid.nelem() == 1)
2184  throw runtime_error(
2185  "Raw data has wrong dimension (1D). "
2186  "You have to use \n"
2187  "AtmFieldsCalcExpand1D instead of AtmFieldsCalc.");
2188 
2189  //Resize variables
2190  t_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2191  z_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2192  vmr_field.resize(
2193  vmr_field_raw.nelem(), p_grid.nelem(), lat_grid.nelem(), 1);
2194  if (nlte_ids.nelem() == nlte_field_raw.nelem())
2195  nlte_field.Data().resize(
2196  nlte_field_raw.nelem(), p_grid.nelem(), lat_grid.nelem(), 1);
2197  else
2198  nlte_field.Data().resize(0, 0, 0, 0);
2199 
2200  // Gridpositions:
2201  ArrayOfGridPosPoly gp_p(p_grid.nelem());
2202  ArrayOfGridPosPoly gp_lat(lat_grid.nelem());
2203 
2204  // Interpolate t_field:
2205 
2206  // Check that interpolation grids are ok (and throw a detailed
2207  // error message if not):
2209  "Raw temperature to p_grid, 2D case", tfr_p_grid, p_grid, interp_order);
2210  chk_interpolation_grids("Raw temperature to lat_grid, 2D case",
2211  tfr_lat_grid,
2212  lat_grid,
2213  interp_order);
2214 
2215  // Calculate grid positions:
2216  p2gridpos_poly(gp_p, tfr_p_grid, p_grid, interp_order);
2217  gridpos_poly(gp_lat, tfr_lat_grid, lat_grid, interp_order);
2218 
2219  // Interpolation weights:
2220  Tensor3 itw(p_grid.nelem(),
2221  lat_grid.nelem(),
2222  (interp_order + 1) * (interp_order + 1));
2223  // (4 interpolation weights are required for example for linear 2D
2224  // interpolation)
2225  interpweights(itw, gp_p, gp_lat);
2226 
2227  // Interpolate:
2228  interp(t_field(joker, joker, 0),
2229  itw,
2230  t_field_raw.data(joker, joker, 0),
2231  gp_p,
2232  gp_lat);
2233 
2234  // Interpolate z_field:
2235 
2236  // Check that interpolation grids are ok (and throw a detailed
2237  // error message if not):
2239  "Raw z to p_grid, 2D case", zfr_p_grid, p_grid, interp_order);
2241  "Raw z to lat_grid, 2D case", zfr_lat_grid, lat_grid, interp_order);
2242 
2243  // Calculate grid positions:
2244  p2gridpos_poly(gp_p, zfr_p_grid, p_grid, interp_order);
2245  gridpos_poly(gp_lat, zfr_lat_grid, lat_grid, interp_order);
2246 
2247  // Interpolation weights:
2248  interpweights(itw, gp_p, gp_lat);
2249 
2250  // Interpolate:
2251  interp(z_field(joker, joker, 0),
2252  itw,
2253  z_field_raw.data(joker, joker, 0),
2254  gp_p,
2255  gp_lat);
2256 
2257  // Interpolate vmr_field.
2258  // Loop over the gaseous species:
2259  for (Index gas_i = 0; gas_i < vmr_field_raw.nelem(); gas_i++) {
2260  ostringstream os;
2261 
2262  if (!(vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID).nelem() !=
2263  1 &&
2264  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LON_GRID).nelem() ==
2265  1)) {
2266  os << "VMR data of the " << gas_i << " the species has "
2267  << "wrong dimension (1D or 3D). \n";
2268  throw runtime_error(os.str());
2269  }
2270 
2271  // Check that interpolation grids are ok (and throw a detailed
2272  // error message if not):
2273  os << "Raw VMR[" << gas_i << "] to p_grid, 2D case";
2275  os.str(),
2276  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_P_GRID),
2277  p_grid,
2278  interp_order);
2279  os.str("");
2280  os << "Raw VMR[" << gas_i << "] to lat_grid, 2D case";
2282  os.str(),
2283  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID),
2284  lat_grid,
2285  interp_order);
2286 
2287  // Calculate grid positions:
2288  p2gridpos_poly(gp_p,
2289  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_P_GRID),
2290  p_grid,
2291  interp_order);
2292  gridpos_poly(gp_lat,
2293  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID),
2294  lat_grid,
2295  interp_order);
2296 
2297  // Interpolation weights:
2298  interpweights(itw, gp_p, gp_lat);
2299 
2300  // Interpolate:
2301  interp(vmr_field(gas_i, joker, joker, 0),
2302  itw,
2303  vmr_field_raw[gas_i].data(joker, joker, 0),
2304  gp_p,
2305  gp_lat);
2306  }
2307 
2308  // Interpolate Non-LTE
2309  for (Index qi_i = 0; qi_i < nlte_field_raw.nelem(); qi_i++) {
2310  ostringstream os;
2311 
2312  if (!(nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_LAT_GRID).nelem() !=
2313  1 &&
2314  nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_LON_GRID).nelem() ==
2315  1)) {
2316  os << "NLTE data of the " << qi_i << " temperature field has "
2317  << "wrong dimension (1D or 3D). \n";
2318  throw std::runtime_error(os.str());
2319  }
2320 
2321  // Check that interpolation grids are ok (and throw a detailed
2322  // error message if not):
2323  os << "Raw NLTE[" << qi_i << "] to p_grid, 2D case";
2325  os.str(),
2326  nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_P_GRID),
2327  p_grid,
2328  interp_order);
2329  os.str("");
2330  os << "Raw NLTE[" << qi_i << "] to lat_grid, 2D case";
2332  os.str(),
2333  nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_LAT_GRID),
2334  lat_grid,
2335  interp_order);
2336 
2337  // Calculate grid positions:
2338  p2gridpos_poly(gp_p,
2339  nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_P_GRID),
2340  p_grid,
2341  interp_order);
2342  gridpos_poly(gp_lat,
2343  nlte_field_raw[qi_i].get_numeric_grid(GFIELD3_LAT_GRID),
2344  lat_grid,
2345  interp_order);
2346 
2347  // Interpolation weights:
2348  interpweights(itw, gp_p, gp_lat);
2349 
2350  // Interpolate:
2351  if (nlte_ids.nelem() == nlte_field_raw.nelem())
2352  interp(nlte_field.Data()(qi_i, joker, joker, 0),
2353  itw,
2354  nlte_field_raw[qi_i].data(joker, joker, 0),
2355  gp_p,
2356  gp_lat);
2357  else
2358  nlte_field.Data().resize(0, 0, 0, 0);
2359  }
2360  }
2361 
2362  //================================================================
2363  // atmosphere_dim = 3
2364  else if (atmosphere_dim == 3) {
2365  if (tfr_lat_grid.nelem() == 1 && tfr_lon_grid.nelem() == 1)
2366  throw runtime_error(
2367  "Raw data has wrong dimension. You have to use \n"
2368  "AtmFieldsCalcExpand1D instead of AtmFieldsCalc.");
2369 
2370  GriddedField3 temp_gfield3;
2371 
2373  temp_gfield3, lat_grid, lon_grid, t_field_raw, interp_order, verbosity);
2375  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
2376  t_field = temp_gfield3.data;
2377 
2379  temp_gfield3, lat_grid, lon_grid, z_field_raw, interp_order, verbosity);
2381  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
2382  z_field = temp_gfield3.data;
2383 
2384  ArrayOfGriddedField3 temp_agfield3;
2385  GriddedFieldLatLonRegrid(temp_agfield3,
2386  lat_grid,
2387  lon_grid,
2388  vmr_field_raw,
2389  interp_order,
2390  verbosity);
2391  try {
2392  GriddedFieldPRegrid(temp_agfield3,
2393  p_grid,
2394  temp_agfield3,
2395  interp_order,
2396  vmr_zeropadding,
2397  verbosity);
2398  } catch (const std::runtime_error& e) {
2399  ostringstream os;
2400  os << e.what() << "\n"
2401  << "Note that you can explicitly set vmr_zeropadding "
2402  << "to 1 in the method call.";
2403  throw runtime_error(os.str());
2404  }
2406  vmr_field, p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
2407 
2408  if (nlte_field_raw.nelem()) {
2409  GriddedFieldLatLonRegrid(temp_agfield3,
2410  lat_grid,
2411  lon_grid,
2412  nlte_field_raw,
2413  interp_order,
2414  verbosity);
2415 
2416  try {
2418  temp_agfield3, p_grid, temp_agfield3, interp_order, 0, verbosity);
2419  } catch (const std::runtime_error& e) {
2420  ostringstream os;
2421  os << e.what() << "\n"
2422  << "Note that you can explicitly set vmr_zeropadding "
2423  << "to 1 in the method call.";
2424  throw runtime_error(os.str());
2425  }
2426 
2427  if (nlte_ids.nelem() == nlte_field_raw.nelem())
2429  nlte_field.Data(), p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
2430  else
2431  nlte_field.Data().resize(0, 0, 0, 0);
2432  }
2433  } else {
2434  // We can never get here, since there was a runtime
2435  // error check for atmosphere_dim at the beginning.
2436  assert(false);
2437  }
2438 
2439  // remove negatives?
2440  if (vmr_nonegative) {
2441  for (Index ib = 0; ib < vmr_field.nbooks(); ib++) {
2442  for (Index ip = 0; ip < vmr_field.npages(); ip++) {
2443  for (Index ir = 0; ir < vmr_field.nrows(); ir++) {
2444  for (Index ic = 0; ic < vmr_field.ncols(); ic++) {
2445  if (vmr_field(ib, ip, ir, ic) < 0) {
2446  vmr_field(ib, ip, ir, ic) = 0;
2447  }
2448  }
2449  }
2450  }
2451  }
2452  }
2453 
2454  // what to do with negative nlte temperatures?
2455  if (nlte_when_negative != -1) {
2456  if (nlte_field_raw.nelem()) {
2457  for (Index ib = 0; ib < nlte_field.Data().nbooks(); ib++) {
2458  for (Index ip = 0; ip < nlte_field.Data().npages(); ip++) {
2459  for (Index ir = 0; ir < nlte_field.Data().nrows(); ir++) {
2460  for (Index ic = 0; ic < nlte_field.Data().ncols(); ic++) {
2461  if (nlte_field.Data()(ib, ip, ir, ic) < 0) {
2462  // Set to atmospheric temperature or to nil.
2463  nlte_field.Data()(ib, ip, ir, ic) =
2464  nlte_when_negative == 1 ? t_field(ip, ir, ic) : 0;
2465  // NOTE: This only makes sense for vibrational NLTE and is bad elsewise
2466  // but since elsewise is bad anyways with negative values, it is
2467  // and acceptable compromise...
2468  }
2469  }
2470  }
2471  }
2472  }
2473  }
2474  }
2475 
2476  nlte_field.ThrowIfNotOK();
2477 }
2478 
2479 /* Workspace method: Doxygen documentation will be auto-generated */
2480 void MagFieldsCalc( //WS Output:
2481  Tensor3& mag_u_field,
2482  Tensor3& mag_v_field,
2483  Tensor3& mag_w_field,
2484  //WS Input
2485  const Vector& p_grid,
2486  const Vector& lat_grid,
2487  const Vector& lon_grid,
2488  const GriddedField3& mag_u_field_raw,
2489  const GriddedField3& mag_v_field_raw,
2490  const GriddedField3& mag_w_field_raw,
2491  const Index& atmosphere_dim,
2492  // WS Generic Input:
2493  const Index& interp_order,
2494  const Verbosity& verbosity) {
2495  CREATE_OUT2;
2496 
2497  const Vector& ufr_p_grid =
2498  mag_u_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2499  const Vector& ufr_lat_grid =
2500  mag_u_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2501  const Vector& ufr_lon_grid =
2502  mag_u_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2503  const Vector& vfr_p_grid =
2504  mag_v_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2505  const Vector& vfr_lat_grid =
2506  mag_v_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2507  const Vector& vfr_lon_grid =
2508  mag_v_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2509  const Vector& wfr_p_grid =
2510  mag_w_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2511  const Vector& wfr_lat_grid =
2512  mag_w_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2513  const Vector& wfr_lon_grid =
2514  mag_w_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2515 
2516  out2 << " Interpolation order: " << interp_order << "\n";
2517 
2518  // Basic checks of input variables
2519  //
2520  // Atmosphere
2521  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2522  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
2523 
2524  //==========================================================================
2525  if (atmosphere_dim == 1) {
2526  if (!(ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1))
2527  throw std::runtime_error(
2528  "Magnetic u field data has wrong dimension (2D or 3D).\n");
2529  if (!(vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1))
2530  throw std::runtime_error(
2531  "Magnetic v field data has wrong dimension (2D or 3D).\n");
2532  if (!(wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1))
2533  throw std::runtime_error(
2534  "Magnetic w field data has wrong dimension (2D or 3D).\n");
2535 
2536  GriddedField3 temp_gfield3;
2537 
2539  temp_gfield3, p_grid, mag_u_field_raw, interp_order, 0, verbosity);
2540  mag_u_field = temp_gfield3.data;
2541 
2543  temp_gfield3, p_grid, mag_v_field_raw, interp_order, 0, verbosity);
2544  mag_v_field = temp_gfield3.data;
2545 
2547  temp_gfield3, p_grid, mag_w_field_raw, interp_order, 0, verbosity);
2548  mag_w_field = temp_gfield3.data;
2549  }
2550 
2551  //=========================================================================
2552  else if (atmosphere_dim == 2) {
2553  if (ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1)
2554  throw std::runtime_error(
2555  "Raw data has wrong dimension (1D). You have to use \n"
2556  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2557  if (vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1)
2558  throw std::runtime_error(
2559  "Raw data has wrong dimension (1D). You have to use \n"
2560  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2561  if (wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1)
2562  throw std::runtime_error(
2563  "Raw data has wrong dimension (1D). You have to use \n"
2564  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2565 
2566  //Resize variables
2567  mag_u_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2568  mag_v_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2569  mag_w_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2570 
2571  // Gridpositions:
2572  ArrayOfGridPosPoly gp_p(p_grid.nelem());
2573  ArrayOfGridPosPoly gp_lat(lat_grid.nelem());
2574 
2575  // Interpolate mag_u_field:
2576 
2577  // Check that interpolation grids are ok (and throw a detailed
2578  // error message if not):
2580  "Raw u field to p_grid, 2D case", ufr_p_grid, p_grid, interp_order);
2581  chk_interpolation_grids("Raw u field to lat_grid, 2D case",
2582  ufr_lat_grid,
2583  lat_grid,
2584  interp_order);
2585 
2586  // Calculate grid positions:
2587  p2gridpos_poly(gp_p, ufr_p_grid, p_grid, interp_order);
2588  gridpos_poly(gp_lat, ufr_lat_grid, lat_grid, interp_order);
2589 
2590  // Interpolation weights:
2591  Tensor3 itwu(p_grid.nelem(),
2592  lat_grid.nelem(),
2593  (interp_order + 1) * (interp_order + 1));
2594  // (4 interpolation weights are required for example for linear 2D
2595  // interpolation)
2596  interpweights(itwu, gp_p, gp_lat);
2597 
2598  // Interpolate:
2599  interp(mag_u_field(joker, joker, 0),
2600  itwu,
2601  mag_u_field_raw.data(joker, joker, 0),
2602  gp_p,
2603  gp_lat);
2604 
2605  // Interpolate mag_v_field:
2606 
2607  // Check that interpolation grids are ok (and throw a detailed
2608  // error message if not):
2610  "Raw v field to p_grid, 2D case", vfr_p_grid, p_grid, interp_order);
2611  chk_interpolation_grids("Raw v field to lat_grid, 2D case",
2612  vfr_lat_grid,
2613  lat_grid,
2614  interp_order);
2615 
2616  // Calculate grid positions:
2617  p2gridpos_poly(gp_p, vfr_p_grid, p_grid, interp_order);
2618  gridpos_poly(gp_lat, vfr_lat_grid, lat_grid, interp_order);
2619 
2620  // Interpolation weights:
2621  Tensor3 itwv(p_grid.nelem(),
2622  lat_grid.nelem(),
2623  (interp_order + 1) * (interp_order + 1));
2624  // (4 interpolation weights are required for example for linear 2D
2625  // interpolation)
2626  interpweights(itwv, gp_p, gp_lat);
2627 
2628  // Interpolate:
2629  interp(mag_v_field(joker, joker, 0),
2630  itwv,
2631  mag_v_field_raw.data(joker, joker, 0),
2632  gp_p,
2633  gp_lat);
2634 
2635  // Interpolate mag_w_field:
2636 
2637  // Check that interpolation grids are ok (and throw a detailed
2638  // error message if not):
2640  "Raw w field to p_grid, 2D case", wfr_p_grid, p_grid, interp_order);
2641  chk_interpolation_grids("Raw w field to lat_grid, 2D case",
2642  wfr_lat_grid,
2643  lat_grid,
2644  interp_order);
2645 
2646  // Calculate grid positions:
2647  p2gridpos_poly(gp_p, wfr_p_grid, p_grid, interp_order);
2648  gridpos_poly(gp_lat, wfr_lat_grid, lat_grid, interp_order);
2649 
2650  // Interpolation weights:
2651  Tensor3 itww(p_grid.nelem(),
2652  lat_grid.nelem(),
2653  (interp_order + 1) * (interp_order + 1));
2654  // (4 interpolation weights are required for example for linear 2D
2655  // interpolation)
2656  interpweights(itww, gp_p, gp_lat);
2657 
2658  // Interpolate:
2659  interp(mag_w_field(joker, joker, 0),
2660  itww,
2661  mag_w_field_raw.data(joker, joker, 0),
2662  gp_p,
2663  gp_lat);
2664  }
2665 
2666  //================================================================
2667  // atmosphere_dim = 3
2668  else if (atmosphere_dim == 3) {
2669  if (ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1)
2670  throw std::runtime_error(
2671  "Raw data has wrong dimension. You have to use \n"
2672  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2673  if (vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1)
2674  throw std::runtime_error(
2675  "Raw data has wrong dimension. You have to use \n"
2676  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2677  if (wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1)
2678  throw std::runtime_error(
2679  "Raw data has wrong dimension. You have to use \n"
2680  "MagFieldsCalcExpand1D instead of MagFieldsCalc.");
2681 
2682  GriddedField3 temp_gfield3;
2683 
2684  GriddedFieldLatLonRegrid(temp_gfield3,
2685  lat_grid,
2686  lon_grid,
2687  mag_u_field_raw,
2688  interp_order,
2689  verbosity);
2691  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
2692  mag_u_field = temp_gfield3.data;
2693 
2694  GriddedFieldLatLonRegrid(temp_gfield3,
2695  lat_grid,
2696  lon_grid,
2697  mag_v_field_raw,
2698  interp_order,
2699  verbosity);
2701  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
2702  mag_v_field = temp_gfield3.data;
2703 
2704  GriddedFieldLatLonRegrid(temp_gfield3,
2705  lat_grid,
2706  lon_grid,
2707  mag_w_field_raw,
2708  interp_order,
2709  verbosity);
2711  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
2712  mag_w_field = temp_gfield3.data;
2713 
2714  } else {
2715  // We can never get here, since there was a runtime
2716  // error check for atmosphere_dim at the beginning.
2717  assert(false);
2718  }
2719 }
2720 
2721 /* Workspace method: Doxygen documentation will be auto-generated */
2723  Tensor3& mag_u_field,
2724  Tensor3& mag_v_field,
2725  Tensor3& mag_w_field,
2726  //WS Input
2727  const Vector& lat_grid,
2728  const Vector& lon_grid,
2729  const Tensor3& z_field,
2730  const GriddedField3& mag_u_field_raw,
2731  const GriddedField3& mag_v_field_raw,
2732  const GriddedField3& mag_w_field_raw,
2733  // WS Generic Input:
2734  const Index& interp_order,
2735  const Numeric& extrapolation_factor,
2736  const Verbosity& verbosity) {
2737  const auto nalt = z_field.npages();
2738  const auto nlat = z_field.nrows();
2739  const auto nlon = z_field.ncols();
2740 
2741  // Check that the fields are correct
2742  for (auto& gf3 : {mag_u_field_raw, mag_v_field_raw, mag_w_field_raw}) {
2743  if (gf3.get_grid_name(0) not_eq "Altitude" or
2744  gf3.get_grid_name(1) not_eq "Latitude" or
2745  gf3.get_grid_name(2) not_eq "Longitude") {
2746  std::ostringstream os;
2747  os << "Grids are bad\n";
2748  os << "Grids must be Altitude, Latitude, Longitude, but are: "
2749  << gf3.get_grid_name(0) << ", " << gf3.get_grid_name(1) << ", "
2750  << gf3.get_grid_name(2) << '\n';
2751  throw std::runtime_error(os.str());
2752  }
2753  }
2754 
2755  // Regrid and rename raw-fields
2756  GriddedField3 u, v, w;
2758  u, lat_grid, lon_grid, mag_u_field_raw, interp_order, verbosity);
2760  v, lat_grid, lon_grid, mag_v_field_raw, interp_order, verbosity);
2762  w, lat_grid, lon_grid, mag_w_field_raw, interp_order, verbosity);
2763 
2764  // Finally interpolate the three fields
2765  mag_u_field.resize(nalt, nlat, nlon);
2766  mag_v_field.resize(nalt, nlat, nlon);
2767  mag_w_field.resize(nalt, nlat, nlon);
2768  for (Index ilat = 0; ilat < nlat; ilat++) {
2769  for (Index ilon = 0; ilon < nlon; ilon++) {
2770  ArrayOfGridPosPoly gp;
2771  Matrix itw;
2772 
2773  chk_interpolation_grids("Magnetic U Field Altitude",
2774  u.get_numeric_grid(0),
2775  z_field(joker, ilat, ilon),
2776  interp_order,
2777  extrapolation_factor,
2778  false);
2779  gp = ArrayOfGridPosPoly(nalt);
2780  gridpos_poly(gp,
2781  u.get_numeric_grid(0),
2782  z_field(joker, ilat, ilon),
2783  interp_order,
2784  extrapolation_factor);
2785  itw = Matrix(nalt, gp[0].w.nelem());
2786  interpweights(itw, gp);
2787  interp(
2788  mag_u_field(joker, ilat, ilon), itw, u.data(joker, ilat, ilon), gp);
2789 
2790  chk_interpolation_grids("Magnetic V Field Altitude",
2791  v.get_numeric_grid(0),
2792  z_field(joker, ilat, ilon),
2793  interp_order,
2794  extrapolation_factor,
2795  false);
2796  gp = ArrayOfGridPosPoly(nalt);
2797  gridpos_poly(gp,
2798  v.get_numeric_grid(0),
2799  z_field(joker, ilat, ilon),
2800  interp_order,
2801  extrapolation_factor);
2802  itw = Matrix(nalt, gp[0].w.nelem());
2803  interpweights(itw, gp);
2804  interp(
2805  mag_v_field(joker, ilat, ilon), itw, v.data(joker, ilat, ilon), gp);
2806 
2807  chk_interpolation_grids("Magnetic W Field Altitude",
2808  w.get_numeric_grid(0),
2809  z_field(joker, ilat, ilon),
2810  interp_order,
2811  extrapolation_factor,
2812  false);
2813  gp = ArrayOfGridPosPoly(nalt);
2814  gridpos_poly(gp,
2815  w.get_numeric_grid(0),
2816  z_field(joker, ilat, ilon),
2817  interp_order,
2818  extrapolation_factor);
2819  itw = Matrix(nalt, gp[0].w.nelem());
2820  interpweights(itw, gp);
2821  interp(
2822  mag_w_field(joker, ilat, ilon), itw, w.data(joker, ilat, ilon), gp);
2823  }
2824  }
2825 }
2826 
2827 /* Workspace method: Doxygen documentation will be auto-generated */
2828 void WindFieldsCalc( //WS Output:
2829  Tensor3& wind_u_field,
2830  Tensor3& wind_v_field,
2831  Tensor3& wind_w_field,
2832  //WS Input
2833  const Vector& p_grid,
2834  const Vector& lat_grid,
2835  const Vector& lon_grid,
2836  const GriddedField3& wind_u_field_raw,
2837  const GriddedField3& wind_v_field_raw,
2838  const GriddedField3& wind_w_field_raw,
2839  const Index& atmosphere_dim,
2840  // WS Generic Input:
2841  const Index& interp_order,
2842  const Verbosity& verbosity) {
2843  CREATE_OUT2;
2844 
2845  const Vector& ufr_p_grid =
2846  wind_u_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2847  const Vector& ufr_lat_grid =
2848  wind_u_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2849  const Vector& ufr_lon_grid =
2850  wind_u_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2851  const Vector& vfr_p_grid =
2852  wind_v_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2853  const Vector& vfr_lat_grid =
2854  wind_v_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2855  const Vector& vfr_lon_grid =
2856  wind_v_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2857  const Vector& wfr_p_grid =
2858  wind_w_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2859  const Vector& wfr_lat_grid =
2860  wind_w_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2861  const Vector& wfr_lon_grid =
2862  wind_w_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2863 
2864  out2 << " Interpolation order: " << interp_order << "\n";
2865 
2866  // Basic checks of input variables
2867  //
2868  // Atmosphere
2869  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2870  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
2871 
2872  //==========================================================================
2873  if (atmosphere_dim == 1) {
2874  if (!(ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1))
2875  throw std::runtime_error(
2876  "Wind u field data has wrong dimension (2D or 3D).\n");
2877  if (!(vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1))
2878  throw std::runtime_error(
2879  "Wind v field data has wrong dimension (2D or 3D).\n");
2880  if (!(wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1))
2881  throw std::runtime_error(
2882  "Wind w field data has wrong dimension (2D or 3D).\n");
2883 
2884  GriddedField3 temp_gfield3;
2885 
2887  temp_gfield3, p_grid, wind_u_field_raw, interp_order, 0, verbosity);
2888  wind_u_field = temp_gfield3.data;
2889 
2891  temp_gfield3, p_grid, wind_v_field_raw, interp_order, 0, verbosity);
2892  wind_v_field = temp_gfield3.data;
2893 
2895  temp_gfield3, p_grid, wind_w_field_raw, interp_order, 0, verbosity);
2896  wind_w_field = temp_gfield3.data;
2897  }
2898 
2899  //=========================================================================
2900  else if (atmosphere_dim == 2) {
2901  if (ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1)
2902  throw std::runtime_error(
2903  "Raw data has wrong dimension (1D). You have to use \n"
2904  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
2905  if (vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1)
2906  throw std::runtime_error(
2907  "Raw data has wrong dimension (1D). You have to use \n"
2908  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
2909  if (wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1)
2910  throw std::runtime_error(
2911  "Raw data has wrong dimension (1D). You have to use \n"
2912  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
2913 
2914  //Resize variables
2915  wind_u_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2916  wind_v_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2917  wind_w_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
2918 
2919  // Gridpositions:
2920  ArrayOfGridPosPoly gp_p(p_grid.nelem());
2921  ArrayOfGridPosPoly gp_lat(lat_grid.nelem());
2922 
2923  // Interpolate wind_u_field:
2924 
2925  // Check that interpolation grids are ok (and throw a detailed
2926  // error message if not):
2928  "Raw u field to p_grid, 2D case", ufr_p_grid, p_grid, interp_order);
2929  chk_interpolation_grids("Raw u field to lat_grid, 2D case",
2930  ufr_lat_grid,
2931  lat_grid,
2932  interp_order);
2933 
2934  // Calculate grid positions:
2935  p2gridpos_poly(gp_p, ufr_p_grid, p_grid, interp_order);
2936  gridpos_poly(gp_lat, ufr_lat_grid, lat_grid, interp_order);
2937 
2938  // Interpolation weights:
2939  Tensor3 itwu(p_grid.nelem(),
2940  lat_grid.nelem(),
2941  (interp_order + 1) * (interp_order + 1));
2942  // (4 interpolation weights are required for example for linear 2D
2943  // interpolation)
2944  interpweights(itwu, gp_p, gp_lat);
2945 
2946  // Interpolate:
2947  interp(wind_u_field(joker, joker, 0),
2948  itwu,
2949  wind_u_field_raw.data(joker, joker, 0),
2950  gp_p,
2951  gp_lat);
2952 
2953  // Interpolate wind_v_field:
2954 
2955  // Check that interpolation grids are ok (and throw a detailed
2956  // error message if not):
2958  "Raw v field to p_grid, 2D case", vfr_p_grid, p_grid, interp_order);
2959  chk_interpolation_grids("Raw v field to lat_grid, 2D case",
2960  vfr_lat_grid,
2961  lat_grid,
2962  interp_order);
2963 
2964  // Calculate grid positions:
2965  p2gridpos_poly(gp_p, vfr_p_grid, p_grid, interp_order);
2966  gridpos_poly(gp_lat, vfr_lat_grid, lat_grid, interp_order);
2967 
2968  // Interpolation weights:
2969  Tensor3 itwv(p_grid.nelem(),
2970  lat_grid.nelem(),
2971  (interp_order + 1) * (interp_order + 1));
2972  // (4 interpolation weights are required for example for linear 2D
2973  // interpolation)
2974  interpweights(itwv, gp_p, gp_lat);
2975 
2976  // Interpolate:
2977  interp(wind_v_field(joker, joker, 0),
2978  itwv,
2979  wind_v_field_raw.data(joker, joker, 0),
2980  gp_p,
2981  gp_lat);
2982 
2983  // Interpolate wind_w_field:
2984 
2985  // Check that interpolation grids are ok (and throw a detailed
2986  // error message if not):
2988  "Raw w field to p_grid, 2D case", wfr_p_grid, p_grid, interp_order);
2989  chk_interpolation_grids("Raw w field to lat_grid, 2D case",
2990  wfr_lat_grid,
2991  lat_grid,
2992  interp_order);
2993 
2994  // Calculate grid positions:
2995  p2gridpos_poly(gp_p, wfr_p_grid, p_grid, interp_order);
2996  gridpos_poly(gp_lat, wfr_lat_grid, lat_grid, interp_order);
2997 
2998  // Interpolation weights:
2999  Tensor3 itww(p_grid.nelem(),
3000  lat_grid.nelem(),
3001  (interp_order + 1) * (interp_order + 1));
3002  // (4 interpolation weights are required for example for linear 2D
3003  // interpolation)
3004  interpweights(itww, gp_p, gp_lat);
3005 
3006  // Interpolate:
3007  interp(wind_w_field(joker, joker, 0),
3008  itww,
3009  wind_w_field_raw.data(joker, joker, 0),
3010  gp_p,
3011  gp_lat);
3012  }
3013 
3014  //================================================================
3015  // atmosphere_dim = 3
3016  else if (atmosphere_dim == 3) {
3017  if (ufr_lat_grid.nelem() == 1 && ufr_lon_grid.nelem() == 1)
3018  throw std::runtime_error(
3019  "Raw data has wrong dimension. You have to use \n"
3020  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
3021  if (vfr_lat_grid.nelem() == 1 && vfr_lon_grid.nelem() == 1)
3022  throw std::runtime_error(
3023  "Raw data has wrong dimension. You have to use \n"
3024  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
3025  if (wfr_lat_grid.nelem() == 1 && wfr_lon_grid.nelem() == 1)
3026  throw std::runtime_error(
3027  "Raw data has wrong dimension. You have to use \n"
3028  "WindFieldsCalcExpand1D instead of WindFieldsCalc.");
3029 
3030  GriddedField3 temp_gfield3;
3031 
3032  GriddedFieldLatLonRegrid(temp_gfield3,
3033  lat_grid,
3034  lon_grid,
3035  wind_u_field_raw,
3036  interp_order,
3037  verbosity);
3039  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
3040  wind_u_field = temp_gfield3.data;
3041 
3042  GriddedFieldLatLonRegrid(temp_gfield3,
3043  lat_grid,
3044  lon_grid,
3045  wind_v_field_raw,
3046  interp_order,
3047  verbosity);
3049  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
3050  wind_v_field = temp_gfield3.data;
3051 
3052  GriddedFieldLatLonRegrid(temp_gfield3,
3053  lat_grid,
3054  lon_grid,
3055  wind_w_field_raw,
3056  interp_order,
3057  verbosity);
3059  temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
3060  wind_w_field = temp_gfield3.data;
3061 
3062  } else {
3063  // We can never get here, since there was a runtime
3064  // error check for atmosphere_dim at the beginning.
3065  assert(false);
3066  }
3067 }
3068 
3069 /* Workspace method: Doxygen documentation will be auto-generated */
3071  Tensor3& z_field,
3072  Tensor4& vmr_field,
3073  EnergyLevelMap& nlte_field,
3074  const Vector& p_grid,
3075  const Vector& lat_grid,
3076  const Vector& lon_grid,
3077  const GriddedField3& t_field_raw,
3078  const GriddedField3& z_field_raw,
3079  const ArrayOfGriddedField3& vmr_field_raw,
3080  const ArrayOfGriddedField3& nlte_field_raw,
3081  const ArrayOfQuantumIdentifier& nlte_ids,
3082  const Vector& nlte_energies,
3083  const Index& atmosphere_dim,
3084  const Index& interp_order,
3085  const Index& vmr_zeropadding,
3086  const Index& vmr_nonegative,
3087  const Index& nlte_when_negative,
3088  const Verbosity& verbosity) {
3089  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3090  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3091 
3092  if (atmosphere_dim == 1) {
3093  throw runtime_error(
3094  "This function is intended for 2D and 3D. For 1D, use *AtmFieldsCalc*.");
3095  }
3096 
3097  // Make 1D interpolation using some dummy variables
3098  Vector vempty(0);
3099  Tensor3 t_temp, z_temp;
3100  Tensor4 vmr_temp;
3101  EnergyLevelMap nlte_temp;
3102  AtmFieldsCalc(t_temp,
3103  z_temp,
3104  vmr_temp,
3105  nlte_temp,
3106  p_grid,
3107  vempty,
3108  vempty,
3109  t_field_raw,
3110  z_field_raw,
3111  vmr_field_raw,
3112  nlte_field_raw,
3113  nlte_ids,
3114  nlte_energies,
3115  1,
3116  interp_order,
3117  vmr_zeropadding,
3118  vmr_nonegative,
3119  nlte_when_negative,
3120  verbosity);
3121 
3122  // Move values from the temporary tensors to the return arguments
3123  const Index np = p_grid.nelem();
3124  const Index nlat = lat_grid.nelem();
3125  Index nlon = lon_grid.nelem();
3126  if (atmosphere_dim == 2) {
3127  nlon = 1;
3128  }
3129  const Index nspecies = vmr_temp.nbooks();
3130  //
3131  assert(t_temp.npages() == np);
3132  //
3133  t_field.resize(np, nlat, nlon);
3134  z_field.resize(np, nlat, nlon);
3135  vmr_field.resize(nspecies, np, nlat, nlon);
3136  if (nlte_field_raw.nelem()) {
3137  nlte_field.Type() = EnergyLevelMapType::Tensor3_t;
3138  nlte_field.Data().resize(nlte_field_raw.nelem(), np, nlat, nlon);
3139  nlte_field.Levels() = nlte_ids;
3140  nlte_field.Energies() = nlte_energies;
3141  }
3142  else
3143  nlte_field = EnergyLevelMap();
3144  //
3145  for (Index ilon = 0; ilon < nlon; ilon++) {
3146  for (Index ilat = 0; ilat < nlat; ilat++) {
3147  for (Index ip = 0; ip < np; ip++) {
3148  t_field(ip, ilat, ilon) = t_temp(ip, 0, 0);
3149  z_field(ip, ilat, ilon) = z_temp(ip, 0, 0);
3150  for (Index is = 0; is < nspecies; is++) {
3151  vmr_field(is, ip, ilat, ilon) = vmr_temp(is, ip, 0, 0);
3152  }
3153  for (Index is = 0; is < nlte_field_raw.nelem(); is++) {
3154  nlte_field.Data()(is, ip, ilat, ilon) = nlte_temp.Data()(is, ip, 0, 0);
3155  }
3156  }
3157  }
3158  }
3159  nlte_field.ThrowIfNotOK();
3160 }
3161 
3162 /* Workspace method: Doxygen documentation will be auto-generated */
3163 void MagFieldsCalcExpand1D(Tensor3& mag_u_field,
3164  Tensor3& mag_v_field,
3165  Tensor3& mag_w_field,
3166  const Vector& p_grid,
3167  const Vector& lat_grid,
3168  const Vector& lon_grid,
3169  const GriddedField3& mag_u_field_raw,
3170  const GriddedField3& mag_v_field_raw,
3171  const GriddedField3& mag_w_field_raw,
3172  const Index& atmosphere_dim,
3173  const Index& interp_order,
3174  const Verbosity& verbosity) {
3175  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3176  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3177 
3178  if (atmosphere_dim == 1)
3179  throw std::runtime_error(
3180  "This function is intended for 2D and 3D. For 1D, use *MagFieldsCalc*.");
3181 
3182  // Make 1D interpolation using some dummy variables
3183  Vector vempty(0);
3184  Tensor3 mag_u_field_temp, mag_v_field_temp, mag_w_field_temp;
3185  MagFieldsCalc(mag_u_field_temp,
3186  mag_v_field_temp,
3187  mag_w_field_temp,
3188  p_grid,
3189  vempty,
3190  vempty,
3191  mag_u_field_raw,
3192  mag_v_field_raw,
3193  mag_w_field_raw,
3194  /*atmosphere_dim = */ 1,
3195  interp_order,
3196  verbosity);
3197 
3198  // Move values from the temporary tensors to the return arguments
3199  const Index np = p_grid.nelem();
3200  const Index nlat = lat_grid.nelem();
3201  Index nlon = lon_grid.nelem();
3202  if (atmosphere_dim == 2) {
3203  nlon = 1;
3204  }
3205  //
3206  assert(mag_u_field_temp.npages() == np);
3207  assert(mag_v_field_temp.npages() == np);
3208  assert(mag_w_field_temp.npages() == np);
3209  //
3210  mag_u_field.resize(np, nlat, nlon);
3211  mag_v_field.resize(np, nlat, nlon);
3212  mag_w_field.resize(np, nlat, nlon);
3213 
3214  for (Index ilon = 0; ilon < nlon; ilon++) {
3215  for (Index ilat = 0; ilat < nlat; ilat++) {
3216  for (Index ip = 0; ip < np; ip++) {
3217  mag_u_field(ip, ilat, ilon) = mag_u_field_temp(ip, 0, 0);
3218  mag_v_field(ip, ilat, ilon) = mag_v_field_temp(ip, 0, 0);
3219  mag_w_field(ip, ilat, ilon) = mag_w_field_temp(ip, 0, 0);
3220  }
3221  }
3222  }
3223 }
3224 
3225 /* Workspace method: Doxygen documentation will be auto-generated */
3226 void WindFieldsCalcExpand1D(Tensor3& wind_u_field,
3227  Tensor3& wind_v_field,
3228  Tensor3& wind_w_field,
3229  const Vector& p_grid,
3230  const Vector& lat_grid,
3231  const Vector& lon_grid,
3232  const GriddedField3& wind_u_field_raw,
3233  const GriddedField3& wind_v_field_raw,
3234  const GriddedField3& wind_w_field_raw,
3235  const Index& atmosphere_dim,
3236  const Index& interp_order,
3237  const Verbosity& verbosity) {
3238  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3239  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3240 
3241  if (atmosphere_dim == 1)
3242  throw std::runtime_error(
3243  "This function is intended for 2D and 3D. For 1D, use *WindFieldsCalc*.");
3244 
3245  // Make 1D interpolation using some dummy variables
3246  Vector vempty(0);
3247  Tensor3 wind_u_field_temp, wind_v_field_temp, wind_w_field_temp;
3248  MagFieldsCalc(wind_u_field_temp,
3249  wind_v_field_temp,
3250  wind_w_field_temp,
3251  p_grid,
3252  vempty,
3253  vempty,
3254  wind_u_field_raw,
3255  wind_v_field_raw,
3256  wind_w_field_raw,
3257  /*atmosphere_dim = */ 1,
3258  interp_order,
3259  verbosity);
3260 
3261  // Move values from the temporary tensors to the return arguments
3262  const Index np = p_grid.nelem();
3263  const Index nlat = lat_grid.nelem();
3264  Index nlon = lon_grid.nelem();
3265  if (atmosphere_dim == 2) {
3266  nlon = 1;
3267  }
3268  //
3269  assert(wind_u_field_temp.npages() == np);
3270  assert(wind_v_field_temp.npages() == np);
3271  assert(wind_w_field_temp.npages() == np);
3272  //
3273  wind_u_field.resize(np, nlat, nlon);
3274  wind_v_field.resize(np, nlat, nlon);
3275  wind_w_field.resize(np, nlat, nlon);
3276 
3277  for (Index ilon = 0; ilon < nlon; ilon++) {
3278  for (Index ilat = 0; ilat < nlat; ilat++) {
3279  for (Index ip = 0; ip < np; ip++) {
3280  wind_u_field(ip, ilat, ilon) = wind_u_field_temp(ip, 0, 0);
3281  wind_v_field(ip, ilat, ilon) = wind_v_field_temp(ip, 0, 0);
3282  wind_w_field(ip, ilat, ilon) = wind_w_field_temp(ip, 0, 0);
3283  }
3284  }
3285  }
3286 }
3287 
3288 /* Workspace method: Doxygen documentation will be auto-generated */
3290  Tensor3& z_field,
3291  Tensor4& vmr_field,
3292  const Vector& p_grid,
3293  const Vector& lat_grid,
3294  const Vector& lon_grid,
3295  const Index& atmosphere_dim,
3296  const Index& chk_vmr_nan,
3297  const Verbosity&) {
3298  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
3299  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3300 
3301  // Sizes
3302  const Index np = p_grid.nelem();
3303  const Index nlat = lat_grid.nelem();
3304  const Index nlon = max(Index(1), lon_grid.nelem());
3305  const Index nspecies = vmr_field.nbooks();
3306 
3307  const bool chknan = chk_vmr_nan;
3308 
3309  if (atmosphere_dim == 1) {
3310  throw runtime_error("No use in calling this method for 1D.");
3311  }
3312  chk_atm_field("t_field", t_field, 1, p_grid, Vector(0), Vector(0));
3313  chk_atm_field("z_field", z_field, 1, p_grid, Vector(0), Vector(0));
3314  if (nspecies)
3315  chk_atm_field("vmr_field",
3316  vmr_field,
3317  1,
3318  nspecies,
3319  p_grid,
3320  Vector(0),
3321  Vector(0),
3322  chknan);
3323 
3324  // Temporary containers
3325  Tensor3 t_temp = t_field, z_temp = z_field;
3326  Tensor4 vmr_temp = vmr_field;
3327 
3328  // Resize and fill
3329  t_field.resize(np, nlat, nlon);
3330  z_field.resize(np, nlat, nlon);
3331  vmr_field.resize(nspecies, np, nlat, nlon);
3332  //
3333  for (Index ilon = 0; ilon < nlon; ilon++) {
3334  for (Index ilat = 0; ilat < nlat; ilat++) {
3335  for (Index ip = 0; ip < np; ip++) {
3336  t_field(ip, ilat, ilon) = t_temp(ip, 0, 0);
3337  z_field(ip, ilat, ilon) = z_temp(ip, 0, 0);
3338  for (Index is = 0; is < nspecies; is++) {
3339  vmr_field(is, ip, ilat, ilon) = vmr_temp(is, ip, 0, 0);
3340  }
3341  }
3342  }
3343  }
3344 }
3345 
3346 /* Workspace method: Doxygen documentation will be auto-generated */
3347 void AtmFieldsExtract1D(Index& atmosphere_dim,
3348  Vector& lat_grid,
3349  Vector& lon_grid,
3350  Tensor3& t_field,
3351  Tensor3& z_field,
3352  Tensor4& vmr_field,
3353  const Index& ilat,
3354  const Index& ilon,
3355  const Verbosity& verbosity) {
3356  if (atmosphere_dim == 1) {
3357  return;
3358  }
3359 
3360  if (ilat < 0 || ilat >= lat_grid.nelem())
3361  throw runtime_error(
3362  "Invalid of *ilat*. It must be >= 0 and less than "
3363  "length of *lat_grid*.");
3364 
3365  if (atmosphere_dim == 2) {
3366  Vector vtmp;
3367  vtmp = t_field(joker, ilat, 0);
3368  t_field.resize(t_field.npages(), 1, 1);
3369  t_field(joker, 0, 0) = vtmp;
3370  vtmp = z_field(joker, ilat, 0);
3371  z_field.resize(z_field.npages(), 1, 1);
3372  z_field(joker, 0, 0) = vtmp;
3373  Matrix mtmp;
3374  mtmp = vmr_field(joker, joker, ilat, 0);
3375  vmr_field.resize(vmr_field.nbooks(), vmr_field.npages(), 1, 1);
3376  vmr_field(joker, joker, 0, 0) = mtmp;
3377  } else if (atmosphere_dim == 3) {
3378  if (ilat < 0 || ilon >= lon_grid.nelem())
3379  throw runtime_error(
3380  "Invalid of *ilon*. It must be >= 0 and less than "
3381  "length of *lon_grid*.");
3382  Vector vtmp;
3383  vtmp = t_field(joker, ilat, ilon);
3384  t_field.resize(t_field.npages(), 1, 1);
3385  t_field(joker, 0, 0) = vtmp;
3386  vtmp = z_field(joker, ilat, ilon);
3387  z_field.resize(z_field.npages(), 1, 1);
3388  z_field(joker, 0, 0) = vtmp;
3389  Matrix mtmp;
3390  mtmp = vmr_field(joker, joker, ilat, ilon);
3391  vmr_field.resize(vmr_field.nbooks(), vmr_field.npages(), 1, 1);
3392  vmr_field(joker, joker, 0, 0) = mtmp;
3393  }
3394 
3395  else {
3396  throw runtime_error("Invalid of *atmosphere_dim*. It must be 1-3.");
3397  }
3398 
3399  AtmosphereSet1D(atmosphere_dim, lat_grid, lon_grid, verbosity);
3400 }
3401 
3402 /* Workspace method: Doxygen documentation will be auto-generated */
3403 void AtmFieldsRefinePgrid( // WS Output:
3404  Vector& p_grid,
3405  Tensor3& t_field,
3406  Tensor3& z_field,
3407  Tensor4& vmr_field,
3408  Index& atmfields_checked,
3409  Index& atmgeom_checked,
3410  Index& cloudbox_checked,
3411  // WS Input:
3412  const Vector& lat_grid,
3413  const Vector& lon_grid,
3414  const Index& atmosphere_dim,
3415  // Control Parameters:
3416  const Numeric& p_step,
3417  const Index& interp_order,
3418  const Verbosity& verbosity) {
3419  // Checks on input parameters:
3420 
3421  // We don't actually need lat_grid and lon_grid, but we have them as
3422  // input variables, so that we can use the standard functions to
3423  // check atmospheric fields and grids. A bit cheesy, but I don't
3424  // want to program all the checks explicitly.
3425 
3426  // Check grids:
3427  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3428 
3429  // Check T field:
3430  chk_atm_field("t_field", t_field, atmosphere_dim, p_grid, lat_grid, lon_grid);
3431 
3432  // Check z field:
3433  chk_atm_field("z_field", z_field, atmosphere_dim, p_grid, lat_grid, lon_grid);
3434 
3435  // Check VMR field (and abs_species):
3436  chk_atm_field("vmr_field",
3437  vmr_field,
3438  atmosphere_dim,
3439  vmr_field.nbooks(),
3440  p_grid,
3441  lat_grid,
3442  lon_grid);
3443 
3444  // Move original p_grid to p_old, freeing p_grid for the refined one.
3445  Vector p_old;
3446  p_old = p_grid;
3447 
3448  p_gridRefine(p_grid,
3449  atmfields_checked,
3450  atmgeom_checked,
3451  cloudbox_checked,
3452  p_old,
3453  p_step,
3454  verbosity);
3455 
3456  AtmFieldPRegrid(z_field, z_field, p_grid, p_old, interp_order, verbosity);
3457  AtmFieldPRegrid(t_field, t_field, p_grid, p_old, interp_order, verbosity);
3458  AtmFieldPRegrid(vmr_field, vmr_field, p_grid, p_old, interp_order, verbosity);
3459 }
3460 
3461 /* Workspace method: Doxygen documentation will be auto-generated */
3462 void AtmRawRead( //WS Output:
3463  GriddedField3& t_field_raw,
3464  GriddedField3& z_field_raw,
3465  ArrayOfGriddedField3& vmr_field_raw,
3466  ArrayOfGriddedField3& nlte_field_raw,
3467  ArrayOfQuantumIdentifier& nlte_quantum_identifiers,
3468  Vector& nlte_vibrational_energies,
3469  //WS Input:
3470  const ArrayOfArrayOfSpeciesTag& abs_species,
3471  //Keyword:
3472  const String& basename,
3473  const Verbosity& verbosity) {
3474  CREATE_OUT3;
3475 
3476  String tmp_basename = basename;
3477  if (basename.length() && basename[basename.length() - 1] != '/')
3478  tmp_basename += ".";
3479 
3480  // Read the temperature field:
3481  String file_name = tmp_basename + "t.xml";
3482  xml_read_from_file(file_name, t_field_raw, verbosity);
3483 
3484  out3 << "Temperature field read from file: " << file_name << "\n";
3485 
3486  // Read geometrical altitude field:
3487  file_name = tmp_basename + "z.xml";
3488  xml_read_from_file(file_name, z_field_raw, verbosity);
3489 
3490  out3 << "Altitude field read from file: " << file_name << "\n";
3491 
3492  // Clear out vmr_field_raw
3493  vmr_field_raw.resize(0);
3494 
3495  // The species lookup data:
3497 
3498  // We need to read one profile for each tag group.
3499  for (Index i = 0; i < abs_species.nelem(); i++) {
3500  // Determine the name.
3501  file_name = tmp_basename +
3502  species_data[abs_species[i][0].Species()].Name() + ".xml";
3503 
3504  // Add an element for this tag group to the vmr profiles:
3505  GriddedField3 vmr_field_data;
3506  vmr_field_raw.push_back(vmr_field_data);
3507 
3508  // Read the VMR:
3510  file_name, vmr_field_raw[vmr_field_raw.nelem() - 1], verbosity);
3511 
3512  // state the source of profile.
3513  out3 << " " << species_data[abs_species[i][0].Species()].Name()
3514  << " profile read from file: " << file_name << "\n";
3515  }
3516 
3517  // NLTE is ignored by doing this
3518  nlte_field_raw.resize(0);
3519  nlte_quantum_identifiers.resize(0);
3520  nlte_vibrational_energies.resize(0);
3521 }
3522 
3523 /* Workspace method: Doxygen documentation will be auto-generated */
3524 void MagRawRead( //WS Output:
3525  GriddedField3& mag_u_field_raw,
3526  GriddedField3& mag_v_field_raw,
3527  GriddedField3& mag_w_field_raw,
3528  //Keyword:
3529  const String& basename,
3530  const Verbosity& verbosity) {
3531  CREATE_OUT3;
3532 
3533  String tmp_basename = basename;
3534  if (basename.length() && basename[basename.length() - 1] != '/')
3535  tmp_basename += ".";
3536 
3537  // Read magfield u component:
3538  String file_name = tmp_basename + "mag_u.xml";
3539  xml_read_from_file(file_name, mag_u_field_raw, verbosity);
3540 
3541  out3 << "Bu field read from file: " << file_name << "\n";
3542 
3543  // Read magfield v component:
3544  file_name = tmp_basename + "mag_v.xml";
3545  xml_read_from_file(file_name, mag_v_field_raw, verbosity);
3546 
3547  out3 << "Bv field read from file: " << file_name << "\n";
3548 
3549  // Read magfield w component:
3550  file_name = tmp_basename + "mag_w.xml";
3551  xml_read_from_file(file_name, mag_w_field_raw, verbosity);
3552 
3553  out3 << "Bw field read from file: " << file_name << "\n";
3554 }
3555 
3556 /* Workspace method: Doxygen documentation will be auto-generated */
3557 void WindRawRead( //WS Output:
3558  GriddedField3& wind_u_field_raw,
3559  GriddedField3& wind_v_field_raw,
3560  GriddedField3& wind_w_field_raw,
3561  //Keyword:
3562  const String& basename,
3563  const Verbosity& verbosity) {
3564  CREATE_OUT3;
3565 
3566  String tmp_basename = basename;
3567  if (basename.length() && basename[basename.length() - 1] != '/')
3568  tmp_basename += ".";
3569 
3570  // Read wind field u component:
3571  String file_name = tmp_basename + "wind_u.xml";
3572  xml_read_from_file(file_name, wind_u_field_raw, verbosity);
3573 
3574  out3 << "Wind u field read from file: " << file_name << "\n";
3575 
3576  // Read wind field u component:
3577  file_name = tmp_basename + "wind_v.xml";
3578  xml_read_from_file(file_name, wind_v_field_raw, verbosity);
3579 
3580  out3 << "Wind v field read from file: " << file_name << "\n";
3581 
3582  // Read wind field u component:
3583  file_name = tmp_basename + "wind_w.xml";
3584  xml_read_from_file(file_name, wind_w_field_raw, verbosity);
3585 
3586  out3 << "Wind w field read from file: " << file_name << "\n";
3587 }
3588 
3589 /* Workspace method: Doxygen documentation will be auto-generated */
3590 void AtmWithNLTERawRead( //WS Output:
3591  GriddedField3& t_field_raw,
3592  GriddedField3& z_field_raw,
3593  ArrayOfGriddedField3& vmr_field_raw,
3594  ArrayOfGriddedField3& nlte_field_raw,
3595  ArrayOfQuantumIdentifier& nlte_quantum_identifiers,
3596  Vector& nlte_vibrational_energies,
3597  //WS Input:
3598  const ArrayOfArrayOfSpeciesTag& abs_species,
3599  //Keyword:
3600  const String& basename,
3601  const Index& expect_vibrational_energies,
3602  const Verbosity& verbosity) {
3603  CREATE_OUT3;
3604 
3605  String tmp_basename = basename;
3606  if (basename.length() && basename[basename.length() - 1] != '/')
3607  tmp_basename += ".";
3608 
3609  // Read the temperature field:
3610  String file_name = tmp_basename + "t.xml";
3611  xml_read_from_file(file_name, t_field_raw, verbosity);
3612 
3613  out3 << "Temperature field read from file: " << file_name << "\n";
3614 
3615  // Read geometrical altitude field:
3616  file_name = tmp_basename + "z.xml";
3617  xml_read_from_file(file_name, z_field_raw, verbosity);
3618 
3619  out3 << "Altitude field read from file: " << file_name << "\n";
3620 
3621  // Clear out vmr_field_raw
3622  vmr_field_raw.resize(0);
3623 
3624  // The species lookup data:
3626 
3627  // We need to read one profile for each tag group.
3628  for (Index i = 0; i < abs_species.nelem(); i++) {
3629  // Determine the name.
3630  file_name = tmp_basename +
3631  species_data[abs_species[i][0].Species()].Name() + ".xml";
3632 
3633  // Add an element for this tag group to the vmr profiles:
3634  GriddedField3 vmr_field_data;
3635  vmr_field_raw.push_back(vmr_field_data);
3636 
3637  // Read the VMR:
3639  file_name, vmr_field_raw[vmr_field_raw.nelem() - 1], verbosity);
3640 
3641  // state the source of profile.
3642  out3 << " " << species_data[abs_species[i][0].Species()].Name()
3643  << " profile read from file: " << file_name << "\n";
3644  }
3645 
3646  // Read each nlte field:
3647  file_name = tmp_basename + "nlte.xml";
3648  xml_read_from_file(file_name, nlte_field_raw, verbosity);
3649 
3650  out3 << "NLTE field array read from file: " << file_name << "\n";
3651 
3652  // Read each nlte identifier field:
3653  file_name = tmp_basename + "qi.xml";
3654  xml_read_from_file(file_name, nlte_quantum_identifiers, verbosity);
3655 
3656  out3 << "NLTE identifier array read from file: " << file_name << "\n";
3657 
3658  if (expect_vibrational_energies) {
3659  // Read each energy level field:
3660  file_name = tmp_basename + "ev.xml";
3661  xml_read_from_file(file_name, nlte_vibrational_energies, verbosity);
3662 
3663  out3 << "NLTE energy levels array read from file: " << file_name << "\n";
3664  }
3665  else {
3666  nlte_vibrational_energies.resize(0);
3667  }
3668 
3669  if (nlte_field_raw.nelem() != nlte_quantum_identifiers.nelem() or
3670  (nlte_field_raw.nelem() != nlte_vibrational_energies.nelem() and
3671  0 != nlte_vibrational_energies.nelem())) {
3672  ostringstream os;
3673  os << "The quantum identifers and the NLTE temperature fields\n"
3674  << "are of different lengths. This should not be the case.\n"
3675  << "please check the qi.xml and t_nlte.xml files under\n"
3676  << basename << std::endl
3677  << "to correct this error.\n";
3678  throw std::runtime_error(os.str());
3679  }
3680 }
3681 
3682 /* Workspace method: Doxygen documentation will be auto-generated */
3684  const Vector& lat_grid,
3685  const Vector& lon_grid,
3686  const String& filename,
3687  const Index& interp_order,
3688  const Index& set_lowest_altitude_to_zero,
3689  const Verbosity& verbosity) {
3690  CREATE_OUT3;
3691 
3692  out3 << "Reading GriddedField2 surface altitude from " << filename << "\n";
3693  GriddedField2 z_surface_field;
3694  xml_read_from_file(filename, z_surface_field, verbosity);
3695 
3696  out3 << "Surface altitude field interpolated back to lat_grid and lon_grid\n";
3697  GriddedFieldLatLonRegrid(z_surface_field,
3698  lat_grid,
3699  lon_grid,
3700  z_surface_field,
3701  interp_order,
3702  verbosity);
3703  z_surface = z_surface_field.data;
3704  if (set_lowest_altitude_to_zero) {
3705  z_surface -= min(z_surface);
3706  }
3707 }
3708 
3709 /* Workspace method: Doxygen documentation will be auto-generated */
3711  const Vector& lat_grid,
3712  const Vector& lon_grid,
3713  const Numeric& altitude,
3714  const Verbosity& verbosity) {
3715  CREATE_OUT3;
3716  out3 << "Setting surface to constant altitude of " << altitude << " m\n";
3717  z_surface = Matrix(lat_grid.nelem() ? lat_grid.nelem() : 1, lon_grid.nelem() ? lon_grid.nelem() : 1, altitude);
3718 }
3719 
3720 /* Workspace method: Doxygen documentation will be auto-generated */
3722  const Index& atmosphere_dim,
3723  const Vector& p_grid,
3724  const Vector& lat_grid,
3725  const Vector& lon_grid,
3726  const Tensor3& z_field,
3727  const Vector& rtp_pos,
3728  const Tensor3& field,
3729  const Verbosity& verbosity) {
3730  // Input checks
3731  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
3732  chk_atm_field("input argument *field*",
3733  field,
3734  atmosphere_dim,
3735  p_grid,
3736  lat_grid,
3737  lon_grid);
3738  chk_rte_pos(atmosphere_dim, rtp_pos);
3739 
3740  // Determine grid positions
3741  GridPos gp_p, gp_lat, gp_lon;
3742  rte_pos2gridpos(gp_p,
3743  gp_lat,
3744  gp_lon,
3745  atmosphere_dim,
3746  p_grid,
3747  lat_grid,
3748  lon_grid,
3749  z_field,
3750  rtp_pos);
3751 
3752  // Interpolate
3753  outvalue = interp_atmfield_by_gp(atmosphere_dim, field, gp_p, gp_lat, gp_lon);
3754 
3755  CREATE_OUT3;
3756  out3 << " Result = " << outvalue << "\n";
3757 }
3758 
3759 /* Workspace method: Doxygen documentation will be auto-generated */
3760 void p_gridDensify( // WS Output:
3761  Vector& p_grid,
3762  Index& atmfields_checked,
3763  Index& atmgeom_checked,
3764  Index& cloudbox_checked,
3765  // WS Input:
3766  const Vector& p_grid_old,
3767  // Control Parameters:
3768  const Index& nfill,
3769  const Verbosity& verbosity) {
3770  // Check that p_grid and p_grid_old are not the same variable (pointing to the
3771  // same memory space). this as p_grid will be overwritten, but we will need
3772  // both data later on for data regridding.
3773  if (&p_grid == &p_grid_old) {
3774  ostringstream os;
3775  os << "The old and new grids (p_grid and p_grid_old) are not allowed\n"
3776  << "to be identical (pointing to same memory space).\n"
3777  << "But they are doing in your case.";
3778  throw runtime_error(os.str());
3779  }
3780 
3781  // as we manipoulate the overall vertical grid (but not simultaneously the
3782  // atmospheric fields), we reset all atmfields related checked WSV to
3783  // unchecked, forcing the user to do the checks again.
3784  atmfields_checked = 0;
3785  atmgeom_checked = 0;
3786  cloudbox_checked = 0;
3787 
3788  // Check the keyword argument:
3789  if (nfill < 0) {
3790  throw runtime_error("Argument *nfill* must be >= 0.");
3791  }
3792 
3793  // Nothing to do if nfill=0
3794  if (nfill > 0) {
3795  // Allocate new size for p_grid
3796  const Index n0 = p_grid_old.nelem();
3797  p_grid.resize((n0 - 1) * (1 + nfill) + 1);
3798 
3799  Index iout = 0;
3800  p_grid[0] = p_grid_old[0];
3801 
3802  for (Index i = 1; i < n0; i++) {
3803  Vector pnew;
3805  pnew, 2 + nfill, p_grid_old[i - 1], p_grid_old[i], verbosity);
3806  for (Index j = 1; j < nfill + 2; j++) {
3807  iout += 1;
3808  p_grid[iout] = pnew[j];
3809  }
3810  }
3811  }
3812 }
3813 
3814 /* Workspace method: Doxygen documentation will be auto-generated */
3815 void p_gridRefine( // WS Output:
3816  Vector& p_grid,
3817  Index& atmfields_checked,
3818  Index& atmgeom_checked,
3819  Index& cloudbox_checked,
3820  // WS Input:
3821  const Vector& p_grid_old,
3822  // Control Parameters:
3823  const Numeric& p_step10,
3824  const Verbosity&) {
3825  // Check that p_grid and p_grid_old are not the same variable (pointing to the
3826  // same memory space). this as p_grid will be overwritten, but we will need
3827  // both data later on for data regridding.
3828  if (&p_grid == &p_grid_old) {
3829  ostringstream os;
3830  os << "The old and new grids (p_grid and p_grid_old) are not allowed\n"
3831  << "to be identical (pointing to same memory space).\n"
3832  << "But they are doing in your case.";
3833  throw runtime_error(os.str());
3834  }
3835 
3836  // as we manipoulate the overall vertical grid (but not simultaneously the
3837  // atmospheric fields), we reset all atmfields related checked WSV to
3838  // unchecked, forcing the user to do the checks again.
3839  atmfields_checked = 0;
3840  atmgeom_checked = 0;
3841  cloudbox_checked = 0;
3842 
3843  // Check the keyword argument:
3844  if (p_step10 <= 0) {
3845  ostringstream os;
3846  os << "The keyword argument p_step must be >0.";
3847  throw runtime_error(os.str());
3848  }
3849 
3850  // For consistency with other code around arts (e.g., correlation
3851  // lengths in atmlab), p_step is given as log10(p[Pa]). However, we
3852  // convert it here to the natural log:
3853  const Numeric p_step = log(pow(10.0, p_step10));
3854 
3855  // Now starting modification of p_grid
3856 
3857  // We will need the log of the pressure grid:
3858  Vector log_p_old(p_grid_old.nelem());
3859  transform(log_p_old, log, p_grid_old);
3860 
3861  // const Numeric epsilon = 0.01 * p_step; // This is the epsilon that
3862  // // we use for comparing p grid spacings.
3863 
3864  // Construct p_grid (new)
3865  // ----------------------
3866 
3867  ArrayOfNumeric log_p_new; // We take log_p_new as an array of Numeric, so
3868  // that we can easily build it up by appending new
3869  // elements to the end.
3870 
3871  // Check whether there are pressure levels that are further apart
3872  // (in log(p)) than p_step, and insert additional levels if
3873  // necessary:
3874 
3875  log_p_new.push_back(log_p_old[0]);
3876 
3877  for (Index i = 1; i < log_p_old.nelem(); ++i) {
3878  const Numeric dp =
3879  log_p_old[i - 1] - log_p_old[i]; // The grid is descending.
3880 
3881  const Numeric dp_by_p_step = dp / p_step;
3882  // cout << "dp_by_p_step: " << dp_by_p_step << "\n";
3883 
3884  // How many times does p_step fit into dp?
3885  const Index n = (Index)ceil(dp_by_p_step);
3886  // n is the number of intervals that we want to have in the
3887  // new grid. The number of additional points to insert is
3888  // n-1. But we have to insert the original point as well.
3889  // cout << n << "\n";
3890 
3891  const Numeric ddp = dp / (Numeric)n;
3892  // cout << "ddp: " << ddp << "\n";
3893 
3894  for (Index j = 1; j <= n; ++j)
3895  log_p_new.push_back(log_p_old[i - 1] - (Numeric)j * ddp);
3896  }
3897 
3898  // Copy ArrayOfNumeric to proper vector.
3899  Vector log_p(log_p_new.nelem());
3900  for (Index i = 0; i < log_p_new.nelem(); ++i) log_p[i] = log_p_new[i];
3901 
3902  // Copy the new grid to abs_p, removing the log:
3903  p_grid.resize(log_p.nelem());
3904  transform(p_grid, exp, log_p);
3905 }
3906 
3907 /* Workspace method: Doxygen documentation will be auto-generated */
3908 void p_gridFromZRaw( //WS Output
3909  Vector& p_grid,
3910  //WS Input
3911  const GriddedField3& z_field_raw,
3912  const Index& no_negZ,
3913  const Verbosity&) {
3914  // original version excludes negative z. not clear, why this is. maybe is
3915  // currently a convention somehwere in ARTS (DOIT?). negative z seem, however,
3916  // to work fine for clear-sky cases. so we make the negative z exclude an
3917  // option (for consistency until unclarities solved, default: do exclude)
3918  Vector p_grid_raw = z_field_raw.get_numeric_grid(GFIELD3_P_GRID);
3919 
3920  Index i;
3921  if (is_increasing(z_field_raw.data(joker, 0, 0))) {
3922  i = 0;
3923  if (no_negZ) {
3924  while (z_field_raw.data(i, 0, 0) < 0.0) i++;
3925  }
3926  p_grid = p_grid_raw[Range(i, joker)];
3927  } else if (is_decreasing(z_field_raw.data(joker, 0, 0))) {
3928  i = z_field_raw.data.npages() - 1;
3929  if (no_negZ) {
3930  while (z_field_raw.data(i, 0, 0) < 0.0) i--;
3931  }
3932  p_grid = p_grid_raw[Range(i, joker, -1)];
3933  } else {
3934  ostringstream os;
3935  os << "z_field_raw needs to be monotonous, but this is not the case.\n";
3936  throw runtime_error(os.str());
3937  }
3938 }
3939 
3940 /* Workspace method: Doxygen documentation will be auto-generated */
3941 void lat_gridFromZRaw( //WS Output
3942  Vector& lat_grid,
3943  //WS Input
3944  const GriddedField3& z_field_raw,
3945  const Verbosity&) {
3946  lat_grid = z_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
3947 }
3948 
3949 /* Workspace method: Doxygen documentation will be auto-generated */
3950 void lon_gridFromZRaw( //WS Output
3951  Vector& lon_grid,
3952  //WS Input
3953  const GriddedField3& z_field_raw,
3954  const Verbosity&) {
3955  lon_grid = z_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
3956 }
3957 
3958 /* Workspace method: Doxygen documentation will be auto-generated */
3959 void atm_gridsFromZRaw( //WS Output
3960  Vector& p_grid,
3961  Vector& lat_grid,
3962  Vector& lon_grid,
3963  //WS Input
3964  const GriddedField3& z_field_raw,
3965  const Index& no_negZ,
3966  const Verbosity& v) {
3967  p_gridFromZRaw(p_grid, z_field_raw, no_negZ, v);
3968  lat_gridFromZRaw(lat_grid, z_field_raw, v);
3969  lon_gridFromZRaw(lon_grid, z_field_raw, v);
3970 }
3971 
3972 /* Workspace method: Doxygen documentation will be auto-generated */
3973 void lat_gridFromRawField( //WS Output
3974  Vector& lat_grid,
3975  //WS Input
3976  const GriddedField3& field_raw,
3977  const Verbosity&) {
3978  lat_grid = field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
3979 }
3980 
3981 /* Workspace method: Doxygen documentation will be auto-generated */
3982 void lon_gridFromRawField( //WS Output
3983  Vector& lon_grid,
3984  //WS Input
3985  const GriddedField3& field_raw,
3986  const Verbosity&) {
3987  lon_grid = field_raw.get_numeric_grid(GFIELD3_LON_GRID);
3988 }
3989 
3990 /* Workspace method: Doxygen documentation will be auto-generated */
3992  const Index& atmosphere_dim,
3993  const Vector& p_grid,
3994  const Vector& lat_grid,
3995  const Vector& lon_grid,
3996  const Vector& refellipsoid,
3997  const Tensor3& z_field,
3998  const Numeric& planet_rotation_period,
3999  const Verbosity&) {
4000  if (atmosphere_dim < 3)
4001  throw runtime_error("No need to use this method for 1D and 2D.");
4002 
4003  const Index np = p_grid.nelem();
4004  const Index na = lat_grid.nelem();
4005  const Index no = lon_grid.nelem();
4006 
4007  chk_atm_field("z_field", z_field, atmosphere_dim, p_grid, lat_grid, lon_grid);
4008  if (wind_u_field.npages() > 0) {
4009  chk_atm_field("wind_u_field",
4010  wind_u_field,
4011  atmosphere_dim,
4012  p_grid,
4013  lat_grid,
4014  lon_grid);
4015  } else {
4016  wind_u_field.resize(np, na, no);
4017  wind_u_field = 0.;
4018  }
4019 
4020  const Numeric k1 = 2 * PI / planet_rotation_period;
4021 
4022  for (Index a = 0; a < na; a++) {
4023  const Numeric k2 = k1 * cos(DEG2RAD * lat_grid[a]);
4024  const Numeric re = refell2r(refellipsoid, lat_grid[a]);
4025 
4026  for (Index o = 0; o < no; o++) {
4027  for (Index p = 0; p < np; p++) {
4028  wind_u_field(p, a, o) += k2 * (re + z_field(p, a, o));
4029  }
4030  }
4031  }
4032 }
4033 
4034 // A small help function
4035 void z2g(Numeric& g, const Numeric& r, const Numeric& g0, const Numeric& z) {
4036  const Numeric x = r / (r + z);
4037  g = g0 * x * x;
4038 }
4039 
4040 /* Workspace method: Doxygen documentation will be auto-generated */
4042  Tensor3& z_field,
4043  const Index& atmosphere_dim,
4044  const Vector& p_grid,
4045  const Vector& lat_grid,
4046  const Vector& lon_grid,
4047  const Vector& lat_true,
4048  const Vector& lon_true,
4049  const ArrayOfArrayOfSpeciesTag& abs_species,
4050  const Tensor3& t_field,
4051  const Tensor4& vmr_field,
4052  const Vector& refellipsoid,
4053  const Matrix& z_surface,
4054  const Index& atmfields_checked,
4055  const Agenda& g0_agenda,
4056  const Numeric& molarmass_dry_air,
4057  const Numeric& p_hse,
4058  const Numeric& z_hse_accuracy,
4059  const Verbosity& verbosity) {
4060  if (atmfields_checked != 1)
4061  throw runtime_error(
4062  "The atmospheric fields must be flagged to have "
4063  "passed a consistency check (atmfields_checked=1).");
4064 
4065  // Some general variables
4066  const Index np = p_grid.nelem();
4067  const Index nlat = t_field.nrows();
4068  const Index nlon = t_field.ncols();
4069  //
4070  const Index firstH2O = find_first_species_tg(
4071  abs_species, species_index_from_species_name("H2O"));
4072 
4073  if (firstH2O < 0) {
4074  ostringstream os;
4075  os << "No water vapour tag group in *abs_species*.\n"
4076  << "Be aware that this leads to significant variations in atmospheres\n"
4077  << "that contain considerable amounts of water vapour (e.g. Earth)!\n";
4078  CREATE_OUT1;
4079  out1 << os.str();
4080  //throw runtime_error(
4081  // "Water vapour is required (must be a tag group in *abs_species*)." );
4082  }
4083  //
4084  if (p_hse > p_grid[0] || p_hse < p_grid[np - 1]) {
4085  ostringstream os;
4086  os << "The value of *p_hse* must be inside the range of *p_grid*:"
4087  << " p_hse = " << p_hse << " Pa\n"
4088  << " p_grid = << p_grid[np-1]"
4089  << " - " << p_grid[0] << " Pa\n";
4090  throw runtime_error(os.str());
4091  }
4092  //
4093  if (z_hse_accuracy <= 0) {
4094  throw runtime_error("The value of *z_hse_accuracy* must be > 0.");
4095  }
4096  //
4097  chk_latlon_true(atmosphere_dim, lat_grid, lat_true, lon_true);
4098 
4099  // Determine interpolation weights for p_hse
4100  //
4101  ArrayOfGridPos gp(1);
4102  Matrix itw(1, 2);
4103  p2gridpos(gp, p_grid, Vector(1, p_hse));
4104  interpweights(itw, gp);
4105 
4106  // // Molecular weight of water vapour
4107  const Numeric mw = 18.016;
4108 
4109  // mw/molarmass_dry_air matches eps in Eq. 3.14 in Wallace&Hobbs:
4110  const Numeric k = 1 - mw / molarmass_dry_air;
4111 
4112  // Gas constant for 1kg dry air:
4113  const Numeric rd = 1e3 * GAS_CONSTANT / molarmass_dry_air;
4114 
4115  // The calculations
4116  //
4117  for (Index ilat = 0; ilat < nlat; ilat++) {
4118  // The reference ellipsoid is already adjusted to internal 1D or 2D
4119  // views, and lat_grid is the relevant grid for *refellipsoid*, also
4120  // for 2D. On the other hand, extraction of g0 requires that the true
4121  // position is determined.
4122 
4123  // Radius of reference ellipsoid
4124  Numeric re;
4125  if (atmosphere_dim == 1) {
4126  re = refellipsoid[0];
4127  } else {
4128  re = refell2r(refellipsoid, lat_grid[ilat]);
4129  }
4130 
4131  for (Index ilon = 0; ilon < nlon; ilon++) {
4132  // Determine true latitude and longitude
4133  Numeric lat, lon;
4134  Vector pos(atmosphere_dim); // pos[0] can be a dummy value
4135  if (atmosphere_dim >= 2) {
4136  pos[1] = lat_grid[ilat];
4137  if (atmosphere_dim == 3) {
4138  pos[2] = lon_grid[ilon];
4139  }
4140  }
4142  lat, lon, atmosphere_dim, lat_grid, lat_true, lon_true, pos);
4143 
4144  // Get g0
4145  Numeric g0;
4146  g0_agendaExecute(ws, g0, lat, lon, g0_agenda);
4147 
4148  // Determine altitude for p_hse
4149  Vector z_hse(1);
4150  interp(z_hse, itw, z_field(joker, ilat, ilon), gp);
4151 
4152  Numeric z_acc = 2 * z_hse_accuracy;
4153 
4154  while (z_acc > z_hse_accuracy) {
4155  z_acc = 0;
4156  Numeric g2 = g0;
4157 
4158  for (Index ip = 0; ip < np - 1; ip++) {
4159  // Calculate average g
4160  if (ip == 0) {
4161  z2g(g2, re, g0, z_field(ip, ilat, ilon));
4162  }
4163  const Numeric g1 = g2;
4164  z2g(g2, re, g0, z_field(ip + 1, ilat, ilon));
4165  //
4166  const Numeric g = (g1 + g2) / 2.0;
4167 
4168  //Average water vapour VMR
4169  Numeric hm;
4170  if (firstH2O < 0) {
4171  hm = 0.0;
4172  } else {
4173  hm = 0.5 * (vmr_field(firstH2O, ip, ilat, ilon) +
4174  vmr_field(firstH2O, ip + 1, ilat, ilon));
4175  }
4176 
4177  // Average virtual temperature (no liquid water)
4178  // (cf. e.g. Eq. 3.16 in Wallace&Hobbs)
4179  const Numeric tv =
4180  (1 / (2 * (1 - hm * k))) *
4181  (t_field(ip, ilat, ilon) + t_field(ip + 1, ilat, ilon));
4182 
4183  // Change in vertical altitude from ip to ip+1
4184  // (cf. e.g. Eq. 3.24 in Wallace&Hobbs)
4185  const Numeric dz = rd * (tv / g) * log(p_grid[ip] / p_grid[ip + 1]);
4186 
4187  // New altitude
4188  Numeric znew = z_field(ip, ilat, ilon) + dz;
4189  z_acc = max(z_acc, fabs(znew - z_field(ip + 1, ilat, ilon)));
4190  z_field(ip + 1, ilat, ilon) = znew;
4191  }
4192 
4193  // Adjust to z_hse
4194  Vector z_tmp(1);
4195  interp(z_tmp, itw, z_field(joker, ilat, ilon), gp);
4196  z_field(joker, ilat, ilon) -= z_tmp[0] - z_hse[0];
4197  }
4198  }
4199  }
4200 
4201  // Check that there is no gap between the surface and lowest pressure
4202  // level
4203  // (This code is copied from *basics_checkedCalc*. Make this to an internal
4204  // function if used in more places.)
4205  for (Index row = 0; row < z_surface.nrows(); row++) {
4206  for (Index col = 0; col < z_surface.ncols(); col++) {
4207  if (z_surface(row, col) < z_field(0, row, col) ||
4208  z_surface(row, col) >= z_field(z_field.npages() - 1, row, col)) {
4209  ostringstream os;
4210  os << "The surface altitude (*z_surface*) cannot be outside "
4211  << "of the altitudes in *z_field*.";
4212  throw runtime_error(os.str());
4213  }
4214  }
4215  }
4216 }
4217 
4218 /* Workspace method: Doxygen documentation will be auto-generated */
4220  const ArrayOfArrayOfSpeciesTag& abs_species,
4221  const String& species,
4222  const Numeric& vmr_value,
4223  const Verbosity&) {
4224  // Check input
4225  chk_if_in_range("vmr_value", vmr_value, 0, 1);
4226  //
4227  if (abs_species.nelem() != vmr_field.nbooks())
4228  throw runtime_error(
4229  "Size of *vmr_field* and length of *abs_species* do not agree.");
4230 
4231  // Find position for this species.
4232  ArrayOfSpeciesTag tag;
4233  array_species_tag_from_string(tag, species);
4234  Index si = chk_contains("species", abs_species, tag);
4235 
4236  // Apply value
4237  vmr_field(si, joker, joker, joker) = vmr_value;
4238 }
4239 
4240 /* Workspace method: Doxygen documentation will be auto-generated */
4242  const ArrayOfArrayOfSpeciesTag& abs_species,
4243  const Vector& vmr_values,
4244  const Verbosity& verbosity) {
4245  CREATE_OUT3;
4246 
4247  const Index nspecies = abs_species.nelem();
4248 
4249  if (vmr_values.nelem() not_eq nspecies)
4250  throw std::runtime_error("Not same number of vmr_values as abs_species.");
4251 
4252  out3 << "Setting all " << nspecies << " species to constant VMR\n";
4253 
4254  for (Index i = 0; i < nspecies; i++) {
4255  const ArrayOfSpeciesTag& a_abs_species = abs_species[i];
4256  const String species_tag_name = get_tag_group_name(a_abs_species);
4258  vmr_field, abs_species, species_tag_name, vmr_values[i], verbosity);
4259  }
4260 }
4261 
4262 /* Workspace method: Doxygen documentation will be auto-generated */
4264  Index& nlte_do,
4265  EnergyLevelMap& nlte_field,
4266  ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
4267  const ArrayOfQuantumIdentifier& nlte_quantum_identifiers,
4268  const SpeciesAuxData& partition_functions,
4269  const Tensor3& t_field,
4270  const Verbosity& verbosity) {
4271  using Constant::h;
4272 
4273  CREATE_OUT2;
4274  const Index nn = nlte_quantum_identifiers.nelem(), np = t_field.npages(),
4275  nlat = t_field.nrows(), nlon = t_field.ncols();
4276  if (nn == 0) return;
4277 
4278  Tensor4 nlte_tensor4(nn, np, nlat, nlon);
4279  nlte_do = 1;
4280  ArrayOfIndex checked(nn, 0);
4281 
4282  for (Index in = 0; in < nn; in++) {
4283  const QuantumIdentifier& qi = nlte_quantum_identifiers[in];
4284  Tensor3View lte = nlte_tensor4(in, joker, joker, joker);
4285 
4286  for (auto& abs_lines : abs_lines_per_species) {
4287  for (auto& band : abs_lines) {
4288  for (Index k=0; k<band.NumLines(); k++) {
4289  if (Absorption::id_in_line_lower(band, qi, k)) {
4291 
4292  if (not checked[in]) {
4293  checked[in] = 1;
4294 
4295  for (Index ip = 0; ip < np; ip++) {
4296  for (Index ilat = 0; ilat < nlat; ilat++) {
4297  for (Index ilon = 0; ilon < nlon; ilon++) {
4298  lte(ip, ilat, ilon) =
4299  boltzman_factor(t_field(ip, ilat, ilon), band.E0(k)) *
4300  band.g_low(k) / single_partition_function(
4301  t_field(ip, ilat, ilon), partition_functions.getParamType(band.Species(),
4302  band.Isotopologue()),
4303  partition_functions.getParam(band.Species(),
4304  band.Isotopologue()));
4305  }
4306  }
4307  }
4308  }
4309  }
4310 
4311  if (Absorption::id_in_line_upper(band, qi, k)) {
4313 
4314  if (not checked[in]) {
4315  checked[in] = 1;
4316  for (Index ip = 0; ip < np; ip++) {
4317  for (Index ilat = 0; ilat < nlat; ilat++) {
4318  for (Index ilon = 0; ilon < nlon; ilon++) {
4319  lte(ip, ilat, ilon) =
4320  boltzman_factor(t_field(ip, ilat, ilon), band.E0(k) + h*band.F0(k)) *
4321  band.g_upp(k) / single_partition_function(
4322  t_field(ip, ilat, ilon), partition_functions.getParamType(band.Species(),
4323  band.Isotopologue()),
4324  partition_functions.getParam(band.Species(),
4325  band.Isotopologue()));
4326  }
4327  }
4328  }
4329  }
4330  }
4331  }
4332  }
4333  }
4334  }
4335 
4336  for (Index in = 0; in < nn; in++) {
4337  if (not checked[in]) {
4338  out2 << "Did not find match among lines for: "
4339  << nlte_quantum_identifiers[in] << "\n";
4340  }
4341  }
4342 
4343  nlte_field = EnergyLevelMap(nlte_tensor4, nlte_quantum_identifiers);
4344 }
4345 
4346 /* Workspace method: Doxygen documentation will be auto-generated */
4348  Index& nlte_do,
4349  EnergyLevelMap& nlte_field,
4350  ArrayOfArrayOfAbsorptionLines& abs_lines_per_species,
4351  const ArrayOfQuantumIdentifier& nlte_quantum_identifiers,
4352  const Tensor3& t_field,
4353  const Verbosity& verbosity) {
4354  using Constant::h;
4355 
4356  CREATE_OUT2;
4357  const Index nn = nlte_quantum_identifiers.nelem(), np = t_field.npages(),
4358  nlat = t_field.nrows(), nlon = t_field.ncols();
4359  if (nn == 0) return;
4360 
4361  // Find where they are positioned and how many different molecules there are for the NLTE fields
4362  ArrayOfIndex part_fun_pos(nn, 0);
4363  Index x = 1;
4364  for (Index in = 1; in < nn; in++) {
4365  bool found = false;
4366  for (Index ix = 0; ix < in; ix++) {
4367  if (nlte_quantum_identifiers[in].Species() ==
4368  nlte_quantum_identifiers[ix].Species() and
4369  nlte_quantum_identifiers[in].Isotopologue() ==
4370  nlte_quantum_identifiers[ix].Isotopologue()) {
4371  part_fun_pos[in] = part_fun_pos[ix];
4372  found = true;
4373  break;
4374  }
4375  }
4376  if (not found) {
4377  part_fun_pos[in] = x;
4378  x++;
4379  }
4380  }
4381 
4382  Tensor4 part_fun(x, np, nlat, nlon, 0.0);
4383  Tensor4 nlte_tensor4(nn, np, nlat, nlon, 0);
4384  nlte_do = 1;
4385  ArrayOfIndex checked(nn, 0);
4386 
4387  for (Index in = 0; in < nn; in++) {
4388  const QuantumIdentifier& qi = nlte_quantum_identifiers[in];
4389  Tensor3View lte = nlte_tensor4(in, joker, joker, joker);
4390 
4391  for (auto& abs_lines : abs_lines_per_species) {
4392  for (auto& band : abs_lines) {
4393  for (Index k=0; k<band.NumLines(); k++) {
4394  if (Absorption::id_in_line_lower(band, qi, k)) {
4396 
4397  if (not checked[in]) {
4398  checked[in] = 1;
4399 
4400  for (Index ip = 0; ip < np; ip++) {
4401  for (Index ilat = 0; ilat < nlat; ilat++) {
4402  for (Index ilon = 0; ilon < nlon; ilon++) {
4403  lte(ip, ilat, ilon) =
4404  boltzman_factor(t_field(ip, ilat, ilon), band.E0(k)) *
4405  band.g_low(k);
4406  part_fun(part_fun_pos[in], ip, ilat, ilon) +=
4407  lte(ip, ilat, ilon);
4408  }
4409  }
4410  }
4411  }
4412  }
4413 
4414  if (Absorption::id_in_line_upper(band, qi, k)) {
4416 
4417  if (not checked[in]) {
4418  checked[in] = 1;
4419 
4420  for (Index ip = 0; ip < np; ip++) {
4421  for (Index ilat = 0; ilat < nlat; ilat++) {
4422  for (Index ilon = 0; ilon < nlon; ilon++) {
4423  lte(ip, ilat, ilon) =
4424  boltzman_factor(t_field(ip, ilat, ilon),
4425  band.E0(k) + h * band.F0(k)) *
4426  band.g_upp(k);
4427  part_fun(part_fun_pos[in], ip, ilat, ilon) +=
4428  lte(ip, ilat, ilon);
4429  }
4430  }
4431  }
4432  }
4433  }
4434  }
4435  }
4436  }
4437  }
4438 
4439  for (Index in = 0; in < nn; in++) {
4440  if (not checked[in]) {
4441  out2 << "Did not find match among lines for: "
4442  << nlte_quantum_identifiers[in] << "\n";
4443  }
4444  }
4445 
4446  for (Index in = 0; in < nn; in++) {
4447  if (checked[in]) {
4448  nlte_tensor4(in, joker, joker, joker) /=
4449  part_fun(part_fun_pos[in], joker, joker, joker);
4450  }
4451  }
4452 
4453  nlte_field = EnergyLevelMap(nlte_tensor4, nlte_quantum_identifiers);
4454 }
const Index GFIELD3_LON_GRID
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void batch_atm_fields_compactCleanup(ArrayOfGriddedField4 &batch_atm_fields_compact, const Numeric &threshold, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactCleanup.
void vmr_fieldSetConstant(Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Numeric &vmr_value, const Verbosity &)
WORKSPACE METHOD: vmr_fieldSetConstant.
bool is_lon_cyclic(ConstVectorView grid, const Numeric &epsilon)
Check if the given longitude grid is cyclic.
Definition: logic.cc:369
void AtmFieldPRegridHelper(Index &ing_min, Index &ing_max, ArrayOfGridPosPoly &gp_p, Matrix &itw, ConstVectorView p_grid_out, ConstVectorView p_grid_in, const Index &interp_order, const Verbosity &verbosity)
Calculate grid positions and interpolations weights for AtmFieldPRegrid.
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const QuantumIdentifier &self, const ArrayOfSpeciesTag &lineshape_species, bool self_in_list, bool bath_in_list, Type type)
Returns a VMR vector for this model&#39;s main calculations.
void lon_gridFromZRaw(Vector &lon_grid, const GriddedField3 &z_field_raw, const Verbosity &)
WORKSPACE METHOD: lon_gridFromZRaw.
void z_surfaceConstantAltitude(Matrix &z_surface, const Vector &lat_grid, const Vector &lon_grid, const Numeric &altitude, const Verbosity &verbosity)
WORKSPACE METHOD: z_surfaceConstantAltitude.
void GriddedFieldLatLonRegrid(GriddedField2 &gfraw_out, const Vector &lat_true, const Vector &lon_true, const GriddedField2 &gfraw_in_orig, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: GriddedFieldLatLonRegrid.
const Index GFIELD3_LAT_GRID
The Tensor4View class.
Definition: matpackIV.h:284
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:195
void batch_atm_fields_compactAddConstant(ArrayOfGriddedField4 &batch_atm_fields_compact, const String &name, const Numeric &value, const Index &prepend, const ArrayOfString &condensibles, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactAddConstant.
QuantumIdentifier::QType Index LowerQuantumNumbers Species
void GriddedFieldZToPRegridHelper(Index &ing_min, Index &ing_max, ArrayOfGridPosPoly &gp_p, Matrix &itw, const GriddedField &gfraw_in, const Index z_grid_index, ConstVectorView z_grid, const Index &interp_order, const Index &zeropadding, const Verbosity &verbosity)
Calculate grid positions and interpolations weights for GriddedFieldZToPRegrid.
void atm_gridsFromZRaw(Vector &p_grid, Vector &lat_grid, Vector &lon_grid, const GriddedField3 &z_field_raw, const Index &no_negZ, const Verbosity &v)
WORKSPACE METHOD: atm_gridsFromZRaw.
Declarations having to do with the four output streams.
void WindRawRead(GriddedField3 &wind_u_field_raw, GriddedField3 &wind_v_field_raw, GriddedField3 &wind_w_field_raw, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: WindRawRead.
void interp_atmfield_by_gp(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Interpolates an atmospheric field given the grid positions.
void parse_atmcompact_speciestype(String &species_type, const String &field_name, const String &delim)
Definition: cloudbox.cc:848
bool checksize() const final
Consistency check.
The Vector class.
Definition: matpackI.h:860
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
const ArrayOfString & get_string_grid(Index i) const
Get a string grid.
const Tensor4 & Data() const noexcept
Energy level type.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
The Tensor4 class.
Definition: matpackIV.h:421
The range class.
Definition: matpackI.h:160
void AtmRawRead(GriddedField3 &t_field_raw, GriddedField3 &z_field_raw, ArrayOfGriddedField3 &vmr_field_raw, ArrayOfGriddedField3 &nlte_field_raw, ArrayOfQuantumIdentifier &nlte_quantum_identifiers, Vector &nlte_vibrational_energies, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: AtmRawRead.
String get_tag_group_name(const ArrayOfSpeciesTag &tg)
Return the name of a tag group as a string.
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:60
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:49
void nlte_fieldSetLteExternalPartitionFunction(Index &nlte_do, EnergyLevelMap &nlte_field, ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const ArrayOfQuantumIdentifier &nlte_quantum_identifiers, const SpeciesAuxData &partition_functions, const Tensor3 &t_field, const Verbosity &verbosity)
WORKSPACE METHOD: nlte_fieldSetLteExternalPartitionFunction.
const Numeric GAS_CONSTANT
Numeric boltzman_factor(Numeric T, Numeric E0)
Computes exp(- E0/kT)
Definition: linescaling.cc:176
const Vector & Energies() const noexcept
Energy level type.
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:680
void nlte_fieldSetLteInternalPartitionFunction(Index &nlte_do, EnergyLevelMap &nlte_field, ArrayOfArrayOfAbsorptionLines &abs_lines_per_species, const ArrayOfQuantumIdentifier &nlte_quantum_identifiers, const Tensor3 &t_field, const Verbosity &verbosity)
WORKSPACE METHOD: nlte_fieldSetLteInternalPartitionFunction.
This file contains basic functions to handle XML data files.
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
const Index GFIELD3_P_GRID
void p_gridRefine(Vector &p_grid, Index &atmfields_checked, Index &atmgeom_checked, Index &cloudbox_checked, const Vector &p_grid_old, const Numeric &p_step10, const Verbosity &)
WORKSPACE METHOD: p_gridRefine.
Header file for interpolation.cc.
#define min(a, b)
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:63
bool checksize() const final
Consistency check.
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
const Array< SpeciesRecord > species_data
Species Data.
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:989
Structure to store a grid position.
Definition: interpolation.h:73
void gridpos_poly_longitudinal(const String &error_msg, ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation on longitudes.
const Numeric PI
void MagRawRead(GriddedField3 &mag_u_field_raw, GriddedField3 &mag_v_field_raw, GriddedField3 &mag_w_field_raw, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: MagRawRead.
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
const AuxType & getParamType(const Index species, const Index isotopologue) const
Return a constant reference to the parameter types.
Definition: absorption.h:276
The Tensor3 class.
Definition: matpackIII.h:339
void GriddedFieldLatLonRegridHelper(ArrayOfGridPosPoly &gp_lat, ArrayOfGridPosPoly &gp_lon, Tensor3 &itw, GriddedField &gfraw_out, const GriddedField &gfraw_in, const Index lat_grid_index, const Index lon_grid_index, ConstVectorView lat_true, ConstVectorView lon_true, const Index &interp_order, const Verbosity &verbosity)
Calculate grid positions and interpolations weights for GriddedFieldLatLonRegrid. ...
The global header file for ARTS.
void p_gridDensify(Vector &p_grid, Index &atmfields_checked, Index &atmgeom_checked, Index &cloudbox_checked, const Vector &p_grid_old, const Index &nfill, const Verbosity &verbosity)
WORKSPACE METHOD: p_gridDensify.
void AtmosphereSet1D(Index &atmosphere_dim, Vector &lat_grid, Vector &lon_grid, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet1D.
EnergyLevelMapType Type() const noexcept
Energy level type.
void GriddedFieldZToPRegrid(GriddedField3 &gfraw_out, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const GriddedField3 &gfraw_in_orig, const Index &interp_order, const Index &zeropadding, const Verbosity &verbosity)
WORKSPACE METHOD: GriddedFieldZToPRegrid.
void chk_interpolation_grids_loose_no_data_check(Index &ing_min, Index &ing_max, const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order)
Check interpolation grids.
Definition: check_input.cc:773
const Index GFIELD4_FIELD_NAMES
void AtmWithNLTERawRead(GriddedField3 &t_field_raw, GriddedField3 &z_field_raw, ArrayOfGriddedField3 &vmr_field_raw, ArrayOfGriddedField3 &nlte_field_raw, ArrayOfQuantumIdentifier &nlte_quantum_identifiers, Vector &nlte_vibrational_energies, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Index &expect_vibrational_energies, const Verbosity &verbosity)
WORKSPACE METHOD: AtmWithNLTERawRead.
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true, ConstVectorView pos)
Determines the true alt and lon for an "ARTS position".
Definition: rte.cc:2314
void p2gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Index order, const Numeric &extpolfac)
p2gridpos_poly
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:351
_CS_string_type str() const
Definition: sstream.h:491
void ThrowIfNotOK() const
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
void atm_fields_compactFromMatrix(GriddedField4 &af, const Index &atmosphere_dim, const Matrix &im, const ArrayOfString &field_names, const Verbosity &)
WORKSPACE METHOD: atm_fields_compactFromMatrix.
void chk_interpolation_pgrids(const String &which_interpolation, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Index order, const Numeric &extpolfac)
Check log pressure interpolation grids.
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:149
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
Declarations for agendas.
The Tensor3View class.
Definition: matpackIII.h:239
void set_grid(Index i, const Vector &g)
Set a numeric grid.
void toupper()
Convert to upper case.
Definition: mystring.h:74
void chk_griddedfield_gridname(const GriddedField &gf, const Index gridindex, const String &gridname)
Check name of grid in GriddedField.
void AtmFieldsCalc(Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, EnergyLevelMap &nlte_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &t_field_raw, const GriddedField3 &z_field_raw, const ArrayOfGriddedField3 &vmr_field_raw, const ArrayOfGriddedField3 &nlte_field_raw, const ArrayOfQuantumIdentifier &nlte_ids, const Vector &nlte_energies, const Index &atmosphere_dim, const Index &interp_order, const Index &vmr_zeropadding, const Index &vmr_nonegative, const Index &nlte_when_negative, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsCalc.
const ArrayOfQuantumIdentifier & Levels() const noexcept
Energy level type.
bool id_in_line_upper(const Lines &band, const QuantumIdentifier &id, size_t line_index)
Checks if the external quantum identifier match a line&#39;s ID.
void atm_fields_compactCreateFromField(GriddedField4 &atm_fields_compact, const String &name, const GriddedField3 &field, const Verbosity &)
WORKSPACE METHOD: atm_fields_compactCreateFromField.
Index get_grid_size(Index i) const
Get the size of a grid.
void AtmosphereSet2D(Index &atmosphere_dim, Vector &lon_grid, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet2D.
QuantumIdentifier::QType Isotopologue
#define CREATE_OUT1
Definition: messages.h:205
void resize(const GriddedField2 &gf)
Make this GriddedField2 the same size as the given one.
void FieldFromGriddedField(Matrix &field_out, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField2 &gfraw_in, const Verbosity &)
WORKSPACE METHOD: FieldFromGriddedField.
const Index GFIELD4_LON_GRID
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
void AtmFieldsExtract1D(Index &atmosphere_dim, Vector &lat_grid, Vector &lon_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, const Index &ilat, const Index &ilon, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsExtract1D.
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void AtmFieldsCalcExpand1D(Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, EnergyLevelMap &nlte_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &t_field_raw, const GriddedField3 &z_field_raw, const ArrayOfGriddedField3 &vmr_field_raw, const ArrayOfGriddedField3 &nlte_field_raw, const ArrayOfQuantumIdentifier &nlte_ids, const Vector &nlte_energies, const Index &atmosphere_dim, const Index &interp_order, const Index &vmr_zeropadding, const Index &vmr_nonegative, const Index &nlte_when_negative, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsCalcExpand1D.
Class to identify and match lines by their quantum numbers.
Definition: quantum.h:390
const Vector & get_numeric_grid(Index i) const
Get a numeric grid.
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
const Joker joker
const ArrayOfGriddedField1 & getParam(const Index species, const Index isotopologue) const
Return a constant reference to the parameters.
Definition: absorption.cc:145
void GriddedFieldPRegrid(GriddedField3 &gfraw_out, const Vector &p_grid, const GriddedField3 &gfraw_in_orig, const Index &interp_order, const Index &zeropadding, const Verbosity &verbosity)
WORKSPACE METHOD: GriddedFieldPRegrid.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
void VectorNLogSpace(Vector &x, const Index &n, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
void AtmFieldsAndParticleBulkPropFieldFromCompact(Vector &p_grid, Vector &lat_grid, Vector &lon_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, Tensor4 &particle_bulkprop_field, ArrayOfString &particle_bulkprop_names, const ArrayOfArrayOfSpeciesTag &abs_species, const GriddedField4 &atm_fields_compact, const Index &atmosphere_dim, const String &delim, const Numeric &p_min, const Index &check_gridnames, const Verbosity &)
WORKSPACE METHOD: AtmFieldsAndParticleBulkPropFieldFromCompact.
const Numeric EPSILON_LON_CYCLIC
Data value accuracy requirement for values at 0 and 360 deg if longitudes are cyclic.
Declarations required for the calculation of absorption coefficients.
void MagFieldsFromAltitudeRawCalc(Tensor3 &mag_u_field, Tensor3 &mag_v_field, Tensor3 &mag_w_field, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const GriddedField3 &mag_u_field_raw, const GriddedField3 &mag_v_field_raw, const GriddedField3 &mag_w_field_raw, const Index &interp_order, const Numeric &extrapolation_factor, const Verbosity &verbosity)
WORKSPACE METHOD: MagFieldsFromAltitudeRawCalc.
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1135
const Index GFIELD4_LAT_GRID
void FieldFromGriddedFieldCheckLatLonHelper(const Vector &lat_grid, const Vector &lon_grid, const Index ilat, const Index ilon, const GriddedField &gfield)
Check for correct grid dimensions.
Header file for special_interp.cc.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:664
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:531
void p_gridFromZRaw(Vector &p_grid, const GriddedField3 &z_field_raw, const Index &no_negZ, const Verbosity &)
WORKSPACE METHOD: p_gridFromZRaw.
void lat_gridFromZRaw(Vector &lat_grid, const GriddedField3 &z_field_raw, const Verbosity &)
WORKSPACE METHOD: lat_gridFromZRaw.
void chk_interpolation_pgrids_loose_no_data_check(Index &ing_min, Index &ing_max, const String &which_interpolation, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Index order)
Check log pressure interpolation grids.
Definition: check_input.cc:890
Array< QuantumIdentifier > ArrayOfQuantumIdentifier
Definition: quantum.h:747
void WindFieldsCalcExpand1D(Tensor3 &wind_u_field, Tensor3 &wind_v_field, Tensor3 &wind_w_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &wind_u_field_raw, const GriddedField3 &wind_v_field_raw, const GriddedField3 &wind_w_field_raw, const Index &atmosphere_dim, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: WindFieldsCalcExpand1D.
basic_ostringstream< char, string_char_traits< char >, alloc > ostringstream
Definition: sstream.h:204
void batch_atm_fields_compactFromArrayOfMatrix(ArrayOfGriddedField4 &batch_atm_fields_compact, const Index &atmosphere_dim, const ArrayOfMatrix &am, const ArrayOfString &field_names, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactFromArrayOfMatrix.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
This can be used to make arrays out of anything.
Definition: array.h:40
void WindFieldsCalc(Tensor3 &wind_u_field, Tensor3 &wind_v_field, Tensor3 &wind_w_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &wind_u_field_raw, const GriddedField3 &wind_v_field_raw, const GriddedField3 &wind_w_field_raw, const Index &atmosphere_dim, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: WindFieldsCalc.
void AtmFieldsExpand1D(Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Index &atmosphere_dim, const Index &chk_vmr_nan, const Verbosity &)
WORKSPACE METHOD: AtmFieldsExpand1D.
Header file for interpolation_poly.cc.
const String & get_grid_name(Index i) const
Get grid name.
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:466
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
const Index GFIELD4_P_GRID
void batch_atm_fields_compactAddSpecies(ArrayOfGriddedField4 &batch_atm_fields_compact, const String &name, const GriddedField3 &species, const Index &prepend, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactAddSpecies.
Constains various line scaling functions.
void InterpAtmFieldToPosition(Numeric &outvalue, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &rtp_pos, const Tensor3 &field, const Verbosity &verbosity)
WORKSPACE METHOD: InterpAtmFieldToPosition.
void atm_fields_compactAddSpecies(GriddedField4 &atm_fields_compact, const String &name, const GriddedField3 &species, const Index &prepend, const Verbosity &verbosity)
WORKSPACE METHOD: atm_fields_compactAddSpecies.
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
#define max(a, b)
void lon_gridFromRawField(Vector &lon_grid, const GriddedField3 &field_raw, const Verbosity &)
WORKSPACE METHOD: lon_gridFromRawField.
void chk_atm_field(const String &x_name, ConstTensor3View x, const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, const bool &chk_lat90)
chk_atm_field (simple fields)
const Numeric DEG2RAD
A constant view of a Vector.
Definition: matpackI.h:476
Header file for stuff related to absorption species tags.
void g0_agendaExecute(Workspace &ws, Numeric &g0, const Numeric lat, const Numeric lon, const Agenda &input_agenda)
Definition: auto_md.cc:23858
Index nelem(const Lines &l)
Number of lines.
void chk_if_equal(const String &x1_name, const String &x2_name, ConstVectorView v1, ConstVectorView v2, Numeric margin)
chk_if_equal
Definition: check_input.cc:381
void GriddedFieldPRegridHelper(Index &ing_min, Index &ing_max, ArrayOfGridPosPoly &gp_p, Matrix &itw, GriddedField &gfraw_out, const GriddedField &gfraw_in, const Index p_grid_index, ConstVectorView p_grid, const Index &interp_order, const Index &zeropadding, const Verbosity &verbosity)
Calculate grid positions and interpolations weights for GriddedFieldPRegrid.
void wind_u_fieldIncludePlanetRotation(Tensor3 &wind_u_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Tensor3 &z_field, const Numeric &planet_rotation_period, const Verbosity &)
WORKSPACE METHOD: wind_u_fieldIncludePlanetRotation.
void GriddedFieldLatLonExpand(GriddedField2 &gfraw_out, const GriddedField2 &gfraw_in_orig, const Verbosity &)
WORKSPACE METHOD: GriddedFieldLatLonExpand.
Workspace class.
Definition: workspace_ng.h:40
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
void set_grid_name(Index i, const String &s)
Set grid name.
void atm_fields_compactAddConstant(GriddedField4 &af, const String &name, const Numeric &value, const Index &prepend, const ArrayOfString &condensibles, const Verbosity &verbosity)
WORKSPACE METHOD: atm_fields_compactAddConstant.
void lat_gridFromRawField(Vector &lat_grid, const GriddedField3 &field_raw, const Verbosity &)
WORKSPACE METHOD: lat_gridFromRawField.
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1476
#define _U_
Definition: config.h:183
Implementation of gridded fields.
void parse_atmcompact_scattype(String &scat_type, const String &field_name, const String &delim)
Definition: cloudbox.cc:912
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
#define CREATE_OUT3
Definition: messages.h:207
Numeric single_partition_function(const Numeric &T, const SpeciesAuxData::AuxType &partition_type, const ArrayOfGriddedField1 &partition_data)
Computes the partition function at one temperature.
Definition: linescaling.cc:72
void z_fieldFromHSE(Workspace &ws, Tensor3 &z_field, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Index &atmfields_checked, const Agenda &g0_agenda, const Numeric &molarmass_dry_air, const Numeric &p_hse, const Numeric &z_hse_accuracy, const Verbosity &verbosity)
WORKSPACE METHOD: z_fieldFromHSE.
Auxiliary data for isotopologues.
Definition: absorption.h:217
Internal cloudbox functions.
void MagFieldsCalcExpand1D(Tensor3 &mag_u_field, Tensor3 &mag_v_field, Tensor3 &mag_w_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &mag_u_field_raw, const GriddedField3 &mag_v_field_raw, const GriddedField3 &mag_w_field_raw, const Index &atmosphere_dim, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MagFieldsCalcExpand1D.
void AtmFieldsRefinePgrid(Vector &p_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, Index &atmfields_checked, Index &atmgeom_checked, Index &cloudbox_checked, const Vector &lat_grid, const Vector &lon_grid, const Index &atmosphere_dim, const Numeric &p_step, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsRefinePgrid.
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:66
void atm_fields_compactCleanup(GriddedField4 &atm_fields_compact, const Numeric &threshold, const Verbosity &)
WORKSPACE METHOD: atm_fields_compactCleanup.
void resize(const GriddedField3 &gf)
Make this GriddedField3 the same size as the given one.
void AtmosphereSet3D(Index &atmosphere_dim, Vector &lat_true, Vector &lon_true, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet3D.
void vmr_fieldSetAllConstant(Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Vector &vmr_values, const Verbosity &verbosity)
WORKSPACE METHOD: vmr_fieldSetAllConstant.
void AtmFieldPRegrid(Tensor3 &atmtensor_out, const Tensor3 &atmtensor_in_orig, const Vector &p_grid_new, const Vector &p_grid_old, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldPRegrid.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void z_surfaceFromFileAndGrid(Matrix &z_surface, const Vector &lat_grid, const Vector &lon_grid, const String &filename, const Index &interp_order, const Index &set_lowest_altitude_to_zero, const Verbosity &verbosity)
WORKSPACE METHOD: z_surfaceFromFileAndGrid.
void parse_atmcompact_speciesname(String &species_name, const String &field_name, const String &delim)
Definition: cloudbox.cc:880
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
#define CREATE_OUT2
Definition: messages.h:206
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
void MagFieldsCalc(Tensor3 &mag_u_field, Tensor3 &mag_v_field, Tensor3 &mag_w_field, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField3 &mag_u_field_raw, const GriddedField3 &mag_v_field_raw, const GriddedField3 &mag_w_field_raw, const Index &atmosphere_dim, const Index &interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MagFieldsCalc.
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
void z2g(Numeric &g, const Numeric &r, const Numeric &g0, const Numeric &z)
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1064
void atm_fields_compactExpand(GriddedField4 &af, Index &nf, const String &name, const Index &prepend, const Verbosity &)
atm_fields_compactExpand
Definition: m_atmosphere.cc:99
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
bool id_in_line_lower(const Lines &band, const QuantumIdentifier &id, size_t line_index)
Checks if the external quantum identifier match a line&#39;s ID.
Array< GridPosPoly > ArrayOfGridPosPoly
An Array of grid positions.