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