ARTS  2.2.66
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 
22 /*===========================================================================
23  === File description
24  ===========================================================================*/
25 
40 /*===========================================================================
41  === External declarations
42  ===========================================================================*/
43 
44 #include <cmath>
45 #include <cfloat>
46 #include "agenda_class.h"
47 #include "arts.h"
48 #include "auto_md.h"
49 #include "absorption.h"
50 #include "abs_species_tags.h"
51 #include "check_input.h"
52 #include "cloudbox.h"
53 #include "geodetic.h"
54 #include "matpackIII.h"
55 #include "messages.h"
56 #include "gridded_fields.h"
57 #include "interpolation.h"
58 #include "interpolation_poly.h"
59 #include "rte.h"
60 #include "special_interp.h"
61 #include "xml_io.h"
62 #include "global_data.h"
63 
64 extern const Index GFIELD3_P_GRID;
65 extern const Index GFIELD3_LAT_GRID;
66 extern const Index GFIELD3_LON_GRID;
67 extern const Index GFIELD4_FIELD_NAMES;
68 extern const Index GFIELD4_P_GRID;
69 extern const Index GFIELD4_LAT_GRID;
70 extern const Index GFIELD4_LON_GRID;
71 
72 extern const Numeric DEG2RAD;
73 extern const Numeric PI;
74 extern const Numeric GAS_CONSTANT;
75 
77 
79 extern const Numeric EPSILON_LON_CYCLIC = 2*DBL_EPSILON;
80 
81 
82 /*===========================================================================
83  *=== Helper functions
84  *===========================================================================*/
85 
87 
102 void atm_fields_compactExpand(GriddedField4& af,
103  Index& nf,
104  const String& name,
105  const Verbosity&)
106 {
107  // Number of fields already present:
108  nf = af.get_string_grid(GFIELD4_FIELD_NAMES).nelem();
109 
110 // Most of the functionality moved from atm_fields_compactAddConstant when
111 // atm_fields_compactAddSpecies was added, in order to share the code.
112 //
113 // Commented out by Gerrit 2011-05-04. I believe this check is not needed.
114 // Oliver agrees. We can still infer the dimensions even if there are zero
115 // fields; e.g., the data might have dimensions (0, 4, 3, 8).
116 //
117 // if (0==nf)
118 // {
119 // ostringstream os;
120 // os << "The *atm_fields_compact* must already contain at least one field,\n"
121 // << "so that we can infer the dimensions from that.";
122 // throw runtime_error( os.str() );
123 // }
124 
125  // Add name of new field to field name list:
126  af.get_string_grid(GFIELD4_FIELD_NAMES).push_back(name);
127 
128  // Save original fields:
129  const Tensor4 dummy = af.data;
130 
131  // Adjust size:
132  af.resize( nf+1, dummy.npages(), dummy.nrows(), dummy.ncols() );
133  nf++; // set to new number of fields
134 
135  // Copy back original field:
136  af.data( Range(0,nf-1), Range(joker), Range(joker), Range(joker) ) = dummy;
137 }
138 
139 
140 
141 /*===========================================================================
142  === The functions (in alphabetical order)
143  ===========================================================================*/
144 
146 /*
147  Helper function for FieldFromGriddedField functions to ensure correct
148  dimension and values of latitude/longitude grids
149 
150  \param[in] lat_grid Latitude grid
151  \param[in] lon_grid Longitude grid
152  \param[in] ilat Latitude grid index in gfield
153  \param[in] ilon Longitude grid index in gfield
154  \param[in] gfield GriddedField
155  */
156 void FieldFromGriddedFieldCheckLatLonHelper(const Vector& lat_grid,
157  const Vector& lon_grid,
158  const Index ilat,
159  const Index ilon,
160  const GriddedField& gfield)
161 {
162  chk_griddedfield_gridname(gfield, ilat, "Latitude");
163  chk_griddedfield_gridname(gfield, ilon, "Longitude");
164 
165  if (lon_grid.nelem() == 0)
166  {
167  chk_size("gfield.lon_grid", gfield.get_numeric_grid(ilon), 1);
168 
169  if (lat_grid.nelem() == 0)
170  chk_size("gfield.lat_grid", gfield.get_numeric_grid(ilat), 1);
171  else
172  chk_if_equal("lat_grid", "gfield.lat_grid", lat_grid, gfield.get_numeric_grid(ilat));
173  }
174  else
175  {
176  chk_if_equal("lat_grid", "gfield.lat_grid", lat_grid,
177  gfield.get_numeric_grid(ilat));
178  chk_if_equal("lon_grid", "gfield.lon_grid", lon_grid,
179  gfield.get_numeric_grid(ilon));
180  }
181 }
182 
183 /* Workspace method: Doxygen documentation will be auto-generated */
184 void FieldFromGriddedField(// WS Generic Output:
185  Matrix& field_out,
186  // WS Input:
187  const Vector& p_grid _U_,
188  const Vector& lat_grid,
189  const Vector& lon_grid,
190  // WS Generic Input:
191  const GriddedField2& gfraw_in,
192  const Verbosity&)
193 {
194  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 0, 1, gfraw_in);
195 
196  field_out = gfraw_in.data;
197 }
198 
199 
200 /* Workspace method: Doxygen documentation will be auto-generated */
201 void FieldFromGriddedField(// WS Generic Output:
202  Tensor3& field_out,
203  // WS Input:
204  const Vector& p_grid,
205  const Vector& lat_grid,
206  const Vector& lon_grid,
207  // WS Generic Input:
208  const GriddedField3& gfraw_in,
209  const Verbosity&)
210 {
211  chk_griddedfield_gridname(gfraw_in, 0, "Pressure");
212  chk_if_equal("p_grid", "gfield.p_grid", p_grid, gfraw_in.get_numeric_grid(0));
213 
214  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 1, 2, gfraw_in);
215 
216  field_out = gfraw_in.data;
217 }
218 
219 
220 /* Workspace method: Doxygen documentation will be auto-generated */
221 void FieldFromGriddedField(// WS Generic Output:
222  Tensor4& field_out,
223  // WS Input:
224  const Vector& p_grid,
225  const Vector& lat_grid,
226  const Vector& lon_grid,
227  // WS Generic Input:
228  const GriddedField4& gfraw_in,
229  const Verbosity&)
230 {
231  chk_griddedfield_gridname(gfraw_in, 1, "Pressure");
232  chk_if_equal("p_grid", "gfield.p_grid", p_grid, gfraw_in.get_numeric_grid(0));
233 
234  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 2, 3, gfraw_in);
235 
236  field_out = gfraw_in.data;
237 }
238 
239 
240 /* Workspace method: Doxygen documentation will be auto-generated */
241 void FieldFromGriddedField(// WS Generic Output:
242  Tensor4& field_out,
243  // WS Input:
244  const Vector& p_grid,
245  const Vector& lat_grid,
246  const Vector& lon_grid,
247  // WS Generic Input:
248  const ArrayOfGriddedField3& gfraw_in,
249  const Verbosity& verbosity)
250 {
251  if (!gfraw_in.nelem())
252  {
253  CREATE_OUT1;
254  out1 << " Warning: gfraw_in is empty, proceeding anyway\n";
255  field_out.resize(0, 0, 0, 0);
256  }
257  else
258  {
259  field_out.resize(gfraw_in.nelem(), p_grid.nelem(),
260  gfraw_in[0].data.nrows(), gfraw_in[0].data.ncols());
261  }
262 
263  for (Index i = 0; i < gfraw_in.nelem(); i++)
264  {
265  chk_griddedfield_gridname(gfraw_in[i], 0, "Pressure");
266  chk_if_equal("p_grid", "gfield.p_grid", p_grid, gfraw_in[i].get_numeric_grid(0));
267 
268  FieldFromGriddedFieldCheckLatLonHelper(lat_grid, lon_grid, 1, 2, gfraw_in[i]);
269 
270  field_out(i, joker, joker, joker) = gfraw_in[i].data;
271  }
272 }
273 
274 
275 /* Workspace method: Doxygen documentation will be auto-generated */
276 void GriddedFieldLatLonExpand(// WS Generic Output:
277  GriddedField2& gfraw_out,
278  // WS Generic Input:
279  const GriddedField2& gfraw_in_orig,
280  const Verbosity&)
281 {
282  const GriddedField2* gfraw_in_pnt;
283  GriddedField2 gfraw_in_copy;
284 
285  if (&gfraw_in_orig == &gfraw_out)
286  {
287  gfraw_in_copy = gfraw_in_orig;
288  gfraw_in_pnt = &gfraw_in_copy;
289  }
290  else
291  gfraw_in_pnt = &gfraw_in_orig;
292 
293  const GriddedField2 &gfraw_in = *gfraw_in_pnt;
294 
295  chk_griddedfield_gridname(gfraw_in, 0, "Latitude");
296  chk_griddedfield_gridname(gfraw_in, 1, "Longitude");
297 
298  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
299  throw runtime_error("Can't expand data because number of Latitudes and Longitudes is greater than 1");
300 
301  gfraw_out.set_grid_name(0, "Latitude");
302  gfraw_out.set_grid_name(1, "Longitude");
303 
304  Vector v(2);
305  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1)
306  {
307  v[0] = -90; v[1] = 90;
308  gfraw_out.set_grid(0, v);
309  gfraw_out.resize(2, gfraw_in.data.ncols());
310 
311  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
312  gfraw_out.data(joker, j) = gfraw_in.data(0, j);
313  }
314  else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1)
315  {
316  v[0] = 0; v[1] = 360;
317  gfraw_out.set_grid(1, v);
318  gfraw_out.resize(gfraw_in.data.nrows(), 2);
319 
320  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
321  gfraw_out.data(j, joker) = gfraw_in.data(j, 0);
322  }
323  else
324  {
325  v[0] = -90; v[1] = 90;
326  gfraw_out.set_grid(0, v);
327  v[0] = 0; v[1] = 360;
328  gfraw_out.set_grid(1, v);
329  gfraw_out.resize(2, 2);
330 
331  gfraw_out.data = gfraw_in.data(0, 0);
332  }
333 }
334 
335 
336 /* Workspace method: Doxygen documentation will be auto-generated */
337 void GriddedFieldLatLonExpand(// WS Generic Output:
338  GriddedField3& gfraw_out,
339  // WS Generic Input:
340  const GriddedField3& gfraw_in_orig,
341  const Verbosity&)
342 {
343  const GriddedField3* gfraw_in_pnt;
344  GriddedField3 gfraw_in_copy;
345 
346  if (&gfraw_in_orig == &gfraw_out)
347  {
348  gfraw_in_copy = gfraw_in_orig;
349  gfraw_in_pnt = &gfraw_in_copy;
350  }
351  else
352  gfraw_in_pnt = &gfraw_in_orig;
353 
354  const GriddedField3 &gfraw_in = *gfraw_in_pnt;
355 
356  chk_griddedfield_gridname(gfraw_in, 1, "Latitude");
357  chk_griddedfield_gridname(gfraw_in, 2, "Longitude");
358 
359  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
360  throw runtime_error("Can't expand data because number of Latitudes and Longitudes is greater than 1");
361 
362  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
363  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
364  gfraw_out.set_grid_name(1, "Latitude");
365  gfraw_out.set_grid_name(2, "Longitude");
366 
367  Vector v(2);
368  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1)
369  {
370  v[0] = -90; v[1] = 90;
371  gfraw_out.set_grid(1, v);
372  gfraw_out.resize(gfraw_in.data.npages(), 2, gfraw_in.data.ncols());
373 
374  for (Index i = 0; i < gfraw_in.data.npages(); i++)
375  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
376  gfraw_out.data(i, joker, j) = gfraw_in.data(i, 0, j);
377  }
378  else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1)
379  {
380  v[0] = 0; v[1] = 360;
381  gfraw_out.set_grid(2, v);
382  gfraw_out.resize(gfraw_in.data.npages(), gfraw_in.data.nrows(), 2);
383 
384  for (Index i = 0; i < gfraw_in.data.npages(); i++)
385  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
386  gfraw_out.data(i, j, joker) = gfraw_in.data(i, j, 0);
387  }
388  else
389  {
390  v[0] = -90; v[1] = 90;
391  gfraw_out.set_grid(1, v);
392  v[0] = 0; v[1] = 360;
393  gfraw_out.set_grid(2, v);
394  gfraw_out.resize(gfraw_in.data.npages(), 2, 2);
395 
396  for (Index i = 0; i < gfraw_in.data.npages(); i++)
397  gfraw_out.data(i, joker, joker) = gfraw_in.data(i, 0, 0);
398  }
399 }
400 
401 /* Workspace method: Doxygen documentation will be auto-generated */
402 void GriddedFieldLatLonExpand(// WS Generic Output:
403  GriddedField4& gfraw_out,
404  // WS Generic Input:
405  const GriddedField4& gfraw_in_orig,
406  const Verbosity&)
407 {
408  const GriddedField4* gfraw_in_pnt;
409  GriddedField4 gfraw_in_copy;
410 
411  if (&gfraw_in_orig == &gfraw_out)
412  {
413  gfraw_in_copy = gfraw_in_orig;
414  gfraw_in_pnt = &gfraw_in_copy;
415  }
416  else
417  gfraw_in_pnt = &gfraw_in_orig;
418 
419  const GriddedField4 &gfraw_in = *gfraw_in_pnt;
420 
421  chk_griddedfield_gridname(gfraw_in, 2, "Latitude");
422  chk_griddedfield_gridname(gfraw_in, 3, "Longitude");
423 
424  if (gfraw_in.data.ncols() != 1 && gfraw_in.data.nrows() != 1)
425  throw runtime_error("Can't expand data because number of Latitudes and Longitudes is greater than 1");
426 
427  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
428  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
429 
430  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
431  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
432 
433  gfraw_out.set_grid_name(2, "Latitude");
434  gfraw_out.set_grid_name(3, "Longitude");
435 
436  Vector v(2);
437  if (gfraw_in.data.nrows() == 1 && gfraw_in.data.ncols() != 1)
438  {
439  v[0] = -90; v[1] = 90;
440  gfraw_out.set_grid(2, v);
441  gfraw_out.resize(gfraw_in.data.nbooks(), gfraw_in.data.npages(), 2, gfraw_in.data.ncols());
442 
443  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
444  for (Index i = 0; i < gfraw_in.data.npages(); i++)
445  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
446  gfraw_out.data(k, i, joker, j) = gfraw_in.data(k, i, 0, j);
447  }
448  else if (gfraw_in.data.nrows() != 1 && gfraw_in.data.ncols() == 1)
449  {
450  v[0] = 0; v[1] = 360;
451  gfraw_out.set_grid(3, v);
452  gfraw_out.resize(gfraw_in.data.nbooks(), gfraw_in.data.npages(), gfraw_in.data.nrows(), 2);
453 
454  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
455  for (Index i = 0; i < gfraw_in.data.npages(); i++)
456  for (Index j = 0; j < gfraw_in.data.nrows(); j++)
457  gfraw_out.data(k, i, j, joker) = gfraw_in.data(k, i, j, 0);
458  }
459  else
460  {
461  v[0] = -90; v[1] = 90;
462  gfraw_out.set_grid(2, v);
463  v[0] = 0; v[1] = 360;
464  gfraw_out.set_grid(3, v);
465  gfraw_out.resize(gfraw_in.data.nbooks(), gfraw_in.data.npages(), 2, 2);
466 
467  for (Index k = 0; k < gfraw_in.data.nbooks(); k++)
468  for (Index i = 0; i < gfraw_in.data.npages(); i++)
469  gfraw_out.data(k, i, joker, joker) = gfraw_in.data(k, i, 0, 0);
470  }
471 }
472 
473 
474 /* Workspace method: Doxygen documentation will be auto-generated */
475 void GriddedFieldLatLonExpand(// WS Generic Output:
476  ArrayOfGriddedField3& gfraw_out,
477  // WS Generic Input:
478  const ArrayOfGriddedField3& gfraw_in,
479  const Verbosity& verbosity)
480 {
481  gfraw_out.resize(gfraw_in.nelem());
482 
483  for (Index i = 0; i < gfraw_in.nelem(); i++)
484  GriddedFieldLatLonExpand(gfraw_out[i], gfraw_in[i], verbosity);
485 }
486 
487 
489 /*
490  This helper function is used by all GriddedFieldLatLonRegrid WSMs to do the common
491  calculation of the grid positions and interpolation weights for latitudes and longitudes.
492 
493  \param[out] ing_min Index in the new grid with first value covered
494  by the old grid.
495  \param[out] ing_max Index in the new grid with last value covered
496  by the old grid.
497  \param[out] gp_p Pressure grid positions
498  \param[out] itw Interpolation weights
499  \param[in,out] gfraw_out Output GriddedField
500  \param[in] gfraw_in Input GriddedField
501  \param[in] p_grid_index Index of pressure grid
502  \param[in] p_grid New pressure grid
503  \param[in] interp_order Interpolation order
504  \param[in] zeropadding Allow zero padding
505  \param[in] verbosity Verbosity levels
506  */
507 void GriddedFieldPRegridHelper(Index& ing_min,
508  Index& ing_max,
509  ArrayOfGridPosPoly& gp_p,
510  Matrix& itw,
511  GriddedField& gfraw_out,
512  const GriddedField& gfraw_in,
513  const Index p_grid_index,
514  ConstVectorView p_grid,
515  const Index& interp_order,
516  const Index& zeropadding,
517  const Verbosity& verbosity)
518 {
519  CREATE_OUT2;
520 
521  chk_griddedfield_gridname(gfraw_in, p_grid_index, "Pressure");
522 
523  out2 << " Interpolation order: " << interp_order << "\n";
524 
525  const ConstVectorView in_p_grid = gfraw_in.get_numeric_grid(p_grid_index);
526 
527  // Initialize output field. Set grids and copy grid names
528  gfraw_out.set_grid(p_grid_index, p_grid);
529  gfraw_out.set_grid_name(p_grid_index, gfraw_in.get_grid_name(p_grid_index));
530 
531  if (zeropadding)
532  {
533  if (in_p_grid[0]<p_grid[p_grid.nelem()-1] ||
534  in_p_grid[in_p_grid.nelem()-1]>p_grid[0])
535  {
536  ing_min = 0;
537  ing_max = ing_min-1;
538  }
539  else
540  chk_interpolation_pgrids_loose_no_data_check(ing_min, ing_max,
541  "Raw field to p_grid",
542  in_p_grid,
543  p_grid,
544  interp_order);
545  }
546  else
547  {
548  ing_min = 0;
549  ing_max = p_grid.nelem()-1;
550  chk_interpolation_pgrids("Raw field to p_grid",
551  in_p_grid,
552  p_grid,
553  interp_order);
554  }
555  Index nelem_in_range = ing_max - ing_min + 1;
556 
557  // Calculate grid positions:
558  if (nelem_in_range>0)
559  {
560  gp_p.resize(nelem_in_range);
561  p2gridpos_poly(gp_p, in_p_grid, p_grid[Range(ing_min, nelem_in_range)],
562  interp_order);
563 
564  // Interpolation weights:
565  itw.resize(nelem_in_range, interp_order+1);
566  interpweights(itw, gp_p);
567  }
568 }
569 
570 
571 /* Workspace method: Doxygen documentation will be auto-generated */
572 void GriddedFieldPRegrid(// WS Generic Output:
573  GriddedField3& gfraw_out,
574  // WS Input:
575  const Vector& p_grid,
576  // WS Generic Input:
577  const GriddedField3& gfraw_in_orig,
578  const Index& interp_order,
579  const Index& zeropadding,
580  const Verbosity& verbosity)
581 {
582  const GriddedField3* gfraw_in_pnt;
583  GriddedField3 gfraw_in_copy;
584 
585  if (&gfraw_in_orig == &gfraw_out)
586  {
587  gfraw_in_copy = gfraw_in_orig;
588  gfraw_in_pnt = &gfraw_in_copy;
589  }
590  else
591  gfraw_in_pnt = &gfraw_in_orig;
592 
593  const GriddedField3 &gfraw_in = *gfraw_in_pnt;
594 
595  const Index p_grid_index = 0;
596 
597  // Resize output GriddedField and copy all non-latitude/longitude grids
598  gfraw_out.resize(p_grid.nelem(), gfraw_in.data.nrows(), gfraw_in.data.ncols());
599  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
600  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
601  gfraw_out.set_grid(2, gfraw_in.get_numeric_grid(2));
602  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
603 
604  ArrayOfGridPosPoly gp_p;
605  Matrix itw;
606 
607  Index ing_min, ing_max;
608 
609  GriddedFieldPRegridHelper(ing_min, ing_max,
610  gp_p, itw, gfraw_out,
611  gfraw_in, p_grid_index,
612  p_grid, interp_order, zeropadding,
613  verbosity);
614 
615  // Interpolate:
616  if (ing_max - ing_min < 0)
617  gfraw_out.data = 0.;
618  else
619  if (ing_max - ing_min + 1 != p_grid.nelem())
620  {
621  gfraw_out.data = 0.;
622  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
623  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
624  {
625 // chk_interpolation_grids_loose_check_data(ing_min, ing_max,
626 // "Raw field to p_grid",
627 // gfraw_in.get_numeric_grid(p_grid_index),
628 // p_grid, gfraw_in.data(joker, i, j));
629  interp(gfraw_out.data(Range(ing_min, ing_max-ing_min+1), i, j),
630  itw, gfraw_in.data(joker, i, j), gp_p);
631  }
632  }
633  else
634  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
635  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
636  interp(gfraw_out.data(joker, i, j),
637  itw, gfraw_in.data(joker, i, j), gp_p);
638 }
639 
640 
641 /* Workspace method: Doxygen documentation will be auto-generated */
642 void GriddedFieldPRegrid(// WS Generic Output:
643  GriddedField4& gfraw_out,
644  // WS Input:
645  const Vector& p_grid,
646  // WS Generic Input:
647  const GriddedField4& gfraw_in_orig,
648  const Index& interp_order,
649  const Index& zeropadding,
650  const Verbosity& verbosity)
651 {
652  const GriddedField4* gfraw_in_pnt;
653  GriddedField4 gfraw_in_copy;
654 
655  if (&gfraw_in_orig == &gfraw_out)
656  {
657  gfraw_in_copy = gfraw_in_orig;
658  gfraw_in_pnt = &gfraw_in_copy;
659  }
660  else
661  gfraw_in_pnt = &gfraw_in_orig;
662 
663  const GriddedField4 &gfraw_in = *gfraw_in_pnt;
664 
665  const Index p_grid_index = 1;
666 
667  // Resize output GriddedField and copy all non-latitude/longitude grids
668  gfraw_out.resize(gfraw_in.data.nbooks(), p_grid.nelem(),
669  gfraw_in.data.nrows(), gfraw_in.data.ncols());
670  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
671  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
672  gfraw_out.set_grid(2, gfraw_in.get_numeric_grid(2));
673  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
674  gfraw_out.set_grid(3, gfraw_in.get_numeric_grid(3));
675  gfraw_out.set_grid_name(3, gfraw_in.get_grid_name(3));
676 
677  ArrayOfGridPosPoly gp_p;
678  Matrix itw;
679 
680  Index ing_min, ing_max;
681 
682  GriddedFieldPRegridHelper(ing_min, ing_max,
683  gp_p, itw, gfraw_out,
684  gfraw_in, p_grid_index,
685  p_grid, interp_order, zeropadding,
686  verbosity);
687 
688  // Interpolate:
689  if (ing_max - ing_min < 0)
690  gfraw_out.data = 0.;
691  else
692  if (ing_max - ing_min + 1 != p_grid.nelem())
693  {
694  gfraw_out.data = 0.;
695  for (Index b = 0; b < gfraw_in.data.nbooks(); b++)
696  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
697  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
698  {
699 // chk_interpolation_grids_loose_check_data(ing_min, ing_max,
700 // "Raw field to p_grid",
701 // gfraw_in.get_numeric_grid(p_grid_index),
702 // p_grid, gfraw_in.data(b, joker, i, j));
703  interp(gfraw_out.data(b, Range(ing_min, ing_max-ing_min), i, j),
704  itw, gfraw_in.data(b, joker, i, j), gp_p);
705  }
706  }
707  else
708  for (Index b = 0; b < gfraw_in.data.nbooks(); b++)
709  for (Index i = 0; i < gfraw_in.data.nrows(); i++)
710  for (Index j = 0; j < gfraw_in.data.ncols(); j++)
711  interp(gfraw_out.data(b, joker, i, j),
712  itw, gfraw_in.data(b, joker, i, j), gp_p);
713 }
714 
715 
716 /* Workspace method: Doxygen documentation will be auto-generated */
717 void GriddedFieldPRegrid(// WS Generic Output:
718  ArrayOfGriddedField3& agfraw_out,
719  // WS Input:
720  const Vector& p_grid,
721  // WS Generic Input:
722  const ArrayOfGriddedField3& agfraw_in,
723  const Index& interp_order,
724  const Index& zeropadding,
725  const Verbosity& verbosity)
726 {
727  agfraw_out.resize(agfraw_in.nelem());
728 
729  for (Index i = 0; i < agfraw_in.nelem(); i++)
730  {
731  GriddedFieldPRegrid(agfraw_out[i], p_grid, agfraw_in[i],
732  interp_order, zeropadding,
733  verbosity);
734  }
735 }
736 
737 
739 /*
740  This helper function is used by all GriddedFieldLatLonRegrid WSMs to do the common
741  calculation of the grid positions and interpolation weights for latitudes and longitudes.
742 
743  \param[out] gp_lat Latitude grid positions
744  \param[out] gp_lon Longitude grid positions
745  \param[out] itw Interpolation weights
746  \param[in,out] gfraw_out Output GriddedField
747  \param[in] gfraw_in Input GriddedField
748  \param[in] lat_grid_index Index of latitude grid
749  \param[in] lon_grid_index Index of longitude grid
750  \param[in] lat_true New latitude grid
751  \param[in] lon_true New longitude grid
752  \param[in] interp_order Interpolation order
753  \param[in] verbosity Verbosity levels
754  */
755 void GriddedFieldLatLonRegridHelper(ArrayOfGridPosPoly& gp_lat,
756  ArrayOfGridPosPoly& gp_lon,
757  Tensor3& itw,
758  GriddedField& gfraw_out,
759  const GriddedField& gfraw_in,
760  const Index lat_grid_index,
761  const Index lon_grid_index,
762  ConstVectorView lat_true,
763  ConstVectorView lon_true,
764  const Index& interp_order,
765  const Verbosity& verbosity)
766 {
767  CREATE_OUT2;
768 
769  if (!lat_true.nelem()) throw runtime_error("The new latitude grid is not allowed to be empty.");
770  if (!lon_true.nelem()) throw runtime_error("The new longitude grid is not allowed to be empty.");
771 
772  chk_griddedfield_gridname(gfraw_in, lat_grid_index, "Latitude");
773  chk_griddedfield_gridname(gfraw_in, lon_grid_index, "Longitude");
774  if (gfraw_in.get_grid_size(lat_grid_index) == 1 || gfraw_in.get_grid_size(lon_grid_index) == 1)
775  throw runtime_error("Raw data has to be true 3D data (nlat>1 and nlon>1).");
776 
777  out2 << " Interpolation order: " << interp_order << "\n";
778 
779  const ConstVectorView in_lat_grid = gfraw_in.get_numeric_grid(lat_grid_index);
780  const ConstVectorView in_lon_grid = gfraw_in.get_numeric_grid(lon_grid_index);
781 
782  // Initialize output field. Set grids and copy grid names
783  gfraw_out.set_grid(lat_grid_index, lat_true);
784  gfraw_out.set_grid_name(lat_grid_index, gfraw_in.get_grid_name(lat_grid_index));
785  gfraw_out.set_grid(lon_grid_index, lon_true);
786  gfraw_out.set_grid_name(lon_grid_index, gfraw_in.get_grid_name(lon_grid_index));
787 
788  chk_interpolation_grids("Raw field to lat_grid, 3D case",
789  in_lat_grid,
790  lat_true,
791  interp_order);
792 
793  // Calculate grid positions:
794  gp_lat.resize(lat_true.nelem());
795  gp_lon.resize(lon_true.nelem());
796  gridpos_poly(gp_lat, in_lat_grid, lat_true, interp_order);
797 
798  gridpos_poly_longitudinal("Raw field to lon_grid, 3D case",
799  gp_lon, in_lon_grid, lon_true, interp_order);
800 
801  // Interpolation weights:
802  itw.resize(lat_true.nelem(), lon_true.nelem(), (interp_order+1)*(interp_order+1));
803  interpweights(itw, gp_lat, gp_lon);
804 }
805 
806 
807 /* Workspace method: Doxygen documentation will be auto-generated */
808 void GriddedFieldLatLonRegrid(// WS Generic Output:
809  GriddedField2& gfraw_out,
810  // WS Input:
811  const Vector& lat_true,
812  const Vector& lon_true,
813  // WS Generic Input:
814  const GriddedField2& gfraw_in_orig,
815  const Index& interp_order,
816  const Verbosity& verbosity)
817 {
818  if (!lat_true.nelem()) throw runtime_error("The new latitude grid is not allowed to be empty.");
819  if (!lon_true.nelem()) throw runtime_error("The new longitude grid is not allowed to be empty.");
820 
821  const GriddedField2* gfraw_in_pnt;
822  GriddedField2 gfraw_in_copy;
823 
824  if (&gfraw_in_orig == &gfraw_out)
825  {
826  gfraw_in_copy = gfraw_in_orig;
827  gfraw_in_pnt = &gfraw_in_copy;
828  }
829  else
830  gfraw_in_pnt = &gfraw_in_orig;
831 
832  const GriddedField2 &gfraw_in = *gfraw_in_pnt;
833 
834  const Index lat_grid_index = 0;
835  const Index lon_grid_index = 1;
836 
837  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
838  gfraw_in.get_grid_size(lon_grid_index) < 2)
839  {
840  ostringstream os;
841  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
842  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
843  throw runtime_error(os.str());
844  }
845 
846  // Resize output GriddedField
847  gfraw_out.resize(lat_true.nelem(), lon_true.nelem());
848 
849  ArrayOfGridPosPoly gp_lat;
850  ArrayOfGridPosPoly gp_lon;
851  Tensor3 itw;
852 
853  // If lon grid is cyclic, the data values at 0 and 360 must match
854  const ConstVectorView& in_lat_grid = gfraw_in.get_numeric_grid(lat_grid_index);
855  const ConstVectorView& in_lon_grid = gfraw_in.get_numeric_grid(lon_grid_index);
856 
857  if (is_lon_cyclic(in_lon_grid))
858  {
859  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++)
860  {
861  if (!is_same_within_epsilon(gfraw_in.data(lat, 0),
862  gfraw_in.data(lat, in_lon_grid.nelem()-1),
863  EPSILON_LON_CYCLIC))
864  {
865  ostringstream os;
866  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
867  << "Mismatch at latitude index : "
868  << lat << " (" << in_lat_grid[lat] << " degrees)\n"
869  << "Value at 0 degrees longitude : "
870  << gfraw_in.data(lat, 0) << "\n"
871  << "Value at 360 degrees longitude: "
872  << gfraw_in.data(lat, in_lon_grid.nelem()-1) << "\n"
873  << "Difference : "
874  << gfraw_in.data(lat, in_lon_grid.nelem()-1) - gfraw_in.data(lat, 0) << "\n"
875  << "Allowed difference : " << EPSILON_LON_CYCLIC;
876  throw runtime_error(os.str());
877  }
878  }
879  }
880 
881  GriddedFieldLatLonRegridHelper(gp_lat, gp_lon, itw, gfraw_out,
882  gfraw_in, lat_grid_index, lon_grid_index,
883  lat_true, lon_true, interp_order,
884  verbosity);
885 
886  // Interpolate:
887  interp(gfraw_out.data, itw, gfraw_in.data, gp_lat, gp_lon);
888 }
889 
890 
891 /* Workspace method: Doxygen documentation will be auto-generated */
892 void GriddedFieldLatLonRegrid(// WS Generic Output:
893  GriddedField3& gfraw_out,
894  // WS Input:
895  const Vector& lat_true,
896  const Vector& lon_true,
897  // WS Generic Input:
898  const GriddedField3& gfraw_in_orig,
899  const Index& interp_order,
900  const Verbosity& verbosity)
901 {
902  if (!lat_true.nelem()) throw runtime_error("The new latitude grid is not allowed to be empty.");
903  if (!lon_true.nelem()) throw runtime_error("The new longitude grid is not allowed to be empty.");
904 
905  const GriddedField3* gfraw_in_pnt;
906  GriddedField3 gfraw_in_copy;
907 
908  if (&gfraw_in_orig == &gfraw_out)
909  {
910  gfraw_in_copy = gfraw_in_orig;
911  gfraw_in_pnt = &gfraw_in_copy;
912  }
913  else
914  gfraw_in_pnt = &gfraw_in_orig;
915 
916  const GriddedField3 &gfraw_in = *gfraw_in_pnt;
917 
918  const Index lat_grid_index = 1;
919  const Index lon_grid_index = 2;
920 
921  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
922  gfraw_in.get_grid_size(lon_grid_index) < 2)
923  {
924  ostringstream os;
925  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
926  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
927  throw runtime_error(os.str());
928  }
929 
930  // Resize output GriddedField and copy all non-latitude/longitude grids
931  gfraw_out.resize(gfraw_in.data.npages(), lat_true.nelem(), lon_true.nelem());
932  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
933  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
934 
935  ArrayOfGridPosPoly gp_lat;
936  ArrayOfGridPosPoly gp_lon;
937  Tensor3 itw;
938 
939  // If lon grid is cyclic, the data values at 0 and 360 must match
940  const ConstVectorView& in_grid0 = gfraw_in.get_numeric_grid(0);
941  const ConstVectorView& in_lat_grid = gfraw_in.get_numeric_grid(lat_grid_index);
942  const ConstVectorView& in_lon_grid = gfraw_in.get_numeric_grid(lon_grid_index);
943 
944  if (is_lon_cyclic(in_lon_grid))
945  {
946  for (Index g0 = 0; g0 < in_grid0.nelem(); g0++)
947  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++)
948  {
949  if (!is_same_within_epsilon(gfraw_in.data(g0, lat, 0),
950  gfraw_in.data(g0, lat, in_lon_grid.nelem()-1),
951  EPSILON_LON_CYCLIC))
952  {
953  ostringstream os;
954  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
955  << "Mismatch at 1st grid index : "
956  << g0 << " (" << in_grid0[g0] << ")\n"
957  << " at latitude index : "
958  << lat << " (" << in_lat_grid[lat] << " degrees)\n"
959  << "Value at 0 degrees longitude : "
960  << gfraw_in.data(g0, lat, 0) << "\n"
961  << "Value at 360 degrees longitude: "
962  << gfraw_in.data(g0, lat, in_lon_grid.nelem()-1) << "\n"
963  << "Difference : "
964  << gfraw_in.data(g0, lat, in_lon_grid.nelem()-1) - gfraw_in.data(g0, lat, 0) << "\n"
965  << "Allowed difference : " << EPSILON_LON_CYCLIC;
966  throw runtime_error(os.str());
967  }
968  }
969  }
970 
971  GriddedFieldLatLonRegridHelper(gp_lat, gp_lon, itw, gfraw_out,
972  gfraw_in, lat_grid_index, lon_grid_index,
973  lat_true, lon_true, interp_order,
974  verbosity);
975 
976  // Interpolate:
977  for (Index i = 0; i < gfraw_in.data.npages(); i++)
978  interp(gfraw_out.data(i, joker, joker),
979  itw, gfraw_in.data(i, joker, joker), gp_lat, gp_lon);
980 }
981 
982 
983 /* Workspace method: Doxygen documentation will be auto-generated */
984 void GriddedFieldLatLonRegrid(// WS Generic Output:
985  GriddedField4& gfraw_out,
986  // WS Input:
987  const Vector& lat_true,
988  const Vector& lon_true,
989  // WS Generic Input:
990  const GriddedField4& gfraw_in_orig,
991  const Index& interp_order,
992  const Verbosity& verbosity)
993 {
994  if (!lat_true.nelem()) throw runtime_error("The new latitude grid is not allowed to be empty.");
995  if (!lon_true.nelem()) throw runtime_error("The new longitude grid is not allowed to be empty.");
996 
997  const GriddedField4* gfraw_in_pnt;
998  GriddedField4 gfraw_in_copy;
999 
1000  if (&gfraw_in_orig == &gfraw_out)
1001  {
1002  gfraw_in_copy = gfraw_in_orig;
1003  gfraw_in_pnt = &gfraw_in_copy;
1004  }
1005  else
1006  gfraw_in_pnt = &gfraw_in_orig;
1007 
1008  const GriddedField4 &gfraw_in = *gfraw_in_pnt;
1009 
1010  const Index lat_grid_index = 2;
1011  const Index lon_grid_index = 3;
1012 
1013  if (gfraw_in.get_grid_size(lat_grid_index) < 2 ||
1014  gfraw_in.get_grid_size(lon_grid_index) < 2)
1015  {
1016  ostringstream os;
1017  os << "Raw data has to be true 3D data (nlat>1 and nlon>1).\n"
1018  << "Use GriddedFieldLatLonExpand to convert 1D or 2D data to 3D!\n";
1019  throw runtime_error(os.str());
1020  }
1021 
1022  // Resize output GriddedField and copy all non-latitude/longitude grids
1023  gfraw_out.resize(gfraw_in.data.nbooks(), gfraw_in.data.npages(),
1024  lat_true.nelem(), lon_true.nelem());
1025  gfraw_out.set_grid(0, gfraw_in.get_numeric_grid(0));
1026  gfraw_out.set_grid_name(0, gfraw_in.get_grid_name(0));
1027  gfraw_out.set_grid(1, gfraw_in.get_numeric_grid(1));
1028  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
1029 
1030  ArrayOfGridPosPoly gp_lat;
1031  ArrayOfGridPosPoly gp_lon;
1032  Tensor3 itw;
1033 
1034  GriddedFieldLatLonRegridHelper(gp_lat, gp_lon, itw, gfraw_out,
1035  gfraw_in, lat_grid_index, lon_grid_index,
1036  lat_true, lon_true, interp_order,
1037  verbosity);
1038 
1039  // If lon grid is cyclic, the data values at 0 and 360 must match
1040  const ConstVectorView& in_grid0 = gfraw_in.get_numeric_grid(0);
1041  const ConstVectorView& in_grid1 = gfraw_in.get_numeric_grid(1);
1042  const ConstVectorView& in_lat_grid = gfraw_in.get_numeric_grid(lat_grid_index);
1043  const ConstVectorView& in_lon_grid = gfraw_in.get_numeric_grid(lon_grid_index);
1044 
1045  if (is_lon_cyclic(in_lon_grid))
1046  {
1047  for (Index g0 = 0; g0 < in_grid0.nelem(); g0++)
1048  for (Index g1 = 0; g1 < in_grid1.nelem(); g1++)
1049  for (Index lat = 0; lat < in_lat_grid.nelem(); lat++)
1050  {
1051  if (!is_same_within_epsilon(gfraw_in.data(g0, g1, lat, 0),
1052  gfraw_in.data(g0, g1, lat, in_lon_grid.nelem()-1),
1053  EPSILON_LON_CYCLIC))
1054  {
1055  ostringstream os;
1056  os << "Data values at 0 and 360 degrees for a cyclic longitude grid must match: \n"
1057  << "Mismatch at 1st grid index : "
1058  << g0 << " (" << in_grid0[g0] << ")\n"
1059  << " at 2nd grid index : "
1060  << g1 << " (" << in_grid1[g1] << ")\n"
1061  << " at latitude index : "
1062  << lat << " (" << in_lat_grid[lat] << " degrees)\n"
1063  << "Value at 0 degrees longitude : "
1064  << gfraw_in.data(g0, g1, lat, 0) << "\n"
1065  << "Value at 360 degrees longitude: "
1066  << gfraw_in.data(g0, g1, lat, in_lon_grid.nelem()-1) << "\n"
1067  << "Difference : "
1068  << gfraw_in.data(g0, g1, lat, in_lon_grid.nelem()-1) - gfraw_in.data(g0, g1, lat, 0) << "\n"
1069  << "Allowed difference : " << EPSILON_LON_CYCLIC;
1070  throw runtime_error(os.str());
1071  }
1072  }
1073  }
1074 
1075  // Interpolate:
1076  for (Index i = 0; i < gfraw_in.data.nbooks(); i++)
1077  for (Index j = 0; j < gfraw_in.data.npages(); j++)
1078  interp(gfraw_out.data(i, j, joker, joker),
1079  itw, gfraw_in.data(i, j, joker, joker), gp_lat, gp_lon);
1080 }
1081 
1082 
1083 /* Workspace method: Doxygen documentation will be auto-generated */
1084 void GriddedFieldLatLonRegrid(// WS Generic Output:
1085  ArrayOfGriddedField3& agfraw_out,
1086  // WS Input:
1087  const Vector& lat_true,
1088  const Vector& lon_true,
1089  // WS Generic Input:
1090  const ArrayOfGriddedField3& agfraw_in,
1091  const Index& interp_order,
1092  const Verbosity& verbosity)
1093 {
1094  agfraw_out.resize(agfraw_in.nelem());
1095 
1096  for (Index i = 0; i < agfraw_in.nelem(); i++)
1097  {
1098  GriddedFieldLatLonRegrid(agfraw_out[i], lat_true, lon_true, agfraw_in[i],
1099  interp_order, verbosity);
1100  }
1101 }
1102 
1103 
1105 /*
1106  This helper function is used by GriddedFieldZToPRegrid WSM to do the common
1107  calculation of the grid positions and interpolation weights for latitudes and longitudes.
1108 
1109  \param[out] ing_min Index in the new grid with first value covered
1110  by the old grid.
1111  \param[out] ing_max Index in the new grid with last value covered
1112  by the old grid.
1113  \param[out] gp_p Altitude grid positions
1114  \param[out] itw Interpolation weights
1115  \param[in] gfraw_in Input GriddedField
1116  \param[in] z_grid_index Index of altitude grid
1117  \param[in] z_grid New z_grid grid
1118  \param[in] interp_order Interpolation order
1119  \param[in] zeropadding Allow zero padding
1120  \param[in] verbosity Verbosity levels
1121  */
1122 void GriddedFieldZToPRegridHelper(Index& ing_min,
1123  Index& ing_max,
1124  ArrayOfGridPosPoly& gp_p,
1125  Matrix& itw,
1126  const GriddedField& gfraw_in,
1127  const Index z_grid_index,
1128  ConstVectorView z_grid,
1129  const Index& interp_order,
1130  const Index& zeropadding,
1131  const Verbosity& verbosity)
1132 {
1133  CREATE_OUT2;
1134 
1135  chk_griddedfield_gridname(gfraw_in, z_grid_index, "Altitude");
1136 
1137  out2 << " Interpolation order: " << interp_order << "\n";
1138 
1139  const ConstVectorView in_z_grid = gfraw_in.get_numeric_grid(z_grid_index);
1140 
1141  if (zeropadding)
1142  {
1143  if (in_z_grid[0] > z_grid[z_grid.nelem()-1] ||
1144  in_z_grid[in_z_grid.nelem()-1] < z_grid[0])
1145  {
1146  ing_min = 0;
1147  ing_max = ing_min-1;
1148  }
1149  else
1150  chk_interpolation_grids_loose_no_data_check(ing_min, ing_max,
1151  "Raw field to z_grid",
1152  in_z_grid,
1153  z_grid,
1154  interp_order);
1155  }
1156  else
1157  {
1158  ing_min = 0;
1159  ing_max = z_grid.nelem()-1;
1160  chk_interpolation_pgrids("Raw field to p_grid",
1161  in_z_grid,
1162  z_grid,
1163  interp_order);
1164  }
1165 
1166  Index nelem_in_range = ing_max - ing_min + 1;
1167 
1168  // Calculate grid positions:
1169  if (nelem_in_range>0)
1170  {
1171  gp_p.resize(nelem_in_range);
1172  gridpos_poly(gp_p, in_z_grid, z_grid[Range(ing_min, nelem_in_range)],
1173  interp_order);
1174 
1175  // Interpolation weights:
1176  itw.resize(nelem_in_range, interp_order+1);
1177  interpweights(itw, gp_p);
1178  }
1179 }
1180 
1181 
1182 /* Workspace method: Doxygen documentation will be auto-generated */
1183 void GriddedFieldZToPRegrid(// WS Generic Output:
1184  GriddedField3& gfraw_out, //grid in P
1185  // WS Input:
1186  const Vector& p_grid,
1187  const Vector& lat_grid,
1188  const Vector& lon_grid,
1189  const Tensor3& z_field,
1190  // WS Generic Input:
1191  const GriddedField3& gfraw_in_orig,//grid in Z
1192  const Index& interp_order, // Only linear interpolation allowed
1193  const Index& zeropadding,
1194  const Verbosity& verbosity)
1195 {
1196 
1197  // z_field must be of the same size as its grids
1198  if(!((z_field.npages()==p_grid.nelem() && z_field.nrows() == lat_grid.nelem()) && z_field.ncols() == lon_grid.nelem()))
1199  throw std::runtime_error("*z_field* must be of the same size as *p_grid*, *lat_grid*, and *lon_grid* in *GriddedFieldZToPRegrid*.");
1200 
1201  // Must name the dimension "Altitude" to ensure user is aware of what they are doing.
1202  chk_griddedfield_gridname(gfraw_in_orig, 0, "Altitude");
1203 
1204  // Each lat and lon grid must be identical between z_field and in field
1205  const Vector& lat_in = gfraw_in_orig.get_numeric_grid(1);
1206  const Vector& lon_in = gfraw_in_orig.get_numeric_grid(2);
1207 
1208  if(lat_grid.nelem()!=lat_in.nelem() || lon_grid.nelem()!=lon_in.nelem())
1209  throw std::runtime_error("Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude to be on the same grid as *z_field*.");
1210  else
1211  {
1212  for(Index ii=0;ii<lat_grid.nelem();ii++)
1213  if(lat_grid[ii]!=lat_in[ii])
1214  throw std::runtime_error("Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude to be the same as for *z_field*.");
1215  for(Index ii=0;ii<lon_grid.nelem();ii++)
1216  if(lon_grid[ii]!=lon_in[ii])
1217  throw std::runtime_error("Gridding of field to regrid is bad.\n*GriddedFieldZToPRegrid* requires latitude and longitude to be the same as for *z_field*.");
1218  }
1219 
1220  // Pointer in case output is input variable (same memory allocated)
1221  const GriddedField3* gfraw_in_pnt;
1222  GriddedField3 gfraw_in_copy;
1223 
1224  if (&gfraw_in_orig == &gfraw_out)
1225  {
1226  gfraw_in_copy = gfraw_in_orig;
1227  gfraw_in_pnt = &gfraw_in_copy;
1228  }
1229  else
1230  gfraw_in_pnt = &gfraw_in_orig;
1231 
1232  // Now output and input are separate variables (not allocating the same memory)
1233  const GriddedField3 &gfraw_in = *gfraw_in_pnt;
1234 
1235  // Right size and order
1236  gfraw_out.resize(p_grid.nelem(),lat_grid.nelem(),lon_grid.nelem());
1237  gfraw_out.set_grid(0, p_grid);
1238  gfraw_out.set_grid_name(0, "Pressure");
1239  gfraw_out.set_grid(1, lat_grid);
1240  gfraw_out.set_grid_name(1, gfraw_in.get_grid_name(1));
1241  gfraw_out.set_grid(2, lon_grid);
1242  gfraw_out.set_grid_name(2, gfraw_in.get_grid_name(2));
1243  gfraw_out.data = 0.;
1244 
1245  ArrayOfGridPosPoly gp_p;
1246  Matrix itw;
1247 
1248  Index ing_min, ing_max;
1249 
1250  for(Index lat_index = 0;lat_index<lat_grid.nelem();lat_index++)
1251  {
1252  for(Index lon_index = 0;lon_index<lon_grid.nelem();lon_index++)
1253  {
1254  const Vector z_out = z_field(joker,lat_index,lon_index);
1255 
1256  GriddedFieldZToPRegridHelper(ing_min, ing_max,
1257  gp_p, itw,
1258  gfraw_in, 0,
1259  z_out, interp_order, zeropadding,
1260  verbosity);
1261 
1262  if (ing_max - ing_min >= 0)
1263  {
1264  Range r = joker;
1265 
1266  if (ing_max - ing_min + 1 != z_out.nelem())
1267  {
1268  r = Range(ing_min, ing_max-ing_min+1);
1269  }
1270 
1271  interp(gfraw_out.data(r, lat_index, lon_index),
1272  itw, gfraw_in.data(joker, lat_index, lon_index), gp_p);
1273  }
1274  }
1275  }
1276 }
1277 
1278 
1279 // Workspace method, doxygen header will be auto-generated.
1280 // 2007-07-25 Stefan Buehler
1281 void atm_fields_compactFromMatrix(// WS Output:
1282  GriddedField4& af, // atm_fields_compact
1283  // WS Input:
1284  const Index& atmosphere_dim,
1285  // WS Generic Input:
1286  const Matrix& im,
1287  // Control Parameters:
1288  const ArrayOfString& field_names,
1289  const Verbosity&)
1290 {
1291  if (1!=atmosphere_dim)
1292  {
1293  ostringstream os;
1294  os << "Atmospheric dimension must be one.";
1295  throw runtime_error( os.str() );
1296  }
1297 
1298  const Index np = im.nrows(); // Number of pressure levels.
1299  const Index nf = im.ncols()-1; // Total number of fields.
1300  ArrayOfIndex f_1; // indices of non-ignored fields
1301  // All fields called "ignore" will be ignored.
1302  String fn_upper; // Temporary variable to hold upper case field_names.
1303 
1304  if (field_names.nelem()!=nf)
1305  {
1306  ostringstream os;
1307  os << "Cannot extract fields from Matrix.\n"
1308  << "*field_names* must have one element less than there are\n"
1309  << "matrix columns.";
1310  throw runtime_error( os.str() );
1311  }
1312 
1313 
1314  // Remove additional fields from the field_names. All fields that
1315  // are flagged by 'ignore' in the field names, small or large letters,
1316  // are removed.
1317  for(Index f = 0; f < field_names.nelem(); f++)
1318  {
1319  fn_upper = field_names[f];
1320  fn_upper.toupper();
1321  //cout << "fieldname[" << f << "]: " << fn_upper;
1322  if(fn_upper != "IGNORE")
1323  {
1324  f_1.push_back(f);
1325  //cout << " put as element " << f_1.size()-1 << " into selection\n";
1326  }
1327  /*else
1328  {
1329  cout << " not selected\n";
1330  }*/
1331  }
1332 
1333  // Copy required field_names to a new variable called field_names_1
1334  Index nf_1=f_1.size(); // Number of required fields.
1335  ArrayOfString field_names_1(nf_1);
1336  for (Index f=0; f<nf_1; f++)
1337  {
1338  field_names_1[f] = field_names[f_1[f]];
1339  //cout << "new fieldname[" << f << "] (formerly [" << f_1[f] << "]): "
1340  // << field_names_1[f] << "\n";
1341  }
1342 
1343 
1344  // out3 << "Copying *" << im_name << "* to *atm_fields_compact*.\n";
1345 
1346  af.set_grid(GFIELD4_FIELD_NAMES, field_names_1);
1347 
1348  af.set_grid(GFIELD4_P_GRID, im(Range(joker),0));
1349 
1350  af.set_grid(GFIELD4_LAT_GRID, Vector());
1351  af.set_grid(GFIELD4_LON_GRID, Vector());
1352 
1353  af.resize(nf_1,np,1,1); // Resize it according to the required fields
1354  for (Index f = 0; f < nf_1; f++)
1355  af.data(f,Range(joker),0,0) = im(Range(joker),f_1[f]+1);
1356 }
1357 
1358 
1359 // Workspace method, doxygen header is auto-generated.
1360 // 2007-07-31 Stefan Buehler
1361 // 2011-05-04 Adapted by Gerrit Holl
1362 void atm_fields_compactAddConstant(// WS Output:
1363  GriddedField4& af,
1364  // Control Parameters:
1365  const String& name,
1366  const Numeric& value,
1367  const Verbosity& verbosity)
1368 {
1369  Index nf; // Will hold new size
1370 
1371  // Add book
1372  atm_fields_compactExpand(af, nf, name, verbosity);
1373 
1374  // Add the constant value:
1375  af.data( nf-1, Range(joker), Range(joker), Range(joker) ) = value;
1376 }
1377 
1378 // Workspace method, doxygen header is auto-generated
1379 // 2011-05-02 Gerrit Holl
1380 void atm_fields_compactAddSpecies(// WS Output:
1381  GriddedField4& atm_fields_compact,
1382  // WS Generic Input:
1383  const String& name,
1384  const GriddedField3& species,
1385  const Verbosity& verbosity)
1386 {
1387 
1388  assert(atm_fields_compact.checksize());
1389  assert(species.checksize());
1390 
1391  ConstVectorView af_p_grid = atm_fields_compact.get_numeric_grid(GFIELD4_P_GRID);
1392  ConstVectorView af_lat_grid = atm_fields_compact.get_numeric_grid(GFIELD4_LAT_GRID);
1393  ConstVectorView af_lon_grid = atm_fields_compact.get_numeric_grid(GFIELD4_LON_GRID);
1394  ConstVectorView sp_p_grid = species.get_numeric_grid(GFIELD3_P_GRID);
1395  ConstVectorView sp_lat_grid = species.get_numeric_grid(GFIELD3_LAT_GRID);
1396  ConstVectorView sp_lon_grid = species.get_numeric_grid(GFIELD3_LON_GRID);
1397 
1398  Index new_n_fields; // To be set in next line
1399  atm_fields_compactExpand(atm_fields_compact, new_n_fields, name, verbosity);
1400 
1401 
1402  // Interpolate species to atm_fields_compact
1403 
1404  // Common for all dim
1405  chk_interpolation_grids("species p_grid to atm_fields_compact p_grid",
1406  sp_p_grid,
1407  af_p_grid);
1408  ArrayOfGridPos p_gridpos(af_p_grid.nelem());
1409  // gridpos(p_gridpos, sp_p_grid, af_p_grid);
1410  p2gridpos(p_gridpos, sp_p_grid, af_p_grid);
1411 
1412  if (sp_lat_grid.nelem() > 1)
1413  {
1414  // Common for all dim>=2
1415  chk_interpolation_grids("species lat_grid to atm_fields_compact lat_grid",
1416  sp_lat_grid,
1417  af_lat_grid);
1418  ArrayOfGridPos lat_gridpos(af_lat_grid.nelem());
1419  gridpos(lat_gridpos, sp_lat_grid, af_lat_grid);
1420 
1421  if (sp_lon_grid.nelem() > 1)
1422  { // 3D-case
1423  chk_interpolation_grids("species lon_grid to atm_fields_compact lon_grid",
1424  sp_lon_grid,
1425  af_lon_grid);
1426  ArrayOfGridPos lon_gridpos(af_lon_grid.nelem());
1427  gridpos(lon_gridpos, sp_lon_grid, af_lon_grid);
1428 
1429  Tensor4 itw(p_gridpos.nelem(), lat_gridpos.nelem(), lon_gridpos.nelem(), 8);
1430  interpweights(itw, p_gridpos, lat_gridpos, lon_gridpos);
1431 
1432  Tensor3 newfield(af_p_grid.nelem(), af_lat_grid.nelem(), af_lon_grid.nelem());
1433  interp(newfield, itw, species.data, p_gridpos, lat_gridpos, lon_gridpos);
1434 
1435  atm_fields_compact.data(new_n_fields-1, joker, joker, joker) = newfield;
1436  } else { // 2D-case
1437 
1438  Tensor3 itw(p_gridpos.nelem(), lat_gridpos.nelem(), 4);
1439  interpweights(itw, p_gridpos, lat_gridpos);
1440 
1441  Matrix newfield(af_p_grid.nelem(), af_lat_grid.nelem());
1442  interp(newfield, itw, species.data(joker, joker, 0), p_gridpos, lat_gridpos);
1443 
1444  atm_fields_compact.data(new_n_fields-1, joker, joker, 0) = newfield;
1445  }
1446  } else { // 1D-case
1447  Matrix itw(p_gridpos.nelem(), 2);
1448  interpweights(itw, p_gridpos);
1449 
1450  Vector newfield(af_p_grid.nelem());
1451  interp(newfield, itw, species.data(joker, 0, 0), p_gridpos);
1452 
1453  atm_fields_compact.data(new_n_fields-1, joker, 0, 0) = newfield;
1454  }
1455 
1456 }
1457 
1458 
1459 
1460 // Workspace method, doxygen header is auto-generated
1461 // 2011-05-11 Gerrit Holl
1462 void batch_atm_fields_compactAddConstant(// WS Output:
1463  ArrayOfGriddedField4& batch_atm_fields_compact,
1464  // WS Generic Input:
1465  const String& name,
1466  const Numeric& value,
1467  const Verbosity& verbosity)
1468 {
1469  for (Index i=0; i<batch_atm_fields_compact.nelem(); i++)
1470  {
1471  atm_fields_compactAddConstant(batch_atm_fields_compact[i], name, value, verbosity);
1472  }
1473 
1474 }
1475 
1476 // Workspace method, doxygen header is auto-generated
1477 // 2011-05-09 Gerrit Holl
1478 void batch_atm_fields_compactAddSpecies(// WS Output:
1479  ArrayOfGriddedField4& batch_atm_fields_compact,
1480  // WS Generic Input:
1481  const String& name,
1482  const GriddedField3& species,
1483  const Verbosity& verbosity)
1484 {
1485  const Index nelem = batch_atm_fields_compact.nelem();
1486 
1487  String fail_msg;
1488  bool failed = false;
1489 
1490  // Parallelise this for-loop (some interpolation is being done, so it may
1491  // be beneficial)
1492 #pragma omp parallel for \
1493  if (!arts_omp_in_parallel() \
1494  && nelem >= arts_omp_get_max_threads())
1495  for (Index i=0; i<nelem; i++)
1496  {
1497  try {
1498  atm_fields_compactAddSpecies(batch_atm_fields_compact[i], name, species, verbosity);
1499  }
1500  catch (runtime_error e)
1501  {
1502 #pragma omp critical (batch_atm_fields_compactAddSpecies_fail)
1503  { fail_msg = e.what(); failed = true; }
1504  }
1505  }
1506 
1507  if (failed) throw runtime_error(fail_msg);
1508 }
1509 
1510 // Workspace method, doxygen header is auto-generated.
1511 void batch_atm_fields_compactFromArrayOfMatrix(// WS Output:
1512  ArrayOfGriddedField4& batch_atm_fields_compact,
1513  // WS Input:
1514  const Index& atmosphere_dim,
1515  // WS Generic Input:
1516  const ArrayOfMatrix& am,
1517  // Control Parameters:
1518  const ArrayOfString& field_names,
1519  const ArrayOfString& extra_field_names,
1520  const Vector& extra_field_values,
1521  const Verbosity& verbosity)
1522 {
1523  const Index amnelem = am.nelem();
1524 
1525  if (amnelem == 0)
1526  {
1527  ostringstream os;
1528  os << "No elements in atmospheric scenario batch.\n"
1529  << "Check, whether any batch atmosphere file has been read!";
1530  throw runtime_error( os.str() );
1531  }
1532 
1533  // We use the existing WSMs atm_fields_compactFromMatrix and
1534  // atm_fields_compactAddConstant to do most of the work.
1535 
1536  // Check that extra_field_names and extra_field_values have matching
1537  // dimensions:
1538  if (extra_field_names.nelem() != extra_field_values.nelem())
1539  {
1540  ostringstream os;
1541  os << "The keyword arguments extra_field_names and\n"
1542  << "extra_field_values must have matching dimensions.";
1543  throw runtime_error( os.str() );
1544  }
1545 
1546  // Make output variable the proper size:
1547  batch_atm_fields_compact.resize(amnelem);
1548 
1549  String fail_msg;
1550  bool failed = false;
1551 
1552  // Loop the batch cases:
1553 #pragma omp parallel for \
1554  if (!arts_omp_in_parallel() \
1555  && amnelem >= arts_omp_get_max_threads())
1556  for (Index i=0; i<amnelem; ++i)
1557  {
1558  // Skip remaining iterations if an error occurred
1559  if (failed) continue;
1560 
1561  // All the input variables are visible here, despite the
1562  // "default(none)". The reason is that they are return by
1563  // reference arguments of this function, which are shared
1564  // automatically.
1565 
1566  // The try block here is necessary to correctly handle
1567  // exceptions inside the parallel region.
1568  try
1569  {
1570  atm_fields_compactFromMatrix(batch_atm_fields_compact[i],
1571  atmosphere_dim,
1572  am[i],
1573  field_names,
1574  verbosity);
1575 
1576  for (Index j=0; j<extra_field_names.nelem(); ++j)
1577  atm_fields_compactAddConstant(batch_atm_fields_compact[i],
1578  extra_field_names[j],
1579  extra_field_values[j],
1580  verbosity);
1581  }
1582  catch (runtime_error e)
1583  {
1584 #pragma omp critical (batch_atm_fields_compactFromArrayOfMatrix_fail)
1585  { fail_msg = e.what(); failed = true; }
1586  }
1587  }
1588 
1589  if (failed) throw runtime_error(fail_msg);
1590 }
1591 
1592 
1593 void AtmFieldsFromCompact(// WS Output:
1594  Vector& p_grid,
1595  Vector& lat_grid,
1596  Vector& lon_grid,
1597  Tensor3& t_field,
1598  Tensor3& z_field,
1599  Tensor4& vmr_field,
1600  Tensor4& massdensity_field,
1601  // WS Input:
1602  const ArrayOfArrayOfSpeciesTag& abs_species,
1603  const ArrayOfString& part_species,
1604  const GriddedField4& atm_fields_compact,
1605  const Index& atmosphere_dim,
1606  const String& delim,
1607  const Verbosity&)
1608 {
1609  // Make a handle on atm_fields_compact to save typing:
1610  const GriddedField4& c = atm_fields_compact;
1611 
1612  // Check if the grids in our data match atmosphere_dim
1613  // (throws an error if the dimensionality is not correct):
1614  chk_atm_grids(atmosphere_dim,
1615  c.get_numeric_grid(GFIELD4_P_GRID),
1616  c.get_numeric_grid(GFIELD4_LAT_GRID),
1617  c.get_numeric_grid(GFIELD4_LON_GRID));
1618 
1619  const Index nf = c.get_grid_size(GFIELD4_FIELD_NAMES);
1620  const Index np = c.get_grid_size(GFIELD4_P_GRID);
1621  const Index nlat = c.get_grid_size(GFIELD4_LAT_GRID);
1622  const Index nlon = c.get_grid_size(GFIELD4_LON_GRID);
1623 
1624  // Grids:
1625  p_grid = c.get_numeric_grid(GFIELD4_P_GRID);
1626  lat_grid = c.get_numeric_grid(GFIELD4_LAT_GRID);
1627  lon_grid = c.get_numeric_grid(GFIELD4_LON_GRID);
1628 
1629  // The order of the fields is:
1630  // T[K] z[m] Particle_1[?] ... Particle_m[?] VMR_1[1] ... VMR_n[1]
1631 
1632  // Number of Particle&VMR species:
1633  const Index nsp = part_species.nelem();
1634  const Index nsa = nf-nsp-2;
1635 
1636  // Check that there is at least one VMR species:
1637  if (nsa<1)
1638  {
1639  ostringstream os;
1640  os << "There must be at least three fields in *atm_fields_compact*.\n"
1641  << "T, z, and at least one VMR.";
1642  throw runtime_error( os.str() );
1643  }
1644 
1645  // Check that first field is T:
1646  if (c.get_string_grid(GFIELD4_FIELD_NAMES)[0] != "T")
1647  {
1648  ostringstream os;
1649  os << "The first field must be \"T\", but it is:"
1650  << c.get_string_grid(GFIELD4_FIELD_NAMES)[0];
1651  throw runtime_error( os.str() );
1652  }
1653 
1654  // Check that second field is z:
1655  if (c.get_string_grid(GFIELD4_FIELD_NAMES)[1] != "z")
1656  {
1657  ostringstream os;
1658  os << "The second field must be \"z\"*, but it is:"
1659  << c.get_string_grid(GFIELD4_FIELD_NAMES)[1];
1660  throw runtime_error( os.str() );
1661  }
1662 
1663  // Check that the (supposed) particle fields match part_species:
1664  for (Index i=0; i<nsp; ++i)
1665  {
1666  const String tf_species = c.get_string_grid(GFIELD4_FIELD_NAMES)[2+i];
1667  String ps_species;
1668  parse_partfield_name(ps_species,part_species[i],delim);
1669  if (tf_species != ps_species)
1670  {
1671  ostringstream os;
1672  os << "Field name for particle field not valid: "
1673  << tf_species << "\n"
1674  << "Based on *part_species*, the field name should be: "
1675  << ps_species;
1676  throw runtime_error( os.str() );
1677  }
1678  }
1679 
1680  // Check that the (supposed) VMR fields match abs_species:
1681  for (Index i=0; i<nsa; ++i)
1682  {
1683  if (i >= abs_species.nelem())
1684  {
1685  ostringstream os;
1686  os << "Based on the field names for absorption species,\n"
1687  << "these species are missing in *abs_species*:\n";
1688  for (Index itf_species = i; itf_species < nsa; ++itf_species)
1689  os << c.get_string_grid(GFIELD4_FIELD_NAMES)[2+nsp+itf_species] << " ";
1690  throw runtime_error(os.str());
1691  }
1692 
1693  const String tf_species = c.get_string_grid(GFIELD4_FIELD_NAMES)[2+nsp+i];
1694 
1695  // Get name of species from abs_species:
1696  using global_data::species_data; // The species lookup data:
1697  const String as_species = species_data[abs_species[i][0].Species()].Name();
1698 
1699  // Species in field name and abs_species should be the same:
1700  if (tf_species != as_species)
1701  {
1702  ostringstream os;
1703  os << "Field name for absorption species field not valid: "
1704  << tf_species << "\n"
1705  << "Based on *abs_species*, the field name should be: "
1706  << as_species;
1707  throw runtime_error( os.str() );
1708  }
1709  }
1710 
1711  // Temperature field (first field):
1712  t_field.resize(np,nlat,nlon);
1713  t_field = c.data(0,Range(joker),Range(joker),Range(joker));
1714 
1715  // Altitude profile (second field):
1716  z_field.resize(np,nlat,nlon);
1717  z_field = c.data(1,Range(joker),Range(joker),Range(joker));
1718 
1719  //Particle profiles:
1720  massdensity_field.resize(nsp,np,nlat,nlon);
1721  massdensity_field = c.data(Range(2,nsp),Range(joker),Range(joker),Range(joker));
1722 
1723  // VMR profiles (remaining fields):
1724  vmr_field.resize(nsa,np,nlat,nlon);
1725  vmr_field = c.data(Range(2+nsp,nsa),Range(joker),Range(joker),Range(joker));
1726 }
1727 
1728 
1729 /* Workspace method: Doxygen documentation will be auto-generated */
1730 void AtmosphereSet1D(// WS Output:
1731  Index& atmosphere_dim,
1732  Vector& lat_grid,
1733  Vector& lon_grid,
1734  const Verbosity& verbosity)
1735 {
1736  CREATE_OUT2;
1737  CREATE_OUT3;
1738 
1739  out2 << " Sets the atmospheric dimensionality to 1.\n";
1740  out3 << " atmosphere_dim = 1\n";
1741  out3 << " lat_grid is set to be an empty vector\n";
1742  out3 << " lon_grid is set to be an empty vector\n";
1743 
1744  atmosphere_dim = 1;
1745  lat_grid.resize(0);
1746  lon_grid.resize(0);
1747 
1748 }
1749 
1750 
1751 /* Workspace method: Doxygen documentation will be auto-generated */
1752 void AtmosphereSet2D(// WS Output:
1753  Index& atmosphere_dim,
1754  Vector& lon_grid,
1755  const Verbosity& verbosity)
1756 {
1757  CREATE_OUT2;
1758  CREATE_OUT3;
1759 
1760  out2 << " Sets the atmospheric dimensionality to 2.\n";
1761  out3 << " atmosphere_dim = 2\n";
1762  out3 << " lon_grid is set to be an empty vector\n";
1763 
1764  atmosphere_dim = 2;
1765  lon_grid.resize(0);
1766 }
1767 
1768 
1769 /* Workspace method: Doxygen documentation will be auto-generated */
1770 void AtmosphereSet3D(// WS Output:
1771  Index& atmosphere_dim,
1772  const Verbosity& verbosity)
1773 {
1774  CREATE_OUT2;
1775  CREATE_OUT3;
1776 
1777  out2 << " Sets the atmospheric dimensionality to 3.\n";
1778  out3 << " atmosphere_dim = 3\n";
1779 
1780  atmosphere_dim = 3;
1781 }
1782 
1783 
1784 
1785 /* Workspace method: Doxygen documentation will be auto-generated */
1786 void AtmFieldsCalc(//WS Output:
1787  Tensor3& t_field,
1788  Tensor3& z_field,
1789  Tensor4& vmr_field,
1790  //WS Input
1791  const Vector& p_grid,
1792  const Vector& lat_grid,
1793  const Vector& lon_grid,
1794  const GriddedField3& t_field_raw,
1795  const GriddedField3& z_field_raw,
1796  const ArrayOfGriddedField3& vmr_field_raw,
1797  const Index& atmosphere_dim,
1798  // WS Generic Input:
1799  const Index& interp_order,
1800  const Index& vmr_zeropadding,
1801  const Index& vmr_nonegative,
1802  const Verbosity& verbosity)
1803 {
1804  CREATE_OUT2;
1805 
1806  const ConstVectorView tfr_p_grid = t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
1807  const ConstVectorView tfr_lat_grid = t_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
1808  const ConstVectorView tfr_lon_grid = t_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
1809  const ConstVectorView zfr_p_grid = z_field_raw.get_numeric_grid(GFIELD3_P_GRID);
1810  const ConstVectorView zfr_lat_grid = z_field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
1811  const ConstVectorView zfr_lon_grid = z_field_raw.get_numeric_grid(GFIELD3_LON_GRID);
1812 
1813  out2 << " Interpolation order: " << interp_order << "\n";
1814 
1815  // Basic checks of input variables
1816  //
1817  // Atmosphere
1818  chk_if_in_range( "atmosphere_dim", atmosphere_dim, 1, 3 );
1819  chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1820 
1821 
1822  //==========================================================================
1823  if ( atmosphere_dim == 1)
1824  {
1825  if( !( tfr_lat_grid.nelem() == 1 &&
1826  tfr_lon_grid.nelem() == 1 ))
1827  throw runtime_error(
1828  "Temperature data (T_field) has wrong dimension "
1829  "(2D or 3D).\n"
1830  );
1831 
1832  if( !( zfr_lat_grid.nelem() == 1 &&
1833  zfr_lon_grid.nelem() == 1 ))
1834  throw runtime_error(
1835  "Altitude data (z_field) has wrong dimension "
1836  "(2D or 3D).\n"
1837  );
1838 
1839  GriddedField3 temp_gfield3;
1840 
1841  GriddedFieldPRegrid(temp_gfield3, p_grid, t_field_raw, interp_order, 0, verbosity);
1842  t_field = temp_gfield3.data;
1843 
1844  GriddedFieldPRegrid(temp_gfield3, p_grid, z_field_raw, interp_order, 0, verbosity);
1845  z_field = temp_gfield3.data;
1846 
1847  ArrayOfGriddedField3 temp_agfield3;
1848  try {
1849  GriddedFieldPRegrid(temp_agfield3, p_grid, vmr_field_raw, interp_order,
1850  vmr_zeropadding, verbosity);
1851  } catch (runtime_error e) {
1852  ostringstream os;
1853  os << e.what() << "\n"
1854  << "Note that you can explicitly set vmr_zeropadding to 1 in the method call.";
1855  throw runtime_error(os.str());
1856 
1857  }
1858  FieldFromGriddedField(vmr_field, p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
1859  }
1860 
1861  //=========================================================================
1862  else if(atmosphere_dim == 2)
1863  {
1864  if( tfr_lat_grid.nelem() == 1 &&
1865  tfr_lon_grid.nelem() == 1 )
1866  throw runtime_error(
1867  "Raw data has wrong dimension (1D). "
1868  "You have to use \n"
1869  "AtmFieldsCalcExpand1D instead of AtmFieldsCalc."
1870  );
1871 
1872  //Resize variables
1873  t_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
1874  z_field.resize(p_grid.nelem(), lat_grid.nelem(), 1);
1875  vmr_field.resize(vmr_field_raw.nelem(), p_grid.nelem(), lat_grid.nelem(),
1876  1);
1877 
1878 
1879  // Gridpositions:
1880  ArrayOfGridPosPoly gp_p(p_grid.nelem());
1881  ArrayOfGridPosPoly gp_lat(lat_grid.nelem());
1882 
1883  // Interpolate t_field:
1884 
1885  // Check that interpolation grids are ok (and throw a detailed
1886  // error message if not):
1887  chk_interpolation_pgrids("Raw temperature to p_grid, 2D case",
1888  tfr_p_grid,
1889  p_grid,
1890  interp_order);
1891  chk_interpolation_grids("Raw temperature to lat_grid, 2D case",
1892  tfr_lat_grid,
1893  lat_grid,
1894  interp_order);
1895 
1896  // Calculate grid positions:
1897  p2gridpos_poly( gp_p, tfr_p_grid, p_grid, interp_order );
1898  gridpos_poly( gp_lat, tfr_lat_grid, lat_grid, interp_order );
1899 
1900  // Interpolation weights:
1901  Tensor3 itw(p_grid.nelem(), lat_grid.nelem(), (interp_order+1)*(interp_order+1));
1902  // (4 interpolation weights are required for example for linear 2D interpolation)
1903  interpweights( itw, gp_p, gp_lat);
1904 
1905  // Interpolate:
1906  interp( t_field(joker, joker, 0 ), itw,
1907  t_field_raw.data(joker, joker, 0), gp_p, gp_lat);
1908 
1909 
1910  // Interpolate z_field:
1911 
1912  // Check that interpolation grids are ok (and throw a detailed
1913  // error message if not):
1914  chk_interpolation_pgrids("Raw z to p_grid, 2D case",
1915  zfr_p_grid,
1916  p_grid,
1917  interp_order);
1918  chk_interpolation_grids("Raw z to lat_grid, 2D case",
1919  zfr_lat_grid,
1920  lat_grid,
1921  interp_order);
1922 
1923  // Calculate grid positions:
1924  p2gridpos_poly( gp_p, zfr_p_grid, p_grid, interp_order );
1925  gridpos_poly( gp_lat, zfr_lat_grid, lat_grid, interp_order );
1926 
1927  // Interpolation weights:
1928  interpweights( itw, gp_p, gp_lat);
1929 
1930  // Interpolate:
1931  interp( z_field(joker, joker, 0), itw,
1932  z_field_raw.data(joker, joker, 0), gp_p, gp_lat);
1933 
1934 
1935  // Interpolate vmr_field.
1936  // Loop over the gaseous species:
1937  for (Index gas_i = 0; gas_i < vmr_field_raw.nelem(); gas_i++)
1938  {
1939  ostringstream os;
1940 
1941  if( !( vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID).nelem() != 1 &&
1942  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LON_GRID).nelem() == 1 ))
1943  {
1944  os << "VMR data of the " << gas_i << "th species has "
1945  << "wrong dimension (1D or 3D). \n";
1946  throw runtime_error( os.str() );
1947  }
1948 
1949  // Check that interpolation grids are ok (and throw a detailed
1950  // error message if not):
1951  os << "Raw VMR[" << gas_i << "] to p_grid, 2D case";
1952  chk_interpolation_pgrids(os.str(),
1953  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_P_GRID),
1954  p_grid,
1955  interp_order);
1956  os.str("");
1957  os << "Raw VMR[" << gas_i << "] to lat_grid, 2D case";
1958  chk_interpolation_grids(os.str(),
1959  vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID),
1960  lat_grid,
1961  interp_order);
1962 
1963  // Calculate grid positions:
1964  p2gridpos_poly(gp_p, vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_P_GRID),
1965  p_grid, interp_order);
1966  gridpos_poly(gp_lat, vmr_field_raw[gas_i].get_numeric_grid(GFIELD3_LAT_GRID),
1967  lat_grid, interp_order);
1968 
1969  // Interpolation weights:
1970  interpweights( itw, gp_p, gp_lat);
1971 
1972  // Interpolate:
1973  interp( vmr_field(gas_i, joker, joker, 0),
1974  itw, vmr_field_raw[gas_i].data(joker, joker, 0),
1975  gp_p, gp_lat);
1976  }
1977  }
1978 
1979  //================================================================
1980  // atmosphere_dim = 3
1981  else if(atmosphere_dim == 3)
1982  {
1983  if( tfr_lat_grid.nelem() == 1 &&
1984  tfr_lon_grid.nelem() == 1 )
1985  throw runtime_error(
1986  "Raw data has wrong dimension. You have to use \n"
1987  "AtmFieldsCalcExpand1D instead of AtmFieldsCalc."
1988  );
1989 
1990  GriddedField3 temp_gfield3;
1991 
1992  GriddedFieldLatLonRegrid(temp_gfield3, lat_grid, lon_grid, t_field_raw, interp_order, verbosity);
1993  GriddedFieldPRegrid(temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
1994  t_field = temp_gfield3.data;
1995 
1996  GriddedFieldLatLonRegrid(temp_gfield3, lat_grid, lon_grid, z_field_raw, interp_order, verbosity);
1997  GriddedFieldPRegrid(temp_gfield3, p_grid, temp_gfield3, interp_order, 0, verbosity);
1998  z_field = temp_gfield3.data;
1999 
2000  ArrayOfGriddedField3 temp_agfield3;
2001  GriddedFieldLatLonRegrid(temp_agfield3, lat_grid, lon_grid, vmr_field_raw, interp_order, verbosity);
2002  try {
2003  GriddedFieldPRegrid(temp_agfield3, p_grid, temp_agfield3, interp_order,
2004  vmr_zeropadding, verbosity);
2005  } catch (runtime_error e) {
2006  ostringstream os;
2007  os << e.what() << "\n"
2008  << "Note that you can explicitly set vmr_zeropadding to 1 in the method call.";
2009  throw runtime_error(os.str());
2010 
2011  }
2012  FieldFromGriddedField(vmr_field, p_grid, lat_grid, lon_grid, temp_agfield3, verbosity);
2013  }
2014  else
2015  {
2016  // We can never get here, since there was a runtime
2017  // error check for atmosphere_dim at the beginning.
2018  assert(false);
2019  }
2020 
2021  // remove negatives?
2022  if( vmr_nonegative )
2023  {
2024  for( Index ib=0; ib<vmr_field.nbooks(); ib++ )
2025  {
2026  for( Index ip=0; ip<vmr_field.npages(); ip++ )
2027  {
2028  for( Index ir=0; ir<vmr_field.nrows(); ir++ )
2029  {
2030  for( Index ic=0; ic<vmr_field.ncols(); ic++ )
2031  {
2032  if( vmr_field(ib,ip,ir,ic) < 0 )
2033  { vmr_field(ib,ip,ir,ic) = 0; }
2034  } } } }
2035  }
2036 }
2037 
2038 
2039 /* Workspace method: Doxygen documentation will be auto-generated */
2040 void AtmFieldsCalcExpand1D(Tensor3& t_field,
2041  Tensor3& z_field,
2042  Tensor4& vmr_field,
2043  const Vector& p_grid,
2044  const Vector& lat_grid,
2045  const Vector& lon_grid,
2046  const GriddedField3& t_field_raw,
2047  const GriddedField3& z_field_raw,
2048  const ArrayOfGriddedField3& vmr_field_raw,
2049  const Index& atmosphere_dim,
2050  const Index& interp_order,
2051  const Index& vmr_zeropadding,
2052  const Index& vmr_nonegative,
2053  const Verbosity& verbosity)
2054 {
2055  chk_if_in_range( "atmosphere_dim", atmosphere_dim, 1, 3 );
2056  chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
2057 
2058  if( atmosphere_dim == 1 )
2059  { throw runtime_error(
2060  "This function is intended for 2D and 3D. For 1D, use *AtmFieldsCalc*.");}
2061 
2062  // Make 1D interpolation using some dummy variables
2063  Vector vempty(0);
2064  Tensor3 t_temp, z_temp;
2065  Tensor4 vmr_temp;
2066  AtmFieldsCalc(t_temp, z_temp, vmr_temp, p_grid, vempty, vempty,
2067  t_field_raw, z_field_raw, vmr_field_raw, 1, interp_order,
2068  vmr_zeropadding, vmr_nonegative, verbosity);
2069 
2070  // Move values from the temporary tensors to the return arguments
2071  const Index np = p_grid.nelem();
2072  const Index nlat = lat_grid.nelem();
2073  Index nlon = lon_grid.nelem();
2074  if( atmosphere_dim == 2 )
2075  { nlon = 1; }
2076  const Index nspecies = vmr_temp.nbooks();
2077  //
2078  assert( t_temp.npages() == np );
2079  //
2080  t_field.resize( np, nlat, nlon );
2081  z_field.resize( np, nlat, nlon );
2082  vmr_field.resize( nspecies, np, nlat, nlon );
2083  //
2084  for( Index ilon=0; ilon<nlon; ilon++ )
2085  {
2086  for( Index ilat=0; ilat<nlat; ilat++ )
2087  {
2088  for( Index ip=0; ip<np; ip++ )
2089  {
2090  t_field(ip,ilat,ilon) = t_temp(ip,0,0);
2091  z_field(ip,ilat,ilon) = z_temp(ip,0,0);
2092  for( Index is=0; is<nspecies; is++ )
2093  { vmr_field(is,ip,ilat,ilon) = vmr_temp(is,ip,0,0); }
2094  }
2095  }
2096  }
2097 }
2098 
2099 
2100 /* Workspace method: Doxygen documentation will be auto-generated */
2101 void AtmFieldsExpand1D(Tensor3& t_field,
2102  Tensor3& z_field,
2103  Tensor4& vmr_field,
2104  const Vector& p_grid,
2105  const Vector& lat_grid,
2106  const Vector& lon_grid,
2107  const Index& atmosphere_dim,
2108  const Verbosity&)
2109 {
2110  chk_if_in_range( "atmosphere_dim", atmosphere_dim, 1, 3 );
2111  chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
2112 
2113  // Sizes
2114  const Index np = p_grid.nelem();
2115  const Index nlat = lat_grid.nelem();
2116  const Index nlon = max( Index(1), lon_grid.nelem() );
2117  const Index nspecies = vmr_field.nbooks();
2118 
2119  if( atmosphere_dim == 1 )
2120  { throw runtime_error( "No use in calling this method for 1D.");}
2121  chk_atm_field( "t_field", t_field, 1, p_grid, Vector(0), Vector(0) );
2122  chk_atm_field( "z_field", z_field, 1, p_grid, Vector(0), Vector(0) );
2123  if( nspecies )
2124  chk_atm_field( "vmr_field", vmr_field, 1, nspecies, p_grid, Vector(0),
2125  Vector(0) );
2126 
2127  // Temporary containers
2128  Tensor3 t_temp = t_field, z_temp = z_field;
2129  Tensor4 vmr_temp = vmr_field;
2130 
2131  // Resize and fill
2132  t_field.resize( np, nlat, nlon );
2133  z_field.resize( np, nlat, nlon );
2134  vmr_field.resize( nspecies, np, nlat, nlon );
2135  //
2136  for( Index ilon=0; ilon<nlon; ilon++ )
2137  {
2138  for( Index ilat=0; ilat<nlat; ilat++ )
2139  {
2140  for( Index ip=0; ip<np; ip++ )
2141  {
2142  t_field(ip,ilat,ilon) = t_temp(ip,0,0);
2143  z_field(ip,ilat,ilon) = z_temp(ip,0,0);
2144  for( Index is=0; is<nspecies; is++ )
2145  { vmr_field(is,ip,ilat,ilon) = vmr_temp(is,ip,0,0); }
2146  }
2147  }
2148  }
2149 }
2150 
2151 
2152 
2153 /* Workspace method: Doxygen documentation will be auto-generated */
2154 void AtmFieldsRefinePgrid(// WS Output:
2155  Vector& p_grid,
2156  Tensor3& t_field,
2157  Tensor3& z_field,
2158  Tensor4& vmr_field,
2159  // WS Input:
2160  const Vector& lat_grid,
2161  const Vector& lon_grid,
2162  const Index& atmosphere_dim,
2163  // Control Parameters:
2164  const Numeric& p_step,
2165  const Verbosity&)
2166 {
2167  // Checks on input parameters:
2168 
2169  // We don't actually need lat_grid and lon_grid, but we have them as
2170  // input variables, so that we can use the standard functions to
2171  // check atmospheric fields and grids. A bit cheesy, but I don't
2172  // want to program all the checks explicitly.
2173 
2174  // Check grids:
2175  chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
2176 
2177  // Check T field:
2178  chk_atm_field("t_field", t_field, atmosphere_dim,
2179  p_grid, lat_grid, lon_grid);
2180 
2181  // Check z field:
2182  chk_atm_field("z_field", z_field, atmosphere_dim,
2183  p_grid, lat_grid, lon_grid);
2184 
2185  // Check VMR field (and abs_species):
2186  chk_atm_field("vmr_field", vmr_field, atmosphere_dim,
2187  vmr_field.nbooks(), p_grid, lat_grid, lon_grid);
2188 
2189  // Check the keyword arguments:
2190  if ( p_step <= 0 )
2191  {
2192  ostringstream os;
2193  os << "The keyword argument p_step must be >0.";
2194  throw runtime_error( os.str() );
2195  }
2196 
2197  // Ok, all input parameters seem to be reasonable.
2198 
2199  // We will need the log of the pressure grid:
2200  Vector log_p_grid(p_grid.nelem());
2201  transform(log_p_grid, log, p_grid);
2202 
2203  // const Numeric epsilon = 0.01 * p_step; // This is the epsilon that
2204  // // we use for comparing p grid spacings.
2205 
2206  // Construct abs_p
2207  // ---------------
2208 
2209  ArrayOfNumeric log_abs_p_a; // We take log_abs_p_a as an array of
2210  // Numeric, so that we can easily
2211  // build it up by appending new elements to the end.
2212 
2213  // Check whether there are pressure levels that are further apart
2214  // (in log(p)) than p_step, and insert additional levels if
2215  // necessary:
2216 
2217  log_abs_p_a.push_back(log_p_grid[0]);
2218 
2219  for (Index i=1; i<log_p_grid.nelem(); ++i)
2220  {
2221  const Numeric dp = log_p_grid[i-1] - log_p_grid[i]; // The grid is descending.
2222 
2223  const Numeric dp_by_p_step = dp/p_step;
2224  // cout << "dp_by_p_step: " << dp_by_p_step << "\n";
2225 
2226  // How many times does p_step fit into dp?
2227  const Index n = (Index) ceil(dp_by_p_step);
2228  // n is the number of intervals that we want to have in the
2229  // new grid. The number of additional points to insert is
2230  // n-1. But we have to insert the original point as well.
2231  // cout << n << "\n";
2232 
2233  const Numeric ddp = dp/(Numeric)n;
2234  // cout << "ddp: " << ddp << "\n";
2235 
2236  for (Index j=1; j<=n; ++j)
2237  log_abs_p_a.push_back(log_p_grid[i-1] - (Numeric)j*ddp);
2238  }
2239 
2240  // Copy to a proper vector, we need this also later for
2241  // interpolation:
2242  Vector log_abs_p(log_abs_p_a.nelem());
2243  for (Index i=0; i<log_abs_p_a.nelem(); ++i)
2244  log_abs_p[i] = log_abs_p_a[i];
2245 
2246  // Copy the new grid to abs_p, removing the log:
2247  Vector abs_p(log_abs_p.nelem());
2248  transform(abs_p, exp, log_abs_p);
2249 
2250 
2251  // We will also have to interpolate T and VMR profiles to the new
2252  // pressure grid. We interpolate in log(p), as usual in ARTS.
2253 
2254  // Grid positions:
2255  ArrayOfGridPos gp(log_abs_p.nelem());
2256  gridpos(gp, log_p_grid, log_abs_p);
2257 
2258  // Interpolation weights:
2259  Matrix itw(gp.nelem(),2);
2260  interpweights(itw,gp);
2261 
2262  // Extent of latitude and longitude grids:
2263  Index nlat = lat_grid.nelem();
2264  Index nlon = lon_grid.nelem();
2265 
2266  // This here is needed so that the method works for 1D and 2D
2267  // atmospheres as well, not just for 3D:
2268  if (0 == nlat) nlat = 1;
2269  if (0 == nlon) nlon = 1;
2270 
2271  // Define variables for output fields:
2272  Tensor3 abs_t(log_abs_p.nelem(), nlat, nlon);
2273  Tensor3 abs_z(log_abs_p.nelem(), nlat, nlon);
2274  Tensor4 abs_vmr(vmr_field.nbooks(),
2275  log_abs_p.nelem(), nlat, nlon);
2276 
2277  for (Index ilat=0; ilat<nlat; ++ilat)
2278  for (Index ilon=0; ilon<nlon; ++ilon)
2279  {
2280  interp(abs_t(joker, ilat, ilon),
2281  itw,
2282  t_field(joker, ilat, ilon),
2283  gp);
2284 
2285  interp(abs_z(joker, ilat, ilon),
2286  itw,
2287  z_field(joker, ilat, ilon),
2288  gp);
2289 
2290  for (Index ivmr=0; ivmr<vmr_field.nbooks(); ++ivmr)
2291  interp(abs_vmr(ivmr, joker, ilat, ilon),
2292  itw,
2293  vmr_field(ivmr, joker, ilat, ilon),
2294  gp);
2295  }
2296 
2297 
2298  // Copy back the new fields to p_grid, t_field, z_field, vmr_field:
2299  p_grid = abs_p;
2300  t_field = abs_t;
2301  z_field = abs_z;
2302  vmr_field = abs_vmr;
2303 }
2304 
2305 
2306 
2307 /* Workspace method: Doxygen documentation will be auto-generated */
2308 void AtmRawRead(//WS Output:
2309  GriddedField3& t_field_raw,
2310  GriddedField3& z_field_raw,
2311  ArrayOfGriddedField3& vmr_field_raw,
2312  //WS Input:
2313  const ArrayOfArrayOfSpeciesTag& abs_species,
2314  //Keyword:
2315  const String& basename,
2316  const Verbosity& verbosity)
2317 {
2318  CREATE_OUT3;
2319 
2320  String tmp_basename = basename;
2321  if (basename.length() && basename[basename.length()-1] != '/')
2322  tmp_basename += ".";
2323 
2324  // Read the temperature field:
2325  String file_name = tmp_basename + "t.xml";
2326  xml_read_from_file( file_name, t_field_raw, verbosity);
2327 
2328  out3 << "Temperature field read from file: " << file_name << "\n";
2329 
2330  // Read geometrical altitude field:
2331  file_name = tmp_basename + "z.xml";
2332  xml_read_from_file( file_name, z_field_raw, verbosity);
2333 
2334  out3 << "Altitude field read from file: " << file_name << "\n";
2335 
2336  // Clear out vmr_field_raw
2337  vmr_field_raw.resize(0);
2338 
2339  // The species lookup data:
2340  using global_data::species_data;
2341 
2342  // We need to read one profile for each tag group.
2343  for ( Index i=0; i<abs_species.nelem(); i ++)
2344  {
2345  // Determine the name.
2346  file_name = tmp_basename +
2347  species_data[abs_species[i][0].Species()].Name() + ".xml";
2348 
2349  // Add an element for this tag group to the vmr profiles:
2350  GriddedField3 vmr_field_data;
2351  vmr_field_raw.push_back(vmr_field_data);
2352 
2353  // Read the VMR:
2354  xml_read_from_file( file_name, vmr_field_raw[vmr_field_raw.nelem()-1], verbosity);
2355 
2356  // state the source of profile.
2357  out3 << " " << species_data[abs_species[i][0].Species()].Name()
2358  << " profile read from file: " << file_name << "\n";
2359  }
2360 }
2361 
2362 
2363 
2364 /* Workspace method: Doxygen documentation will be auto-generated */
2365 void InterpAtmFieldToPosition(
2366  Numeric& outvalue,
2367  const Index& atmosphere_dim,
2368  const Vector& p_grid,
2369  const Vector& lat_grid,
2370  const Vector& lon_grid,
2371  const Tensor3& z_field,
2372  const Vector& rtp_pos,
2373  const Tensor3& field,
2374  const Verbosity& verbosity)
2375 {
2376  // Input checks
2377  chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
2378  chk_atm_field( "input argument *field*", field, atmosphere_dim,
2379  p_grid, lat_grid, lon_grid );
2380  chk_rte_pos( atmosphere_dim, rtp_pos );
2381 
2382  // Determine grid positions
2383  GridPos gp_p, gp_lat, gp_lon;
2384  rte_pos2gridpos( gp_p, gp_lat, gp_lon, atmosphere_dim,
2385  p_grid, lat_grid, lon_grid, z_field, rtp_pos );
2386 
2387  // Interpolate
2388  outvalue = interp_atmfield_by_gp( atmosphere_dim, field, gp_p, gp_lat,
2389  gp_lon );
2390 
2391  CREATE_OUT3;
2392  out3 << " Result = " << outvalue << "\n";
2393 }
2394 
2395 
2396 
2397 /* Workspace method: Doxygen documentation will be auto-generated */
2398 void p_gridDensify(
2399  Vector& p_grid,
2400  const Index& nfill,
2401  const Verbosity& verbosity )
2402 {
2403  if( nfill < 0 )
2404  { throw runtime_error( "Argument *nfill* must be >= 0." ); }
2405 
2406  // Nothing to do if nfill=0
2407  if( nfill > 0 )
2408  {
2409  // Make copy of _pgrid and allocate new size
2410  const Vector p0( p_grid);
2411  const Index n0 = p0.nelem();
2412  //
2413  p_grid.resize( (n0-1)*(1+nfill) + 1 );
2414 
2415  Index iout = 0;
2416  p_grid[0] = p0[0];
2417 
2418  for( Index i=1; i<n0; i++ )
2419  {
2420  Vector pnew;
2421  VectorNLogSpace( pnew, 2+nfill, p0[i-1], p0[i], verbosity );
2422  for( Index j=1; j<nfill+2; j++ )
2423  {
2424  iout += 1;
2425  p_grid[iout] = pnew[j];
2426  }
2427  }
2428  }
2429 }
2430 
2431 
2432 
2433 /* Workspace method: Doxygen documentation will be auto-generated */
2434 void p_gridFromZRaw(//WS Output
2435  Vector& p_grid,
2436  //WS Input
2437  const GriddedField3& z_field_raw,
2438  const Index& no_negZ,
2439  const Verbosity&)
2440 {
2441  // original version excludes negative z. not clear, why this is. maybe is
2442  // currently a convention somehwere in ARTS (DOIT?). negative z seem, however,
2443  // to work fine for clear-sky cases. so we make the negative z exclude an
2444  // option (for consistency until unclarities solved, default: do exclude)
2445  Vector p_grid_raw=z_field_raw.get_numeric_grid(GFIELD3_P_GRID);
2446 
2447  Index i;
2448  if (is_increasing(z_field_raw.data(joker,0,0)))
2449  {
2450  i=0;
2451  if (no_negZ)
2452  {
2453  while ( z_field_raw.data(i,0,0)< 0.0 ) i++;
2454  }
2455  p_grid=p_grid_raw[Range(i,joker)];
2456  }
2457  else if (is_decreasing(z_field_raw.data(joker,0,0)))
2458  {
2459  i=z_field_raw.data.npages()-1;
2460  if (no_negZ)
2461  {
2462  while ( z_field_raw.data(i,0,0)< 0.0 ) i--;
2463  }
2464  p_grid=p_grid_raw[Range(i,joker,-1)];
2465  }
2466  else
2467  {
2468  ostringstream os;
2469  os << "z_field_raw needs to be monotonous, but this is not the case.\n";
2470  throw runtime_error( os.str() );
2471  }
2472 }
2473 
2474 
2475 /* Workspace method: Doxygen documentation will be auto-generated */
2476 void lat_gridFromRawField(//WS Output
2477  Vector& lat_grid,
2478  //WS Input
2479  const GriddedField3& field_raw,
2480  const Verbosity&)
2481 {
2482  lat_grid = field_raw.get_numeric_grid(GFIELD3_LAT_GRID);
2483 }
2484 
2485 
2486 /* Workspace method: Doxygen documentation will be auto-generated */
2487 void lon_gridFromRawField(//WS Output
2488  Vector& lon_grid,
2489  //WS Input
2490  const GriddedField3& field_raw,
2491  const Verbosity&)
2492 {
2493  lon_grid = field_raw.get_numeric_grid(GFIELD3_LON_GRID);
2494 }
2495 
2496 
2497 /* Workspace method: Doxygen documentation will be auto-generated */
2498 void wind_u_fieldIncludePlanetRotation(
2499  Tensor3& wind_u_field,
2500  const Index& atmosphere_dim,
2501  const Vector& p_grid,
2502  const Vector& lat_grid,
2503  const Vector& lon_grid,
2504  const Vector& refellipsoid,
2505  const Tensor3& z_field,
2506  const Numeric& planet_rotation_period,
2507  const Verbosity&)
2508 {
2509  if( atmosphere_dim < 3 )
2510  throw runtime_error( "No need to use this method for 1D and 2D." );
2511 
2512  const Index np = p_grid.nelem();
2513  const Index na = lat_grid.nelem();
2514  const Index no = lon_grid.nelem();
2515 
2516  chk_atm_field( "z_field", z_field, atmosphere_dim,
2517  p_grid, lat_grid, lon_grid );
2518  if( wind_u_field.npages() > 0 )
2519  { chk_atm_field( "wind_u_field", wind_u_field, atmosphere_dim,
2520  p_grid, lat_grid, lon_grid );}
2521  else
2522  {
2523  wind_u_field.resize( np, na, no );
2524  wind_u_field = 0.;
2525  }
2526 
2527  const Numeric k1 = 2 * PI / planet_rotation_period;
2528 
2529 
2530  for( Index a=0; a<na; a++ )
2531  {
2532  const Numeric k2 = k1 * cos( DEG2RAD*lat_grid[a] );
2533  const Numeric re = refell2r( refellipsoid, lat_grid[a] );
2534 
2535  for( Index o=0; o<no; o++ )
2536  {
2537  for( Index p=0; p<np; p++ )
2538  {
2539  wind_u_field(p,a,o) += k2 * ( re + z_field(p,a,o) );
2540  }
2541  }
2542  }
2543 }
2544 
2545 
2546 
2547 
2548 // A small help function
2549 void z2g(
2550  Numeric& g,
2551  const Numeric& r,
2552  const Numeric& g0,
2553  const Numeric& z )
2554 {
2555  const Numeric x = r/(r+z);
2556  g = g0 * x*x;
2557 }
2558 
2559 /* Workspace method: Doxygen documentation will be auto-generated */
2560 void z_fieldFromHSE(
2561  Workspace& ws,
2562  Tensor3& z_field,
2563  const Index& atmosphere_dim,
2564  const Vector& p_grid,
2565  const Vector& lat_grid,
2566  const Vector& lon_grid,
2567  const Vector& lat_true,
2568  const Vector& lon_true,
2569  const ArrayOfArrayOfSpeciesTag& abs_species,
2570  const Tensor3& t_field,
2571  const Tensor4& vmr_field,
2572  const Vector& refellipsoid,
2573  const Matrix& z_surface,
2574  const Index& atmfields_checked,
2575  const Agenda& g0_agenda,
2576  const Numeric& molarmass_dry_air,
2577  const Numeric& p_hse,
2578  const Numeric& z_hse_accuracy,
2579  const Verbosity& verbosity)
2580 {
2581  if( atmfields_checked != 1 )
2582  throw runtime_error( "The atmospheric fields must be flagged to have "
2583  "passed a consistency check (atmfields_checked=1)." );
2584 
2585  // Some general variables
2586  const Index np = p_grid.nelem();
2587  const Index nlat = t_field.nrows();
2588  const Index nlon = t_field.ncols();
2589  //
2590  const Index firstH2O = find_first_species_tg( abs_species,
2591  species_index_from_species_name("H2O") );
2592 
2593  if( firstH2O < 0 )
2594  {
2595  ostringstream os;
2596  os << "No water vapour tag group in *abs_species*.\n"
2597  << "Be aware that this leads to significant variations in atmospheres\n"
2598  << "that contain considerable amounts of water vapour (e.g. Earth)!\n";
2599  CREATE_OUT1;
2600  out1 << os.str();
2601  //throw runtime_error(
2602  // "Water vapour is required (must be a tag group in *abs_species*)." );
2603  }
2604  //
2605  if( p_hse>p_grid[0] || p_hse < p_grid[np-1] )
2606  {
2607  ostringstream os;
2608  os << "The value of *p_hse* must be inside the range of *p_grid*:"
2609  << " p_hse = " << p_hse << " Pa\n"
2610  << " p_grid = << p_grid[np-1]" << " - " << p_grid[0] << " Pa\n";
2611  throw runtime_error( os.str() );
2612  }
2613  //
2614  if( z_hse_accuracy <= 0 )
2615  { throw runtime_error( "The value of *z_hse_accuracy* must be > 0." ); }
2616  //
2617  chk_latlon_true( atmosphere_dim, lat_grid, lat_true, lon_true );
2618 
2619 
2620  // Determine interpolation weights for p_hse
2621  //
2622  ArrayOfGridPos gp(1);
2623  Matrix itw( 1, 2);
2624  p2gridpos( gp, p_grid, Vector(1,p_hse) );
2625  interpweights ( itw, gp );
2626 
2627 
2628  // // Molecular weight of water vapour
2629  const Numeric mw = 18.016;
2630 
2631  // mw/molarmass_dry_air matches eps in Eq. 3.14 in Wallace&Hobbs:
2632  const Numeric k = 1 - mw/molarmass_dry_air;
2633 
2634  // Gas constant for 1kg dry air:
2635  const Numeric rd = 1e3*GAS_CONSTANT/molarmass_dry_air;
2636 
2637  // The calculations
2638  //
2639  for( Index ilat=0; ilat<nlat; ilat++ )
2640  {
2641  // The reference ellipsoid is already adjusted to internal 1D or 2D
2642  // views, and lat_grid is the relevant grid for *refellipsoid*, also
2643  // for 2D. On the other hand, extraction of g0 requires that the true
2644  // position is determined.
2645 
2646  // Radius of reference ellipsoid
2647  Numeric re;
2648  if( atmosphere_dim == 1 )
2649  { re = refellipsoid[0]; }
2650  else
2651  { re = refell2r( refellipsoid, lat_grid[ilat] ); }
2652 
2653  for( Index ilon=0; ilon<nlon; ilon++ )
2654  {
2655  // Determine true latitude and longitude
2656  Numeric lat, lon;
2657  Vector pos(atmosphere_dim); // pos[0] can be a dummy value
2658  if( atmosphere_dim >= 2 )
2659  {
2660  pos[1] = lat_grid[ilat];
2661  if( atmosphere_dim == 3 )
2662  { pos[2] = lon_grid[ilon]; }
2663  }
2664  pos2true_latlon( lat, lon, atmosphere_dim, lat_grid, lat_true,
2665  lon_true, pos );
2666 
2667  // Get g0
2668  Numeric g0;
2669  g0_agendaExecute( ws, g0, lat, lon, g0_agenda );
2670 
2671  // Determine altitude for p_hse
2672  Vector z_hse(1);
2673  interp( z_hse, itw, z_field(joker,ilat,ilon), gp );
2674 
2675  Numeric z_acc = 2 * z_hse_accuracy;
2676 
2677  while( z_acc > z_hse_accuracy )
2678  {
2679  z_acc = 0;
2680  Numeric g2=g0;
2681 
2682  for( Index ip=0; ip<np-1; ip++ )
2683  {
2684  // Calculate average g
2685  if( ip == 0 )
2686  { z2g( g2, re, g0, z_field(ip,ilat,ilon) ); }
2687  const Numeric g1 = g2;
2688  z2g( g2, re, g0, z_field(ip+1,ilat,ilon) );
2689  //
2690  const Numeric g = ( g1 + g2 ) / 2.0;
2691 
2692  //Average water vapour VMR
2693  Numeric hm;
2694  if( firstH2O < 0 )
2695  {
2696  hm = 0.0;
2697  }
2698  else
2699  {
2700  hm = 0.5* ( vmr_field(firstH2O,ip,ilat,ilon)
2701  + vmr_field(firstH2O,ip+1,ilat,ilon) );
2702  }
2703 
2704  // Average virtual temperature (no liquid water)
2705  // (cf. e.g. Eq. 3.16 in Wallace&Hobbs)
2706  const Numeric tv = ( 1 / ( 2 * (1-hm*k) ) ) *
2707  ( t_field(ip,ilat,ilon) + t_field(ip+1,ilat,ilon) );
2708 
2709  // Change in vertical altitude from ip to ip+1
2710  // (cf. e.g. Eq. 3.24 in Wallace&Hobbs)
2711  const Numeric dz = rd*(tv/g)*log( p_grid[ip]/p_grid[ip+1] );
2712 
2713  // New altitude
2714  Numeric znew = z_field(ip,ilat,ilon) + dz;
2715  z_acc = max( z_acc, fabs(znew-z_field(ip+1,ilat,ilon)) );
2716  z_field(ip+1,ilat,ilon) = znew;
2717  }
2718 
2719  // Adjust to z_hse
2720  Vector z_tmp(1);
2721  interp( z_tmp, itw, z_field(joker,ilat,ilon), gp );
2722  z_field(joker,ilat,ilon) -= z_tmp[0] - z_hse[0];
2723  }
2724  }
2725  }
2726 
2727  // Check that there is no gap between the surface and lowest pressure
2728  // level
2729  // (This code is copied from *basics_checkedCalc*. Make this to an internal
2730  // function if used in more places.)
2731  for( Index row=0; row<z_surface.nrows(); row++ )
2732  {
2733  for( Index col=0; col<z_surface.ncols(); col++ )
2734  {
2735  if( z_surface(row,col)<z_field(0,row,col) ||
2736  z_surface(row,col)>=z_field(z_field.npages()-1,row,col) )
2737  {
2738  ostringstream os;
2739  os << "The surface altitude (*z_surface*) cannot be outside "
2740  << "of the altitudes in *z_field*.";
2741  throw runtime_error( os.str() );
2742  }
2743  }
2744  }
2745 
2746 }
2747 
2748 
2749 
2750 
GFIELD3_LON_GRID
const Index GFIELD3_LON_GRID
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
AtmFieldsExpand1D
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 Verbosity &)
WORKSPACE METHOD: AtmFieldsExpand1D.
Definition: m_atmosphere.cc:2101
is_lon_cyclic
bool is_lon_cyclic(ConstVectorView grid, const Numeric &epsilon)
Check if the given longitude grid is cyclic.
Definition: logic.cc:466
GriddedFieldLatLonRegrid
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.
Definition: m_atmosphere.cc:808
GFIELD3_LAT_GRID
const Index GFIELD3_LAT_GRID
Agenda
The Agenda class.
Definition: agenda_class.h:44
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:176
AtmFieldsFromCompact
void AtmFieldsFromCompact(Vector &p_grid, Vector &lat_grid, Vector &lon_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, Tensor4 &massdensity_field, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfString &part_species, const GriddedField4 &atm_fields_compact, const Index &atmosphere_dim, const String &delim, const Verbosity &)
WORKSPACE METHOD: AtmFieldsFromCompact.
Definition: m_atmosphere.cc:1593
GriddedFieldZToPRegridHelper
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.
Definition: m_atmosphere.cc:1122
messages.h
Declarations having to do with the four output streams.
interp_atmfield_by_gp
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)
interp_atmfield_by_gp
Definition: special_interp.cc:215
GriddedField::get_numeric_grid
ConstVectorView get_numeric_grid(Index i) const
Get a numeric grid.
Definition: gridded_fields.cc:93
Vector
The Vector class.
Definition: matpackI.h:556
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:1046
batch_atm_fields_compactFromArrayOfMatrix
void batch_atm_fields_compactFromArrayOfMatrix(ArrayOfGriddedField4 &batch_atm_fields_compact, const Index &atmosphere_dim, const ArrayOfMatrix &am, const ArrayOfString &field_names, const ArrayOfString &extra_field_names, const Vector &extra_field_values, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactFromArrayOfMatrix.
Definition: m_atmosphere.cc:1511
GriddedField::get_string_grid
const ArrayOfString & get_string_grid(Index i) const
Get a string grid.
Definition: gridded_fields.cc:159
p_gridDensify
void p_gridDensify(Vector &p_grid, const Index &nfill, const Verbosity &verbosity)
WORKSPACE METHOD: p_gridDensify.
Definition: m_atmosphere.cc:2398
is_increasing
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:275
Tensor4
The Tensor4 class.
Definition: matpackIV.h:383
GriddedField4::checksize
virtual bool checksize() const
Consistency check.
Definition: gridded_fields.h:382
Range
The range class.
Definition: matpackI.h:148
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
GAS_CONSTANT
const Numeric GAS_CONSTANT
xml_io.h
This file contains basic functions to handle XML data files.
p2gridpos
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
p2gridpos
Definition: special_interp.cc:608
GFIELD3_P_GRID
const Index GFIELD3_P_GRID
interpolation.h
Header file for interpolation.cc.
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:146
GriddedField4::resize
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
Definition: gridded_fields.h:417
ConstTensor4View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:75
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
global_data::species_data
const Array< SpeciesRecord > species_data
Species Data.
Definition: partition_function_data.cc:43
chk_interpolation_grids
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:1094
GridPos
Structure to store a grid position.
Definition: interpolation.h:74
gridpos_poly_longitudinal
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.
Definition: interpolation_poly.cc:302
AtmFieldsCalc
void AtmFieldsCalc(Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_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 Index &atmosphere_dim, const Index &interp_order, const Index &vmr_zeropadding, const Index &vmr_nonegative, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsCalc.
Definition: m_atmosphere.cc:1786
my_basic_string
The implementation for String, the ARTS string class.
Definition: mystring.h:63
PI
const Numeric PI
AtmFieldsCalcExpand1D
void AtmFieldsCalcExpand1D(Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_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 Index &atmosphere_dim, const Index &interp_order, const Index &vmr_zeropadding, const Index &vmr_nonegative, const Verbosity &verbosity)
WORKSPACE METHOD: AtmFieldsCalcExpand1D.
Definition: m_atmosphere.cc:2040
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
GriddedField3::checksize
virtual bool checksize() const
Consistency check.
Definition: gridded_fields.h:322
Tensor3
The Tensor3 class.
Definition: matpackIII.h:348
GriddedFieldLatLonRegridHelper
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. ...
Definition: m_atmosphere.cc:755
arts.h
The global header file for ARTS.
AtmosphereSet1D
void AtmosphereSet1D(Index &atmosphere_dim, Vector &lat_grid, Vector &lon_grid, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet1D.
Definition: m_atmosphere.cc:1730
batch_atm_fields_compactAddSpecies
void batch_atm_fields_compactAddSpecies(ArrayOfGriddedField4 &batch_atm_fields_compact, const String &name, const GriddedField3 &species, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactAddSpecies.
Definition: m_atmosphere.cc:1478
GriddedFieldZToPRegrid
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.
Definition: m_atmosphere.cc:1183
chk_interpolation_grids_loose_no_data_check
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:859
GFIELD4_FIELD_NAMES
const Index GFIELD4_FIELD_NAMES
pos2true_latlon
void pos2true_latlon(Numeric &lat, Numeric &lon, const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true, ConstVectorView pos)
pos2true_latlon
Definition: rte.cc:2438
p2gridpos_poly
void p2gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Index order, const Numeric &extpolfac)
p2gridpos_poly
Definition: special_interp.cc:643
is_same_within_epsilon
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:446
chk_latlon_true
void chk_latlon_true(const Index &atmosphere_dim, ConstVectorView lat_grid, ConstVectorView lat_true, ConstVectorView lon_true)
chk_latlon_true
Definition: check_input.cc:1819
atm_fields_compactFromMatrix
void atm_fields_compactFromMatrix(GriddedField4 &af, const Index &atmosphere_dim, const Matrix &im, const ArrayOfString &field_names, const Verbosity &)
WORKSPACE METHOD: atm_fields_compactFromMatrix.
Definition: m_atmosphere.cc:1281
chk_interpolation_pgrids
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.
Definition: check_input.cc:1266
chk_rte_pos
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
Definition: check_input.cc:1961
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
agenda_class.h
Declarations for agendas.
max
#define max(a, b)
Definition: continua.cc:20461
GriddedField::set_grid
void set_grid(Index i, const Vector &g)
Set a numeric grid.
Definition: gridded_fields.cc:225
my_basic_string::toupper
void toupper()
Convert to upper case.
Definition: mystring.h:87
chk_griddedfield_gridname
void chk_griddedfield_gridname(const GriddedField &gf, const Index gridindex, const String &gridname)
Check name of grid in GriddedField.
Definition: check_input.cc:2096
atm_fields_compactAddSpecies
void atm_fields_compactAddSpecies(GriddedField4 &atm_fields_compact, const String &name, const GriddedField3 &species, const Verbosity &verbosity)
WORKSPACE METHOD: atm_fields_compactAddSpecies.
Definition: m_atmosphere.cc:1380
GriddedField::get_grid_size
Index get_grid_size(Index i) const
Get the size of a grid.
Definition: gridded_fields.h:120
AtmosphereSet2D
void AtmosphereSet2D(Index &atmosphere_dim, Vector &lon_grid, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet2D.
Definition: m_atmosphere.cc:1752
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:212
GriddedField2::resize
void resize(const GriddedField2 &gf)
Make this GriddedField2 the same size as the given one.
Definition: gridded_fields.h:288
FieldFromGriddedField
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.
Definition: m_atmosphere.cc:184
GFIELD4_LON_GRID
const Index GFIELD4_LON_GRID
find_first_species_tg
Index find_first_species_tg(const ArrayOfArrayOfSpeciesTag &tgs, const Index &spec)
Find first occurrence of species in tag groups.
Definition: abs_species_tags.cc:590
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:167
xml_read_from_file
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:831
joker
const Joker joker
GriddedFieldPRegrid
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.
Definition: m_atmosphere.cc:572
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Matrix
The Matrix class.
Definition: matpackI.h:788
auto_md.h
VectorNLogSpace
void VectorNLogSpace(Vector &x, const Index &n, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
Definition: m_basic_types.cc:1040
EPSILON_LON_CYCLIC
const Numeric EPSILON_LON_CYCLIC
Data value accuracy requirement for values at 0 and 360 deg if longitudes are cyclic.
absorption.h
Declarations required for the calculation of absorption coefficients.
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
Definition: interpolation_poly.cc:127
atm_fields_compactExpand
void atm_fields_compactExpand(GriddedField4 &af, Index &nf, const String &name, const Verbosity &)
atm_fields_compactExpand
Definition: m_atmosphere.cc:102
refell2r
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
GFIELD4_LAT_GRID
const Index GFIELD4_LAT_GRID
FieldFromGriddedFieldCheckLatLonHelper
void FieldFromGriddedFieldCheckLatLonHelper(const Vector &lat_grid, const Vector &lon_grid, const Index ilat, const Index ilon, const GriddedField &gfield)
Check for correct grid dimensions.
Definition: m_atmosphere.cc:156
special_interp.h
Header file for special_interp.cc.
Tensor3::resize
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
species_index_from_species_name
Index species_index_from_species_name(String name)
Return species index for given species name.
Definition: absorption.cc:1260
AtmRawRead
void AtmRawRead(GriddedField3 &t_field_raw, GriddedField3 &z_field_raw, ArrayOfGriddedField3 &vmr_field_raw, const ArrayOfArrayOfSpeciesTag &abs_species, const String &basename, const Verbosity &verbosity)
WORKSPACE METHOD: AtmRawRead.
Definition: m_atmosphere.cc:2308
p_gridFromZRaw
void p_gridFromZRaw(Vector &p_grid, const GriddedField3 &z_field_raw, const Index &no_negZ, const Verbosity &)
WORKSPACE METHOD: p_gridFromZRaw.
Definition: m_atmosphere.cc:2434
chk_interpolation_pgrids_loose_no_data_check
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:989
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
Array
This can be used to make arrays out of anything.
Definition: array.h:40
AtmFieldsRefinePgrid
void AtmFieldsRefinePgrid(Vector &p_grid, Tensor3 &t_field, Tensor3 &z_field, Tensor4 &vmr_field, const Vector &lat_grid, const Vector &lon_grid, const Index &atmosphere_dim, const Numeric &p_step, const Verbosity &)
WORKSPACE METHOD: AtmFieldsRefinePgrid.
Definition: m_atmosphere.cc:2154
interpolation_poly.h
Header file for interpolation_poly.cc.
GriddedField::get_grid_name
const String & get_grid_name(Index i) const
Get grid name.
Definition: gridded_fields.h:111
chk_size
void chk_size(const String &x_name, ConstVectorView x, const Index &c)
Runtime check for size of Vector.
Definition: check_input.cc:465
rte_pos2gridpos
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)
rte_pos2gridpos
Definition: special_interp.cc:688
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:101
GFIELD4_P_GRID
const Index GFIELD4_P_GRID
parse_partfield_name
void parse_partfield_name(String &partfield_name, const String &part_string, const String &delim)
Definition: cloudbox.cc:3700
InterpAtmFieldToPosition
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.
Definition: m_atmosphere.cc:2365
Vector::resize
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
lon_gridFromRawField
void lon_gridFromRawField(Vector &lon_grid, const GriddedField3 &field_raw, const Verbosity &)
WORKSPACE METHOD: lon_gridFromRawField.
Definition: m_atmosphere.cc:2487
chk_atm_field
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)
Definition: check_input.cc:1385
DEG2RAD
const Numeric DEG2RAD
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:292
atm_fields_compactAddConstant
void atm_fields_compactAddConstant(GriddedField4 &af, const String &name, const Numeric &value, const Verbosity &verbosity)
WORKSPACE METHOD: atm_fields_compactAddConstant.
Definition: m_atmosphere.cc:1362
abs_species_tags.h
Header file for stuff related to absorption species tags.
g0_agendaExecute
void g0_agendaExecute(Workspace &ws, Numeric &g0, const Numeric lat, const Numeric lon, const Agenda &input_agenda)
Definition: auto_md.cc:13631
chk_if_equal
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:356
GriddedFieldPRegridHelper
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.
Definition: m_atmosphere.cc:507
wind_u_fieldIncludePlanetRotation
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.
Definition: m_atmosphere.cc:2498
GriddedFieldLatLonExpand
void GriddedFieldLatLonExpand(GriddedField2 &gfraw_out, const GriddedField2 &gfraw_in_orig, const Verbosity &)
WORKSPACE METHOD: GriddedFieldLatLonExpand.
Definition: m_atmosphere.cc:276
Workspace
Workspace class.
Definition: workspace_ng.h:47
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
GriddedField::set_grid_name
void set_grid_name(Index i, const String &s)
Set grid name.
Definition: gridded_fields.h:166
lat_gridFromRawField
void lat_gridFromRawField(Vector &lat_grid, const GriddedField3 &field_raw, const Verbosity &)
WORKSPACE METHOD: lat_gridFromRawField.
Definition: m_atmosphere.cc:2476
transform
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1838
_U_
#define _U_
Definition: config.h:167
batch_atm_fields_compactAddConstant
void batch_atm_fields_compactAddConstant(ArrayOfGriddedField4 &batch_atm_fields_compact, const String &name, const Numeric &value, const Verbosity &verbosity)
WORKSPACE METHOD: batch_atm_fields_compactAddConstant.
Definition: m_atmosphere.cc:1462
gridded_fields.h
Implementation of gridded fields.
is_decreasing
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:324
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:214
AtmosphereSet3D
void AtmosphereSet3D(Index &atmosphere_dim, const Verbosity &verbosity)
WORKSPACE METHOD: AtmosphereSet3D.
Definition: m_atmosphere.cc:1770
z_fieldFromHSE
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.
Definition: m_atmosphere.cc:2560
cloudbox.h
Internal cloudbox functions.
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
GriddedField3::resize
void resize(const GriddedField3 &gf)
Make this GriddedField3 the same size as the given one.
Definition: gridded_fields.h:354
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:802
chk_atm_grids
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
Definition: check_input.cc:1305
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:213
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
z2g
void z2g(Numeric &g, const Numeric &r, const Numeric &g0, const Numeric &z)
Definition: m_atmosphere.cc:2549
Tensor4::resize
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1403
global_data.h
rte.h
Declaration of functions in rte.cc.
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
check_input.h