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