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