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