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