ARTS 2.5.9 (git: 825fa5f2)
m_fluxes.cc
Go to the documentation of this file.
1/* Copyright (C) 2018
2 Manfred Brath <manfred.brath@uni-hamburg.de>
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
19/*===========================================================================
20 === File description
21 ===========================================================================*/
22#include <iostream>
23#include <stdexcept>
24#include "absorption.h"
25#include "agenda_class.h"
26#include "agenda_set.h"
27#include "arts_constants.h"
28#include "arts_conversions.h"
29#include "auto_md.h"
30#include "check_input.h"
31#include "legendre.h"
32#include "math_funcs.h"
33#include "matpackVII.h"
34#include "messages.h"
35#include "sorting.h"
36#include "surface.h"
37#include "workspace_ng.h"
38#include "check_input.h"
39#include "global_data.h"
40
52inline constexpr Numeric PI=Constant::pi;
54
55/*===========================================================================
56 === The functions
57 ===========================================================================*/
58
59/* Workspace method: Doxygen documentation will be auto-generated */
61 Vector& aa_grid,
62 Vector& za_grid_weights,
63 // Keywords:
64 const Index& N_za_grid,
65 const Index& N_aa_grid,
66 const String& za_grid_type,
67 const Verbosity&) {
68 // Azimuth angle grid
69 if (N_aa_grid > 1)
70 nlinspace(aa_grid, -180, 180, N_aa_grid);
71 else if (N_aa_grid < 1) {
72 ostringstream os;
73 os << "N_aa_grid must be > 0 (even for 1D).";
74 throw std::runtime_error(os.str());
75 } else {
76 aa_grid.resize(1);
77 aa_grid[0] = 0.;
78 }
79
80 if (N_za_grid % 2 == 1) {
81 ostringstream os;
82 os << "N_za_grid must be even.";
83 throw runtime_error(os.str());
84 }
85
86 Index nph = N_za_grid / 2;
87
88 //calculate zenith angle grid
89 za_grid.resize(N_za_grid);
90 za_grid = 0.;
91 za_grid_weights.resize(N_za_grid);
92 za_grid_weights = 0;
93
94 if (za_grid_type == "double_gauss") {
95 Vector x;
96 Vector w;
97 Vector xtemp;
98 Vector wtemp;
99 //Numeric theta;
100
101 //calculate legendre weights and evaluation points
102 gsl_integration_glfixed_table_alloc(xtemp, wtemp, nph);
103
104 x.resize(nph);
105 w.resize(nph);
106
107 // reorder and expand weights and abscissa vectors
108 // transform abscissa vector from cos(theta)-space to theta-space
109 // and adjust the domain and weights
110 if (nph % 2 == 1) {
111 x[xtemp.nelem() - 1] = acos((xtemp[0] + 1) / 2) / DEG2RAD;
112 w[wtemp.nelem() - 1] = wtemp[0] / 2;
113
114 for (Index i = 0; i < xtemp.nelem() - 1; i++) {
115 x[i] = acos((xtemp[xtemp.nelem() - 1 - i] + 1) / 2.) / DEG2RAD;
116 x[xtemp.nelem() + i] = acos(1 - (xtemp[i + 1] + 1) / 2.) / DEG2RAD;
117
118 w[i] = wtemp[wtemp.nelem() - 1 - i] / 2;
119 w[wtemp.nelem() + i] = wtemp[i + 1] / 2;
120 }
121 } else {
122 for (Index i = 0; i < xtemp.nelem(); i++) {
123 x[i] = acos((xtemp[xtemp.nelem() - 1 - i] + 1) / 2.) / DEG2RAD;
124 x[xtemp.nelem() + i] = acos(1 - (xtemp[i] + 1) / 2.) / DEG2RAD;
125
126 w[i] = wtemp[wtemp.nelem() - 1 - i] / 2;
127 w[wtemp.nelem() + i] = wtemp[i] / 2;
128 }
129 }
130
131 for (Index i = 0; i < nph; i++) {
132 //set the angles
133 //theta=x[i];//acos((x[i]+1)/2)/DEG2RAD;
134 za_grid[i] = x[i];
135 za_grid[za_grid.nelem() - 1 - i] = 180 - x[i];
136
137 // set the weights to the right component
138 za_grid_weights[i] = w[i];
139 za_grid_weights[za_grid_weights.nelem() - 1 - i] = w[i];
140 }
141
142 } else if (za_grid_type == "linear") {
143 Vector x;
144 Vector w;
146
147 for (Index i = 0; i < N_za_grid; i++) {
148 za_grid[i] = (x[i] + 1) * 90.;
149
150 // set the weights to the right component
151 // by adjusting the domain, we also have to adjust the weights
152 za_grid_weights[i] = w[i] * sin(za_grid[i] * DEG2RAD);
153 }
154 } else if (za_grid_type == "linear_mu") {
155 Vector x;
156 Vector w;
157
158 //calculate weights in cos(theta) space
160
161 //allocate
162 Vector za_grid_temp;
163 za_grid_temp.resize(x.nelem());
164
165 for (Index i = 0; i < N_za_grid; i++) {
166 za_grid_temp[i] = acos(x[i]) / DEG2RAD;
167 }
168
169 //#sort weights and theta in increasing direction of za_grid
170 za_grid = za_grid_temp[Range(x.nelem() - 1, x.nelem(), -1)];
171 za_grid_weights = w[Range(x.nelem() - 1, x.nelem(), -1)];
172
173 } else {
174 ostringstream os;
175 os << "The selected grid type is not implemented";
176 throw std::runtime_error(os.str());
177 }
178
179 //be sure that the first and the last angle are within the closed interval
180 //between 0 and 180 deg, because ARTS is picky if the angles are due to numerics
181 // slightly below and above,respectively.
182 if (za_grid[0] < 0) {
183 za_grid[0] = 0.;
184 }
185
186 if (za_grid[za_grid.nelem() - 1] > 180) {
187 za_grid[za_grid.nelem() - 1] = 180.;
188 }
189}
190
191/* Workspace method: Doxygen documentation will be auto-generated */
193 const Vector& p_grid,
194 const Tensor4& irradiance_field,
195 const Tensor3& specific_heat_capacity,
196 const Numeric& g0,
197 const Verbosity&) {
198 //allocate
199 heating_rates.resize(irradiance_field.nbooks(),
200 irradiance_field.npages(),
201 irradiance_field.nrows());
202 heating_rates = 0;
203
204 // allocate some auxiliary variables
205 Numeric net_flux_b; //net_flux bottom
206 Numeric net_flux_c; //net_flux center
207 Numeric net_flux_t; //net_flux top
208 Index idx;
209
210 // calculate heating rates, we skip the upper and lower boundary here, because to achieve the same
211 // second order accuracy for the derivation of the net flux at the edged, we use
212 // a differentiation based on polynomial interpolation
213 for (Index b = 1; b < irradiance_field.nbooks() - 1; b++) {
214 for (Index p = 0; p < irradiance_field.npages(); p++) {
215 for (Index r = 0; r < irradiance_field.nrows(); r++) {
216 net_flux_b = (irradiance_field(b - 1, p, r, 0) +
217 irradiance_field(b - 1, p, r, 1));
218 net_flux_t = (irradiance_field(b + 1, p, r, 0) +
219 irradiance_field(b + 1, p, r, 1));
220
221 heating_rates(b, p, r) = (net_flux_t - net_flux_b) /
222 (p_grid[b + 1] - p_grid[b - 1]) * g0 /
223 specific_heat_capacity(b, p, r);
224 }
225 }
226 }
227
228 idx = irradiance_field.nbooks();
229
230 // now calculate the heating rates for the upper and lower boundary
231 for (Index p = 0; p < irradiance_field.npages(); p++) {
232 for (Index r = 0; r < irradiance_field.nrows(); r++) {
233 // lower boundary
234 net_flux_b =
235 (irradiance_field(0, p, r, 0) + irradiance_field(0, p, r, 1));
236 net_flux_c =
237 (irradiance_field(1, p, r, 0) + irradiance_field(1, p, r, 1));
238 net_flux_t =
239 (irradiance_field(2, p, r, 0) + irradiance_field(0, p, r, 1));
240
241 heating_rates(0, p, r) = (-3 * net_flux_b + 4 * net_flux_c - net_flux_t) /
242 (p_grid[2] - p_grid[0]) * g0 /
243 specific_heat_capacity(0, p, r);
244
245 // upper boundary
246 net_flux_t = (irradiance_field(idx - 1, p, r, 0) +
247 irradiance_field(idx - 1, p, r, 1));
248 net_flux_c = (irradiance_field(idx - 2, p, r, 0) +
249 irradiance_field(idx - 2, p, r, 1));
250 net_flux_b = (irradiance_field(idx - 3, p, r, 0) +
251 irradiance_field(idx - 3, p, r, 1));
252
253 heating_rates(idx - 1, p, r) =
254 -(-3 * net_flux_t + 4 * net_flux_c - net_flux_b) /
255 (p_grid[idx-1] - p_grid[idx-3]) * g0 / specific_heat_capacity(0, p, r);
256 }
257 }
258}
259
260/* Workspace method: Doxygen documentation will be auto-generated */
262 const Tensor5& radiance_field,
263 const Vector& za_grid,
264 const Vector& aa_grid,
265 const Vector& za_grid_weights,
266 const Verbosity&) {
267 // Number of zenith angles.
268 const Index N_scat_za = za_grid.nelem();
269 const Index N_scat_aa = aa_grid.nelem();
270
271 Tensor4 radiance_field_aa_integrated;
272
273 //azimuth integration
274 if (N_scat_aa == 1) //1D case no azimuth dependency
275 {
276 radiance_field_aa_integrated =
277 radiance_field(joker, joker, joker, joker, 0);
278 radiance_field_aa_integrated *= 2 * PI;
279
280 } else //general case with azimuth dependency
281 {
282 radiance_field_aa_integrated.resize(radiance_field.nshelves(),
283 radiance_field.nbooks(),
284 radiance_field.npages(),
285 radiance_field.nrows());
286 radiance_field_aa_integrated = 0.;
287
288 for (Index b = 0; b < radiance_field_aa_integrated.nbooks(); b++) {
289 for (Index p = 0; p < radiance_field_aa_integrated.npages(); p++) {
290 for (Index r = 0; r < radiance_field_aa_integrated.nrows(); r++) {
291 for (Index c = 0; c < radiance_field_aa_integrated.ncols(); c++) {
292 for (Index i = 0; i < N_scat_aa - 1; i++) {
293 radiance_field_aa_integrated(b, p, r, c) +=
294 (radiance_field(b, p, r, c, i) +
295 radiance_field(b, p, r, c, i + 1)) /
296 2. * abs(aa_grid[i + 1] - aa_grid[i]) * DEG2RAD;
297 }
298 }
299 }
300 }
301 }
302 }
303
304 //allocate
305 irradiance_field.resize(radiance_field.nshelves(),
306 radiance_field.nbooks(),
307 radiance_field.npages(),
308 2);
309 irradiance_field = 0;
310
311 // zenith angle integration
312
313 for (Index b = 0; b < irradiance_field.nbooks(); b++) {
314 for (Index p = 0; p < irradiance_field.npages(); p++) {
315 for (Index r = 0; r < irradiance_field.nrows(); r++) {
316 for (Index i = 0; i < N_scat_za; i++) {
317 if (za_grid[i] <= 90.) {
318 irradiance_field(b, p, r, 0) +=
319 radiance_field_aa_integrated(b, p, r, i) *
320 cos(za_grid[i] * DEG2RAD) * (-1.) * za_grid_weights[i];
321 } else {
322 irradiance_field(b, p, r, 1) +=
323 radiance_field_aa_integrated(b, p, r, i) *
324 cos(za_grid[i] * DEG2RAD) * (-1.) * za_grid_weights[i];
325 }
326 }
327 }
328 }
329 }
330}
331
332/* Workspace method: Doxygen documentation will be auto-generated */
334 const Vector& f_grid,
335 const Tensor5& spectral_radiation_field,
336 const Verbosity&) {
337 if (f_grid.nelem() != spectral_radiation_field.nshelves()) {
338 throw runtime_error(
339 "The length of f_grid does not match with\n"
340 " the first dimension of the spectral_radiation_field");
341 }
342
343 //allocate
344 radiation_field.resize(spectral_radiation_field.nbooks(),
345 spectral_radiation_field.npages(),
346 spectral_radiation_field.nrows(),
347 spectral_radiation_field.ncols());
348 radiation_field = 0;
349
350 // frequency integration
351 for (Index i = 0; i < spectral_radiation_field.nshelves() - 1; i++) {
352 const Numeric df = f_grid[i + 1] - f_grid[i];
353
354 for (Index b = 0; b < radiation_field.nbooks(); b++) {
355 for (Index p = 0; p < radiation_field.npages(); p++) {
356 for (Index r = 0; r < radiation_field.nrows(); r++) {
357 for (Index c = 0; c < radiation_field.ncols(); c++) {
358 radiation_field(b, p, r, c) +=
359 (spectral_radiation_field(i + 1, b, p, r, c) +
360 spectral_radiation_field(i, b, p, r, c)) /
361 2 * df;
362 }
363 }
364 }
365 }
366 }
367}
368
369/* Workspace method: Doxygen documentation will be auto-generated */
371 const Vector& f_grid,
372 const Tensor7& spectral_radiation_field,
373 const Verbosity&) {
374 if (f_grid.nelem() != spectral_radiation_field.nlibraries()) {
375 throw runtime_error(
376 "The length of f_grid does not match with\n"
377 " the first dimension of the spectral_radiation_field");
378 }
379
380 //allocate
381 radiation_field.resize(spectral_radiation_field.nvitrines(),
382 spectral_radiation_field.nshelves(),
383 spectral_radiation_field.nbooks(),
384 spectral_radiation_field.npages(),
385 spectral_radiation_field.nrows());
386 radiation_field = 0;
387
388 // frequency integration
389 for (Index i = 0; i < spectral_radiation_field.nlibraries() - 1; i++) {
390 const Numeric df = f_grid[i + 1] - f_grid[i];
391
392 for (Index s = 0; s < radiation_field.nshelves(); s++) {
393 for (Index b = 0; b < radiation_field.nbooks(); b++) {
394 for (Index p = 0; p < radiation_field.npages(); p++) {
395 for (Index r = 0; r < radiation_field.nrows(); r++) {
396 for (Index c = 0; c < radiation_field.ncols(); c++) {
397 radiation_field(s, b, p, r, c) +=
398 (spectral_radiation_field(i + 1, s, b, p, r, c, 0) +
399 spectral_radiation_field(i, s, b, p, r, c, 0)) /
400 2 * df;
401 }
402 }
403 }
404 }
405 }
406 }
407}
408
409/* Workspace method: Doxygen documentation will be auto-generated */
411 Tensor5& spectral_irradiance_field,
412 const Tensor7& spectral_radiance_field,
413 const Vector& za_grid,
414 const Vector& aa_grid,
415 const Vector& za_grid_weights,
416 const Verbosity&) {
417 // Number of zenith angles.
418 const Index N_scat_za = spectral_radiance_field.npages();
419 const Index N_scat_aa = spectral_radiance_field.nrows();
420
421 Tensor5 iy_field_aa_integrated;
422
423 //azimuth integration
424 if (N_scat_aa == 1) //1D case no azimuth dependency
425 {
426 iy_field_aa_integrated =
427 spectral_radiance_field(joker, joker, joker, joker, joker, 0, 0);
428 iy_field_aa_integrated *= 2 * PI;
429
430 } else //general case with azimuth dependency
431 {
432 iy_field_aa_integrated.resize(spectral_radiance_field.nlibraries(),
433 spectral_radiance_field.nvitrines(),
434 spectral_radiance_field.nshelves(),
435 spectral_radiance_field.nbooks(),
436 spectral_radiance_field.npages());
437 iy_field_aa_integrated = 0.;
438
439 for (Index s = 0; s < iy_field_aa_integrated.nshelves(); s++) {
440 for (Index b = 0; b < iy_field_aa_integrated.nbooks(); b++) {
441 for (Index p = 0; p < iy_field_aa_integrated.npages(); p++) {
442 for (Index r = 0; r < iy_field_aa_integrated.nrows(); r++) {
443 for (Index c = 0; c < iy_field_aa_integrated.ncols(); c++) {
444 for (Index i = 0; i < N_scat_aa - 1; i++) {
445 iy_field_aa_integrated(s, b, p, r, c) +=
446 (spectral_radiance_field(s, b, p, r, c, i, 0) +
447 spectral_radiance_field(s, b, p, r, c, i + 1, 0)) /
448 2. * abs(aa_grid[i + 1] - aa_grid[i]) * DEG2RAD;
449 }
450 }
451 }
452 }
453 }
454 }
455 }
456
457 //allocate
458 spectral_irradiance_field.resize(spectral_radiance_field.nlibraries(),
459 spectral_radiance_field.nvitrines(),
460 spectral_radiance_field.nshelves(),
461 spectral_radiance_field.nbooks(),
462 2);
463 spectral_irradiance_field = 0;
464
465 // zenith angle integration
466 for (Index s = 0; s < spectral_irradiance_field.nshelves(); s++) {
467 for (Index b = 0; b < spectral_irradiance_field.nbooks(); b++) {
468 for (Index p = 0; p < spectral_irradiance_field.npages(); p++) {
469 for (Index r = 0; r < spectral_irradiance_field.nrows(); r++) {
470 for (Index i = 0; i < N_scat_za; i++) {
471 if (za_grid[i] <= 90.) {
472 spectral_irradiance_field(s, b, p, r, 0) +=
473 iy_field_aa_integrated(s, b, p, r, i) *
474 cos(za_grid[i] * DEG2RAD) * (-1.) * za_grid_weights[i];
475 } else {
476 spectral_irradiance_field(s, b, p, r, 1) +=
477 iy_field_aa_integrated(s, b, p, r, i) *
478 cos(za_grid[i] * DEG2RAD) * (-1.) * za_grid_weights[i];
479 }
480 }
481 }
482 }
483 }
484 }
485}
486
487/* Workspace method: Doxygen documentation will be auto-generated */
489 Workspace& ws,
490 Tensor7& spectral_radiance_field,
491 Tensor3& trans_field,
492 const Agenda& propmat_clearsky_agenda,
493 const Agenda& water_p_eq_agenda,
494 const Agenda& iy_space_agenda,
495 const Agenda& iy_surface_agenda,
496 const Agenda& iy_cloudbox_agenda,
497 const Index& stokes_dim,
498 const Vector& f_grid,
499 const Index& atmosphere_dim,
500 const Vector& p_grid,
501 const Tensor3& z_field,
502 const Tensor3& t_field,
503 const EnergyLevelMap& nlte_field,
504 const Tensor4& vmr_field,
505 const ArrayOfArrayOfSpeciesTag& abs_species,
506 const Tensor3& wind_u_field,
507 const Tensor3& wind_v_field,
508 const Tensor3& wind_w_field,
509 const Tensor3& mag_u_field,
510 const Tensor3& mag_v_field,
511 const Tensor3& mag_w_field,
512 const Matrix& z_surface,
513 const Numeric& ppath_lmax,
514 const Numeric& rte_alonglos_v,
515 const String& rt_integration_option,
516 const Tensor3& surface_props_data,
517 const Vector& za_grid,
518 const Index& use_parallel_za [[maybe_unused]],
519 const Verbosity& verbosity) {
520 // Check input
521 if (atmosphere_dim != 1)
522 throw runtime_error("This method only works for atmosphere_dim = 1.");
523
524 // Sizes
525 const Index nl = p_grid.nelem();
526 const Index nf = f_grid.nelem();
527 const Index nza = za_grid.nelem();
528
529 // Init spectral_radiance_field and trans_field
530 spectral_radiance_field.resize(nf, nl, 1, 1, nza, 1, stokes_dim);
531 trans_field.resize(nf, nl, nza);
532
533 // De-activate cloudbox
534 const Index cloudbox_on = 0, ppath_inside_cloudbox_do = 0;
535 const ArrayOfIndex cloudbox_limits(0);
536
537 // Various variables
538 const String iy_unit = "1";
539 const ArrayOfString iy_aux_vars(0);
540 const Vector rte_pos2(0);
541 const Index iy_agenda_call1 = 1;
542 const Tensor3 iy_transmittance(0, 0, 0);
543 const Index jacobian_do = 0;
544 const ArrayOfRetrievalQuantity jacobian_quantities(0);
545 // Create one altitude just above TOA
546 const Numeric z_space = z_field(nl - 1, 0, 0) + 10;
547
548 ArrayOfString fail_msg;
549 bool failed = false;
550
551 // Define iy_main_agenda to be consistent with the assumptions of
552 // this method. This definition of iy_main_agenda will be used to when
553 // calculating the the radiation reflected by the surface
554 const Agenda iy_main_agenda=AgendaManip::get_iy_main_agenda(ws, "EmissionPlaneParallel");
555
556 // Index in p_grid where field at surface shall be placed
557 const Index i0 = index_of_zsurface(z_surface(0, 0), z_field(joker, 0, 0));
558
559 // Loop zenith angles
560 //
561 if (nza) {
563
564#pragma omp parallel for if (!arts_omp_in_parallel() && nza > 1 && \
565 use_parallel_za) firstprivate(wss)
566 for (Index i = 0; i < nza; i++) {
567 if (failed) continue;
568 try {
569 // Define output variables
570 Ppath ppath;
571 Vector ppvar_p, ppvar_t;
572 Matrix iy, ppvar_vmr, ppvar_wind, ppvar_mag, ppvar_f;
573 EnergyLevelMap ppvar_nlte;
574 Tensor3 ppvar_iy;
575 Tensor4 ppvar_trans_cumulat, ppvar_trans_partial;
576 ArrayOfMatrix iy_aux;
577 ArrayOfTensor3 diy_dx;
578
579 Index iy_id = i;
580 Vector rte_los(1, za_grid[i]);
581 Vector rte_pos(1, za_grid[i] < 90 ? z_surface(0, 0) : z_space);
582
583 ppathPlaneParallel(ppath,
584 atmosphere_dim,
585 z_field,
586 z_surface,
587 cloudbox_on,
588 cloudbox_limits,
589 ppath_inside_cloudbox_do,
590 rte_pos,
591 rte_los,
592 ppath_lmax,
593 verbosity);
594 ARTS_ASSERT(ppath.gp_p[ppath.np - 1].idx == i0 ||
595 ppath.gp_p[ppath.np - 1].idx == nl - 2);
596
598 iy,
599 iy_aux,
600 diy_dx,
601 ppvar_p,
602 ppvar_t,
603 ppvar_nlte,
604 ppvar_vmr,
605 ppvar_wind,
606 ppvar_mag,
607 ppvar_f,
608 ppvar_iy,
609 ppvar_trans_cumulat,
610 ppvar_trans_partial,
611 iy_id,
612 stokes_dim,
613 f_grid,
614 atmosphere_dim,
615 p_grid,
616 t_field,
617 nlte_field,
618 vmr_field,
619 abs_species,
620 wind_u_field,
621 wind_v_field,
622 wind_w_field,
623 mag_u_field,
624 mag_v_field,
625 mag_w_field,
626 cloudbox_on,
627 iy_unit,
628 iy_aux_vars,
629 jacobian_do,
630 jacobian_quantities,
631 ppath,
632 rte_pos2,
633 propmat_clearsky_agenda,
634 water_p_eq_agenda,
635 rt_integration_option,
636 iy_main_agenda,
637 iy_space_agenda,
638 iy_surface_agenda,
639 iy_cloudbox_agenda,
640 iy_agenda_call1,
641 iy_transmittance,
642 rte_alonglos_v,
643 surface_props_data,
644 verbosity);
645 ARTS_ASSERT(iy.ncols() == stokes_dim);
646
647 // First and last points are most easily handled separately
648 if (za_grid[i] < 90) {
649 spectral_radiance_field(joker, i0, 0, 0, i, 0, joker) =
650 ppvar_iy(joker, joker, 0);
651 spectral_radiance_field(joker, nl - 1, 0, 0, i, 0, joker) =
652 ppvar_iy(joker, joker, ppath.np - 1);
653 trans_field(joker, 0, i) = ppvar_trans_partial(0, joker, 0, 0);
654 trans_field(joker, nl - 1, i) =
655 ppvar_trans_partial(ppath.np - 1, joker, 0, 0);
656 } else {
657 spectral_radiance_field(joker, nl - 1, 0, 0, i, 0, joker) =
658 ppvar_iy(joker, joker, 0);
659 spectral_radiance_field(joker, i0, 0, 0, i, 0, joker) =
660 ppvar_iy(joker, joker, ppath.np - 1);
661 trans_field(joker, nl - 1, i) = ppvar_trans_partial(0, joker, 0, 0);
662 trans_field(joker, 0, i) =
663 ppvar_trans_partial(ppath.np - 1, joker, 0, 0);
664 }
665
666 // Remaining points
667 for (Index p = 1; p < ppath.np - 1; p++) {
668 // We just store values at pressure levels
669 if (ppath.gp_p[p].fd[0] < 1e-2) {
670 spectral_radiance_field(
671 joker, ppath.gp_p[p].idx, 0, 0, i, 0, joker) =
672 ppvar_iy(joker, joker, p);
673 trans_field(joker, ppath.gp_p[p].idx, i) =
674 ppvar_trans_partial(p, joker, 0, 0);
675 }
676 }
677
678 // We don't want undefined values to possibly affect an interpolation,
679 // and to be safe we set the fields for underground levels to equal the
680 // ones at the surface
681 for (Index p = 0; p < i0; p++) {
682 spectral_radiance_field(joker, p, 0, 0, i, 0, joker) =
683 spectral_radiance_field(joker, i0, 0, 0, i, 0, joker);
684 trans_field(joker, p, i) = trans_field(joker, i0, i);
685 }
686
687 } catch (const std::exception& e) {
688#pragma omp critical(planep_setabort)
689 failed = true;
690
691 ostringstream os;
692 os << "Run-time error at nza #" << i << ": \n" << e.what();
693#pragma omp critical(planep_push_fail_msg)
694 fail_msg.push_back(os.str());
695 }
696 }
697 }
698
699 if (fail_msg.nelem()) {
700 ostringstream os;
701 for (auto& msg : fail_msg) os << msg << '\n';
702 throw runtime_error(os.str());
703 }
704}
705
706/* Workspace method: Doxygen documentation will be auto-generated */
708 Tensor7& spectral_radiance_field,
709 const Index& atmosphere_dim,
710 const Vector& p_grid,
711 const Index& cloudbox_on,
712 const ArrayOfIndex& cloudbox_limits,
713 const Tensor7& cloudbox_field,
714 const Verbosity&) {
715 if (atmosphere_dim > 1)
716 throw runtime_error("This method can only be used for 1D calculations.\n");
717 if (!cloudbox_on)
718 throw runtime_error("Cloudbox is off. This is not handled by this method.");
719 if (cloudbox_limits[0] || cloudbox_limits[1] != p_grid.nelem() - 1)
720 throw runtime_error(
721 "The cloudbox must cover all pressure levels "
722 "to use this method.");
723
724 // If all OK, it is just to copy
725 spectral_radiance_field = cloudbox_field;
726}
727
728/* Workspace method: Doxygen documentation will be auto-generated */
730 Workspace& ws,
731 Tensor7& spectral_radiance_field,
732 const Agenda& propmat_clearsky_agenda,
733 const Agenda& water_p_eq_agenda,
734 const Agenda& iy_space_agenda,
735 const Agenda& iy_surface_agenda,
736 const Agenda& iy_cloudbox_agenda,
737 const Index& stokes_dim,
738 const Vector& f_grid,
739 const Index& atmosphere_dim,
740 const Vector& p_grid,
741 const Tensor3& z_field,
742 const Tensor3& t_field,
743 const EnergyLevelMap& nlte_field,
744 const Tensor4& vmr_field,
745 const ArrayOfArrayOfSpeciesTag& abs_species,
746 const Tensor3& wind_u_field,
747 const Tensor3& wind_v_field,
748 const Tensor3& wind_w_field,
749 const Tensor3& mag_u_field,
750 const Tensor3& mag_v_field,
751 const Tensor3& mag_w_field,
752 const Matrix& z_surface,
753 const Index& cloudbox_on,
754 const ArrayOfIndex& cloudbox_limits,
755 const Tensor7& cloudbox_field,
756 const Numeric& ppath_lmax,
757 const Numeric& rte_alonglos_v,
758 const String& rt_integration_option,
759 const Tensor3& surface_props_data,
760 const Vector& za_grid,
761 const Index& use_parallel_za [[maybe_unused]],
762 const Verbosity& verbosity) {
763 // Check input
764 if (atmosphere_dim != 1)
765 throw runtime_error("This method only works for atmosphere_dim = 1.");
766 if (!cloudbox_on)
767 throw runtime_error("No ned to use this method with cloudbox=0.");
768 if (cloudbox_limits[0])
769 throw runtime_error(
770 "The first element of *cloudbox_limits* must be zero "
771 "to use this method.");
772
773 // Sizes
774 const Index nl = p_grid.nelem();
775 const Index nf = f_grid.nelem();
776 const Index nza = za_grid.nelem();
777
778 // Init spectral_radiance_field
779 spectral_radiance_field.resize(nf, nl, 1, 1, nza, 1, stokes_dim);
780
781 // and copy the part taken from cloudbox_field
782 spectral_radiance_field(joker,
783 Range(0, cloudbox_limits[1] + 1),
784 joker,
785 joker,
786 joker,
787 joker,
788 joker) = cloudbox_field;
789
790 // Various variables
791 const Index ppath_inside_cloudbox_do = 0;
792 const String iy_unit = "1";
793 const ArrayOfString iy_aux_vars(0);
794 const Vector rte_pos2(0);
795 const Index iy_agenda_call1 = 1;
796 const Tensor3 iy_transmittance(0, 0, 0);
797 const Index jacobian_do = 0;
798 const ArrayOfRetrievalQuantity jacobian_quantities(0);
799 // Create one altitude just above TOA
800 const Numeric z_space = z_field(nl - 1, 0, 0) + 10;
801
802 ArrayOfString fail_msg;
803 bool failed = false;
804
805 // Define iy_main_agenda to be consistent with the assumptions of
806 // this method (but the agenda will not be used).
807 const Agenda iy_main_agenda=AgendaManip::get_iy_main_agenda(ws, "EmissionPlaneParallel");;
808
809 // Variables related to the top of the cloudbox
810 const Index i0 = cloudbox_limits[1];
811 const Numeric z_top = z_field(i0 + 1, 0, 0); // Note i0+1
812
813 // Loop zenith angles
814 //
815 if (nza) {
817
818#pragma omp parallel for if (!arts_omp_in_parallel() && nza > 1 && \
819 use_parallel_za) firstprivate(wss)
820 for (Index i = 0; i < nza; i++) {
821 if (failed) continue;
822 try {
823 // Define output variables
824 Ppath ppath;
825 Vector ppvar_p, ppvar_t;
826 Matrix iy, ppvar_vmr, ppvar_wind, ppvar_mag, ppvar_f;
827 EnergyLevelMap ppvar_nlte;
828 Tensor3 ppvar_iy;
829 Tensor4 ppvar_trans_cumulat, ppvar_trans_partial;
830 ArrayOfMatrix iy_aux;
831 ArrayOfTensor3 diy_dx;
832
833 Index iy_id = i;
834 Vector rte_los(1, za_grid[i]);
835 Vector rte_pos(1, za_grid[i] < 90 ? z_top : z_space);
836
837 ppathPlaneParallel(ppath,
838 atmosphere_dim,
839 z_field,
840 z_surface,
841 cloudbox_on,
842 cloudbox_limits,
843 ppath_inside_cloudbox_do,
844 rte_pos,
845 rte_los,
846 ppath_lmax,
847 verbosity);
848 ARTS_ASSERT(ppath.gp_p[ppath.np - 1].idx == i0 ||
849 ppath.gp_p[ppath.np - 1].idx == nl - 2);
850
852 iy,
853 iy_aux,
854 diy_dx,
855 ppvar_p,
856 ppvar_t,
857 ppvar_nlte,
858 ppvar_vmr,
859 ppvar_wind,
860 ppvar_mag,
861 ppvar_f,
862 ppvar_iy,
863 ppvar_trans_cumulat,
864 ppvar_trans_partial,
865 iy_id,
866 stokes_dim,
867 f_grid,
868 atmosphere_dim,
869 p_grid,
870 t_field,
871 nlte_field,
872 vmr_field,
873 abs_species,
874 wind_u_field,
875 wind_v_field,
876 wind_w_field,
877 mag_u_field,
878 mag_v_field,
879 mag_w_field,
880 cloudbox_on,
881 iy_unit,
882 iy_aux_vars,
883 jacobian_do,
884 jacobian_quantities,
885 ppath,
886 rte_pos2,
887 propmat_clearsky_agenda,
888 water_p_eq_agenda,
889 rt_integration_option,
890 iy_main_agenda,
891 iy_space_agenda,
892 iy_surface_agenda,
893 iy_cloudbox_agenda,
894 iy_agenda_call1,
895 iy_transmittance,
896 rte_alonglos_v,
897 surface_props_data,
898 verbosity);
899 ARTS_ASSERT(iy.ncols() == stokes_dim);
900
901 // First and last points are most easily handled separately
902 // But field at top cloudbox already known from copying above
903 if (za_grid[i] < 90) {
904 spectral_radiance_field(joker, i0 + 1, 0, 0, i, 0, joker) =
905 ppvar_iy(joker, joker, 0);
906 spectral_radiance_field(joker, nl - 1, 0, 0, i, 0, joker) =
907 ppvar_iy(joker, joker, ppath.np - 1);
908 } else {
909 spectral_radiance_field(joker, nl - 1, 0, 0, i, 0, joker) =
910 ppvar_iy(joker, joker, 0);
911 }
912
913 // Remaining points
914 for (Index p = 1; p < ppath.np - 1; p++) {
915 // We just store values at pressure levels
916 if (ppath.gp_p[p].fd[0] < 1e-2) {
917 spectral_radiance_field(
918 joker, ppath.gp_p[p].idx, 0, 0, i, 0, joker) =
919 ppvar_iy(joker, joker, p);
920 }
921 }
922
923 } catch (const std::exception& e) {
924#pragma omp critical(planep_setabort)
925 failed = true;
926
927 ostringstream os;
928 os << "Run-time error at nza #" << i << ": \n" << e.what();
929#pragma omp critical(planep_push_fail_msg)
930 fail_msg.push_back(os.str());
931 }
932 }
933 }
934
935 if (fail_msg.nelem()) {
936 ostringstream os;
937 for (auto& msg : fail_msg) os << msg << '\n';
938 throw runtime_error(os.str());
939 }
940}
Declarations required for the calculation of absorption coefficients.
Declarations for agendas.
Constants of physical expressions as constexpr.
Common ARTS conversions.
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index ncols() const noexcept
Definition: matpackI.h:1080
Index ncols() const noexcept
Definition: matpackIV.h:146
Index nrows() const noexcept
Definition: matpackIV.h:145
Index nbooks() const noexcept
Definition: matpackIV.h:143
Index npages() const noexcept
Definition: matpackIV.h:144
Index nrows() const noexcept
Definition: matpackV.h:157
Index ncols() const noexcept
Definition: matpackV.h:158
Index npages() const noexcept
Definition: matpackV.h:156
Index nbooks() const noexcept
Definition: matpackV.h:155
Index nshelves() const noexcept
Definition: matpackV.h:154
Index npages() const noexcept
Definition: matpackVII.h:162
Index nrows() const noexcept
Definition: matpackVII.h:163
Index nlibraries() const noexcept
Definition: matpackVII.h:158
Index nvitrines() const noexcept
Definition: matpackVII.h:159
Index nshelves() const noexcept
Definition: matpackVII.h:160
Index nbooks() const noexcept
Definition: matpackVII.h:161
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The Matrix class.
Definition: matpackI.h:1285
The range class.
Definition: matpackI.h:160
The Tensor3 class.
Definition: matpackIII.h:352
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
The Tensor4 class.
Definition: matpackIV.h:435
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1054
The Tensor5 class.
Definition: matpackV.h:524
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1741
The Tensor7 class.
Definition: matpackVII.h:2407
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5488
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Workspace class.
Definition: workspace_ng.h:53
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
bool gsl_integration_glfixed_table_alloc(Vector &x, Vector &w, Index n)
gsl_integration_glfixed_table_alloc
Definition: legendre.cc:1940
Contains the code to calculate Legendre polynomials.
void spectral_irradiance_fieldFromSpectralRadianceField(Tensor5 &spectral_irradiance_field, const Tensor7 &spectral_radiance_field, const Vector &za_grid, const Vector &aa_grid, const Vector &za_grid_weights, const Verbosity &)
WORKSPACE METHOD: spectral_irradiance_fieldFromSpectralRadianceField.
Definition: m_fluxes.cc:410
constexpr Numeric DEG2RAD
Definition: m_fluxes.cc:53
void heating_ratesFromIrradiance(Tensor3 &heating_rates, const Vector &p_grid, const Tensor4 &irradiance_field, const Tensor3 &specific_heat_capacity, const Numeric &g0, const Verbosity &)
WORKSPACE METHOD: heating_ratesFromIrradiance.
Definition: m_fluxes.cc:192
void irradiance_fieldFromRadiance(Tensor4 &irradiance_field, const Tensor5 &radiance_field, const Vector &za_grid, const Vector &aa_grid, const Vector &za_grid_weights, const Verbosity &)
WORKSPACE METHOD: irradiance_fieldFromRadiance.
Definition: m_fluxes.cc:261
void spectral_radiance_fieldClearskyPlaneParallel(Workspace &ws, Tensor7 &spectral_radiance_field, Tensor3 &trans_field, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Matrix &z_surface, const Numeric &ppath_lmax, const Numeric &rte_alonglos_v, const String &rt_integration_option, const Tensor3 &surface_props_data, const Vector &za_grid, const Index &use_parallel_za, const Verbosity &verbosity)
WORKSPACE METHOD: spectral_radiance_fieldClearskyPlaneParallel.
Definition: m_fluxes.cc:488
void spectral_radiance_fieldCopyCloudboxField(Tensor7 &spectral_radiance_field, const Index &atmosphere_dim, const Vector &p_grid, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor7 &cloudbox_field, const Verbosity &)
WORKSPACE METHOD: spectral_radiance_fieldCopyCloudboxField.
Definition: m_fluxes.cc:707
void spectral_radiance_fieldExpandCloudboxField(Workspace &ws, Tensor7 &spectral_radiance_field, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &z_field, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor7 &cloudbox_field, const Numeric &ppath_lmax, const Numeric &rte_alonglos_v, const String &rt_integration_option, const Tensor3 &surface_props_data, const Vector &za_grid, const Index &use_parallel_za, const Verbosity &verbosity)
WORKSPACE METHOD: spectral_radiance_fieldExpandCloudboxField.
Definition: m_fluxes.cc:729
constexpr Numeric PI
Definition: m_fluxes.cc:52
void RadiationFieldSpectralIntegrate(Tensor4 &radiation_field, const Vector &f_grid, const Tensor5 &spectral_radiation_field, const Verbosity &)
WORKSPACE METHOD: RadiationFieldSpectralIntegrate.
Definition: m_fluxes.cc:333
void AngularGridsSetFluxCalc(Vector &za_grid, Vector &aa_grid, Vector &za_grid_weights, const Index &N_za_grid, const Index &N_aa_grid, const String &za_grid_type, const Verbosity &)
WORKSPACE METHOD: AngularGridsSetFluxCalc.
Definition: m_fluxes.cc:60
void ppathPlaneParallel(Ppath &ppath, const Index &atmosphere_dim, const Tensor3 &z_field, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Verbosity &)
WORKSPACE METHOD: ppathPlaneParallel.
Definition: m_ppath.cc:1202
void iyEmissionStandard(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandard.
Definition: m_rte.cc:1388
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:812
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
Agenda get_iy_main_agenda(Workspace &ws, const String &option)
Definition: agenda_set.cc:225
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
Contains sorting routines.
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath_struct.h:51
Index index_of_zsurface(const Numeric &z_surface, ConstVectorView z_profile)
Lccates the surface with respect to pressure levels.
Definition: surface.cc:55
#define w
#define c
#define b
This file contains the Workspace class.