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