ARTS  2.4.0(git:4fb77825)
m_microphysics.cc
Go to the documentation of this file.
1 /* Copyright (C) 2011-2017
2  Jana Mendrok <jana.mendrok@gmail.com>
3  Daniel Kreyling <daniel.kreyling@nict.go.jp>
4  Manfred Brath <manfred.brath@uni-hamburg.de>
5  Patrick Eriksson <patrick.eriksson@chalmers.se>
6 
7  This program is free software; you can redistribute it and/or modify it
8  under the terms of the GNU General Public License as published by the
9  Free Software Foundation; either version 2, or (at your option) any
10  later version.
11 
12  This program is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with this program; if not, write to the Free Software
19  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
20  USA. */
21 
22 /*===========================================================================
23  === File description
24  ===========================================================================*/
25 
37 /*===========================================================================
38  === External declarations
39  ===========================================================================*/
40 #include <cmath>
41 #include <cstdlib>
42 #include <stdexcept>
43 
44 #include "array.h"
45 #include "arts.h"
46 #include "auto_md.h"
47 #include "check_input.h"
48 #include "cloudbox.h"
49 #include "file.h"
50 #include "interpolation.h"
51 #include "lin_alg.h"
52 #include "logic.h"
53 #include "math_funcs.h"
54 #include "messages.h"
55 #include "microphysics.h"
56 #include "optproperties.h"
57 #include "parameters.h"
58 #include "physics_funcs.h"
59 #include "psd.h"
60 #include "rte.h"
61 #include "sorting.h"
62 #include "special_interp.h"
63 #include "xml_io.h"
64 
65 extern const String SCATSPECIES_MAINTAG;
66 
67 /*===========================================================================
68  === The functions (in alphabetical order)
69  ===========================================================================*/
70 
71 /* Workspace method: Doxygen documentation will be auto-generated */
75  const Verbosity&) {
76  const Index np_total = TotalNumberOfElements(scat_meta);
77 
78  particle_masses.resize(np_total, 1);
79 
80  Index i_se_flat = 0;
81  for (Index i_ss = 0; i_ss < scat_meta.nelem(); i_ss++) {
82  for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
83  if (std::isnan(scat_meta[i_ss][i_se].mass) ||
84  scat_meta[i_ss][i_se].mass <= 0 || scat_meta[i_ss][i_se].mass > 1.) {
85  ostringstream os;
86  os << "A presumably incorrect value found for "
87  << "scat_meta[" << i_ss << "][" << i_se << "].mass.\n"
88  << "The value is " << scat_meta[i_ss][i_se].mass;
89  throw std::runtime_error(os.str());
90  }
91 
92  if (scat_meta[i_ss][i_se].diameter_volume_equ <= 0 ||
93  scat_meta[i_ss][i_se].diameter_volume_equ > 0.5) {
94  ostringstream os;
95  os << "A presumably incorrect value found for "
96  << "scat_meta[" << i_ss << "][" << i_se << "].diameter_volume_equ.\n"
97  << "The value is " << scat_meta[i_ss][i_se].diameter_volume_equ;
98  throw std::runtime_error(os.str());
99  }
100 
101  particle_masses(i_se_flat, 0) = scat_meta[i_ss][i_se].diameter_volume_equ;
102 
103  i_se_flat++;
104  }
105  }
106 }
107 
108 /* Workspace method: Doxygen documentation will be auto-generated */
109 void particle_massesFromMetaData( //WS Output:
111  // WS Input:
113  const Verbosity&) {
114  // resize particle_masses to required diemsions and properly initialize values
116  particle_masses = 0.;
117 
118  // calculate and set particle_masses
119  Index i_se_flat = 0;
120  for (Index i_ss = 0; i_ss < scat_meta.nelem(); i_ss++) {
121  for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
122  particle_masses(i_se_flat, i_ss) =
123  scat_meta[i_ss][i_se].diameter_volume_equ;
124  i_se_flat++;
125  }
126  }
127 }
128 
129 /* Workspace method: Doxygen documentation will be auto-generated */
132  const Vector& pnd_size_grid,
133  const Matrix& psd_data,
134  const Vector& psd_size_grid,
135  const Tensor3& dpsd_data_dx,
136  const Index& quad_order,
137  const Verbosity&) {
138  // Some sizes
139  const Index np = psd_data.nrows();
140  const Index ng = psd_size_grid.nelem();
141  Index ndx = 0;
142  const bool do_dx = !dpsd_data_dx.empty();
143 
144  // Checks
145  if (ng < 2)
146  throw runtime_error(
147  "The method requires that length of *psd_size_grid* is >= 2.");
148  if (ng != pnd_size_grid.nelem())
149  throw runtime_error(
150  "So far, the method requires that *psd_size_grid* and "
151  "*pnd_size_grid* have same length.");
152  for (Index i = 0; i < ng; i++) {
153  if (psd_size_grid[i] != pnd_size_grid[i])
154  throw runtime_error(
155  "So far, the method requires that *psd_size_grid* and "
156  "*pnd_size_grid* are identical.");
157  }
158  if (psd_data.ncols() != ng)
159  throw runtime_error(
160  "Number of columns in *psd_data* and length of "
161  "*psd_size_grid* must match.");
162 
163  pnd_data.resize(np, ng);
164  if (do_dx) {
165  if (dpsd_data_dx.ncols() != ng)
166  throw runtime_error(
167  "Number of columns in *dpsd_data_dx* and length of "
168  "*psd_size_grid* must match.");
169  ndx = dpsd_data_dx.npages();
170  dpnd_data_dx.resize(ndx, np, ng);
171  } else {
172  dpnd_data_dx.resize(0, 0, 0);
173  }
174 
175  // Get sorted version of psd_size_grid (and, since pnd_size_grid so far is
176  // identical, of this as well implicitly)
177  ArrayOfIndex intarr;
178  Vector psd_size_grid_sorted(ng);
180  for (Index i = 0; i < ng; i++)
181  psd_size_grid_sorted[i] = psd_size_grid[intarr[i]];
182 
183  if (!is_increasing(psd_size_grid_sorted))
184  throw std::runtime_error(
185  "*psd_size_grid* is not allowed to contain "
186  "duplicate values.");
187 
188  // Calculate pnd by intrgation of psd for given nodes/bins
189  Vector quadweights(ng);
190  bin_quadweights(quadweights, psd_size_grid_sorted, quad_order);
191 
192  for (Index i = 0; i < ng; i++) {
193  for (Index ip = 0; ip < np; ip++) {
194  pnd_data(ip, intarr[i]) = quadweights[i] * psd_data(ip, intarr[i]);
195  }
196 
197  if (do_dx) {
198  for (Index ip = 0; ip < np; ip++) {
199  for (Index ix = 0; ix < ndx; ix++) {
200  dpnd_data_dx(ix, ip, intarr[i]) =
201  quadweights[i] * dpsd_data_dx(ix, ip, intarr[i]);
202  }
203  }
204  }
205  }
206 }
207 
208 /* Workspace method: Doxygen documentation will be auto-generated */
211  const Vector& pnd_size_grid,
212  const Matrix& psd_data,
213  const Vector& psd_size_grid,
214  const Tensor3& dpsd_data_dx,
216  const Vector& f_grid,
217  const Index& scat_data_checked,
218  const Index& quad_order,
219  const Index& scat_index,
220  const Numeric& threshold_rsec,
221  const Numeric& threshold_bext,
222  const Numeric& threshold_rpnd,
223  const Verbosity&) {
224  // Some sizes
225  const Index np = psd_data.nrows();
226  const Index ng = psd_size_grid.nelem();
227  Index ndx = 0;
228  const bool do_dx = !dpsd_data_dx.empty();
229 
230  // Checks
231  if (ng < 3)
232  throw runtime_error(
233  "The method requires that length of *psd_size_grid* is >= 3.");
234  if (ng != pnd_size_grid.nelem())
235  throw runtime_error(
236  "So far, the method requires that *psd_size_grid* and"
237  " *pnd_size_grid* have the same length.");
238  for (Index i = 0; i < ng; i++) {
239  if (psd_size_grid[i] != pnd_size_grid[i])
240  throw runtime_error(
241  "So far, the method requires that *psd_size_grid* and"
242  " *pnd_size_grid* are identical.");
243  }
244  if (psd_data.ncols() != ng)
245  throw runtime_error(
246  "Number of columns in *psd_data* and length of"
247  " *psd_size_grid* must match.");
248 
249  pnd_data.resize(np, ng);
250  if (do_dx) {
251  if (dpsd_data_dx.ncols() != ng)
252  throw runtime_error(
253  "Number of columns in *dpsd_data_dx* and length of"
254  " *psd_size_grid* must match.");
255  ndx = dpsd_data_dx.npages();
256  dpnd_data_dx.resize(ndx, np, ng);
257  } else {
258  dpnd_data_dx.resize(0, 0, 0);
259  }
260 
261  if (!scat_data_checked)
262  throw runtime_error(
263  "*scat_data* must have passed a consistency check"
264  " (scat_data_checked=1).\n"
265  "Alternatively, use *pndFromPsdBasic*.");
266  if (scat_index >= scat_data.nelem())
267  throw runtime_error(
268  "*scat_index* exceeds the number of available"
269  " scattering species.");
270  if (scat_data[scat_index].nelem() != ng)
271  throw runtime_error(
272  "Number of scattering elements in this scattering"
273  " species (*scat_index*) inconsistent with length of"
274  " *pnd_size_grid*.");
275 
276  // Get sorted version of psd_size_grid (and, since pnd_size_grid so far is
277  // identical, of this as well implicitly)
278  ArrayOfIndex intarr;
279  Vector psd_size_grid_sorted(ng);
281  for (Index i = 0; i < ng; i++)
282  psd_size_grid_sorted[i] = psd_size_grid[intarr[i]];
283 
284  if (!is_increasing(psd_size_grid_sorted))
285  throw std::runtime_error(
286  "*psd_size_grid* is not allowed to contain "
287  "duplicate values.");
288 
289  // Calculate pnd by integration of psd for given nodes/bins
290  Vector quadweights(ng);
291  bin_quadweights(quadweights, psd_size_grid_sorted, quad_order);
292 
293  for (Index i = 0; i < ng;
294  i++) //loop over pnd_size_grid aka scattering elements
295  {
296  for (Index ip = 0; ip < np; ip++) //loop over pressure levels
297  {
298  pnd_data(ip, intarr[i]) = quadweights[i] * psd_data(ip, intarr[i]);
299  }
300 
301  if (do_dx) {
302  for (Index ip = 0; ip < np; ip++) {
303  for (Index ix = 0; ix < ndx; ix++) {
304  dpnd_data_dx(ix, ip, intarr[i]) =
305  quadweights[i] * dpsd_data_dx(ix, ip, intarr[i]);
306  }
307  }
308  }
309  }
310 
311  ArrayOfSingleScatteringData sds = scat_data[scat_index];
312  Index fstart = 0;
313  Index nf = f_grid.nelem();
314  Matrix bulkext(np, nf, 0.);
315  Vector ext(nf), ext_s0(nf), ext_s1(nf), ext_l0(nf), ext_l1(nf);
316 
317  // check that pnd_size_grid and scat_data cover (scalar) bulk extinction properly.
318  // (formerly we used mass. but heavy particles might not necessarily
319  // contribute much to optprops. same is valid for total particle number.)
320  //
321  // We do that by
322  // (a) checking that psd*ext decreases at the edges, i.e. to increase the
323  // probability that the main extinction peak is covered and in order to
324  // avoid passing check (b) through making the edge bins tiny
325  // (b) checking that the extinction contributions of the edge bins is small
326  // compared to the total bulk bin.
327  // The tests are somewhat over-sensitive on the small particle side in the
328  // sense that the checks are triggered easily by low mass density cases - which
329  // have their main extinction contributed by small particles, but for which
330  // the total extinction is quite low. Hence, there is also a total extinction
331  // threshold, below which the checks are not applied. Furthermore, such low
332  // densities do typically not occur throughout the whole atmospheric profile,
333  // but higher density layers might exists, meaning the low density layers
334  // contribute very little to the total optical properties, hence should not
335  // dominate the check criteria (no need to be very strict with them if they
336  // hardly contribute anyways). Therefore, checks are also not applied if the
337  // edge bin number density at a given atmospheric level is low compared to the
338  // maximum number of this scattering element in the profile.
339  //
340  // Technically, we need to cover all freqs separately (as optprops vary
341  // strongly with freq). We indeed do that here.
342  // FIXME: Shall a single freq or reduced freq-number option be implemented?
343  // If so, the don't-apply-test criteria need to be checked again though (so it
344  // does not happen that effectively all of the reduced-freq testing is skipped
345  // due to them).
346  //
347  // Technically, we also need to use identical temperatures (and preferably the
348  // ones used for deriving psd). But all scat elements can be on different
349  // T-grids, ie we would need to interpolate. Instead we just use the first in
350  // line.
351  // FIXME: Maybe refine in future allowing to select between low, high, median
352  // grid point or even specify one temp to use.
353  //
354  // And, technically ext might be directional dependent (for oriented
355  // particles). Also this might be different for each scat element. However, at
356  // least for now, we only check the 0-direction.
357  // FIXME: Check further directions?
358 
359  for (Index ise = 0; ise < ng;
360  ise++) //loop over pnd_size_grid aka scattering elements
361  {
362  if (sds[intarr[ise]].ext_mat_data.nshelves() > 1)
363  ext = sds[intarr[ise]].ext_mat_data(joker, 0, 0, 0, 0);
364  else
365  ext = sds[intarr[ise]].ext_mat_data(0, 0, 0, 0, 0);
366 
367  // keep the ext data of the edge bins
368  if (ise == 0)
369  ext_s0 = ext;
370  else if (ise == 1)
371  ext_s1 = ext;
372  else if (ise == ng - 2)
373  ext_l1 = ext;
374  else if (ise == ng - 1)
375  ext_l0 = ext;
376 
377  for (Index ip = 0; ip < np; ip++) //loop over pressure levels
378  if (abs(pnd_data(ip, joker).sum()) > 0.)
379  for (Index f = fstart; f < (fstart + nf); f++)
380  bulkext(ip, f) += pnd_data(ip, intarr[ise]) * ext[f];
381  }
382 
383  Numeric max0 = 0, max1 = 0;
384  for (Index ip = 0; ip < np; ip++) //loop over pressure levels
385  {
386  max0 = max(abs(pnd_data(ip, intarr[0])), max0);
387  max1 = max(abs(pnd_data(ip, intarr[ng - 1])), max1);
388  }
389 
390  Numeric contrib;
391  for (Index ip = 0; ip < np; ip++) //loop over pressure levels
392  {
393  if (abs(pnd_data(ip, joker).sum()) > 0.) {
394  for (Index f = fstart; f < (fstart + nf); f++) {
395  /* for( Index ise=0; ise<ng; ise++ )
396  {
397  if( sds[ise].ext_mat_data.nshelves()>1 )
398  ext = sds[ise].ext_mat_data(joker,0,0,0,0);
399  else
400  ext = sds[ise].ext_mat_data(0,0,0,0,0);
401  cout << " se #" << ise << ": contrib = pnd*ext/bext = "
402  << abs(pnd_data(ip,ise)) << "*" << ext[f] << "/"
403  << abs(bulkext(ip,f)) << " = "
404  << 1e2*abs(pnd_data(ip,ise))*ext[f]/abs(bulkext(ip,f))
405  << "%\n";
406  }*/
407 
408  // check that bin-width normalized extinction (or ext*psd) is
409  // decreasing.
410  if (abs(bulkext(ip, f)) > 1e-2 * threshold_bext) {
411  if (abs(psd_data(ip, intarr[0])) > 0. and
412  ext_s0[f] * abs(psd_data(ip, intarr[0])) >=
413  ext_s1[f] * abs(psd_data(ip, intarr[1]))) {
414  ostringstream os;
415  os << " Bin-width normalized extinction (ext*psd) not decreasing"
416  << " at small size edge\n"
417  << " at atm level #" << ip << " and freq point #" << f << ".\n"
418  << " ext_s0=" << ext_s0[f]
419  << ", psd_s0=" << abs(psd_data(ip, intarr[0]))
420  << ", ext_s0*psd_s0=" << ext_s0[f] * abs(psd_data(ip, intarr[0]))
421  << "\n LARGER EQUAL\n"
422  << " ext_s1=" << ext_s1[f]
423  << ", psd_s1=" << abs(psd_data(ip, intarr[1]))
424  << ", ext_s1*psd_s1=" << ext_s1[f] * abs(psd_data(ip, intarr[1]))
425  << "\n"
426  << " Total bulk ext = " << abs(bulkext(ip, f)) << "\n"
427  << " Need to add smaller sized particles!\n";
428  throw runtime_error(os.str());
429  }
430 
431  if (abs(psd_data(ip, intarr[ng - 1])) > 0. and
432  ext_l0[f] * abs(psd_data(ip, intarr[ng - 1])) >=
433  ext_l1[f] * abs(psd_data(ip, intarr[ng - 2]))) {
434  ostringstream os;
435  os << "Bin-width normalized extinction (ext*psd) not decreasing"
436  << " at large size edge\n"
437  << "at atm level #" << ip << " and freq point #" << f << ".\n"
438  << " ext_l0=" << ext_l0[f]
439  << ", psd_l0=" << abs(psd_data(ip, intarr[ng - 1]))
440  << ", ext_l0*psd_l0="
441  << ext_l0[f] * abs(psd_data(ip, intarr[ng - 1]))
442  << "\n LARGER EQUAL\n"
443  << " ext_l1=" << ext_l1[f]
444  << ", psd_l1=" << abs(psd_data(ip, intarr[ng - 2]))
445  << ", ext_l1*psd_l1="
446  << ext_l1[f] * abs(psd_data(ip, intarr[ng - 2])) << "\n"
447  << " Total bulk ext = " << abs(bulkext(ip, f)) << "\n"
448  << " Need to add larger sized particles!\n";
449  throw runtime_error(os.str());
450  }
451  }
452 
453  // check that contribution of edge bins to total extinction is
454  // sufficiently small
455  if (abs(bulkext(ip, f)) > threshold_bext) {
456  if (abs(pnd_data(ip, intarr[0])) > threshold_rpnd * max0) {
457  contrib =
458  abs(pnd_data(ip, intarr[0])) * ext_s0[f] / abs(bulkext(ip, f));
459  //cout << " small edge contrib = pnd*ext/bext = "
460  // << abs(pnd_data(ip,intarr[0])) << "*" << ext_s0[f] << "/"
461  // << abs(bulkext(ip,f)) << " = " << contrib << "\n";
462  if (abs(pnd_data(ip, intarr[0])) * ext_s0[f] >
463  threshold_rsec * abs(bulkext(ip, f))) {
464  ostringstream os;
465  os << "Contribution of edge bin to total extinction too high"
466  << " (" << contrib * 1e2 << "% of " << abs(bulkext(ip, f))
467  << ") at small size edge\n"
468  << "at atm level #" << ip << " and freq point #" << f << ".\n"
469  << " Need to add smaller sized particles or refine the size"
470  << " grid on the small size edge!\n";
471  throw runtime_error(os.str());
472  }
473  }
474  if (abs(pnd_data(ip, intarr[ng - 1])) > threshold_rpnd * max1) {
475  contrib = abs(pnd_data(ip, intarr[ng - 1])) * ext_l0[f] /
476  abs(bulkext(ip, f));
477  //cout << " large edge contrib = pnd*ext/bext = "
478  // << abs(pnd_data(ip,ng-1)) << "*" << ext_l0[f] << "/"
479  // << abs(bulkext(ip,f)) << " = " << contrib << "\n";
480  if (abs(pnd_data(ip, intarr[ng - 1])) * ext_l0[f] >
481  threshold_rsec * abs(bulkext(ip, f))) {
482  ostringstream os;
483  os << "Contribution of edge bin to total extinction too high"
484  << " (" << contrib * 1e2 << "% of " << abs(bulkext(ip, f))
485  << ") at large size edge\n"
486  << "at atm level #" << ip << " and freq point #" << f << ".\n"
487  << " Need to add larger sized particles or refine the size"
488  << " grid on the large size edge!\n";
489  throw runtime_error(os.str());
490  }
491  }
492  }
493  }
494  }
495  }
496 }
497 
498 /* Workspace method: Doxygen documentation will be auto-generated */
500  Workspace& ws,
503  const Index& atmosphere_dim,
504  const Vector& p_grid,
505  const Vector& lat_grid,
506  const Vector& lon_grid,
507  const Tensor3& t_field,
508  const Index& cloudbox_on,
517  const Index& jacobian_do,
519  const Verbosity&) {
520  // Do nothing if cloudbix ix inactive
521  if (!cloudbox_on) {
522  return;
523  }
524 
525  if (particle_bulkprop_field.empty())
526  throw runtime_error("*particle_bulkprop_field* is empty.");
527 
528  // Checks (not totally complete, but should cover most mistakes)
529  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
532  chk_atm_field("particle_bulkprop_field",
535  particle_bulkprop_names.nelem(),
536  p_grid,
537  lat_grid,
538  lon_grid);
539 
540  // Number of scattering species
541  const Index nss = scat_data.nelem();
542 
543  if (particle_bulkprop_names.nelem() == 0 ||
544  particle_bulkprop_names.nelem() != particle_bulkprop_field.nbooks()) {
545  throw runtime_error(
546  "Number of fields in *particle_bulkprop_field*"
547  " inconsistent with number of names in"
548  " *particle_bulkprop_names*.");
549  }
550  if (particle_bulkprop_field.nbooks() < nss) {
551  throw runtime_error(
552  "At least one field per scattering species required"
553  " in *particle_bulkprop_field*.");
554  }
555 
556  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
557  throw runtime_error(
558  "Length of *cloudbox_limits* incorrect with respect "
559  "to *atmosphere_dim*.");
560  if (cloudbox_limits[1] <= cloudbox_limits[0] || cloudbox_limits[0] < 0 ||
561  cloudbox_limits[1] >= p_grid.nelem())
562  throw runtime_error("Invalid data in pressure part of *cloudbox_limits*.");
563  if (atmosphere_dim > 1) {
564  if (cloudbox_limits[3] <= cloudbox_limits[2] || cloudbox_limits[2] < 0 ||
565  cloudbox_limits[3] >= lat_grid.nelem())
566  throw runtime_error(
567  "Invalid data in latitude part of *cloudbox_limits*.");
568  if (atmosphere_dim > 2) {
569  if (cloudbox_limits[5] <= cloudbox_limits[4] || cloudbox_limits[4] < 0 ||
570  cloudbox_limits[5] >= lon_grid.nelem())
571  throw runtime_error(
572  "Invalid data in longitude part of *cloudbox_limits*.");
573  }
574  }
575 
576  if (nss < 1) throw runtime_error("*scat_data* is empty!.");
577  if (scat_species.nelem() != nss)
578  throw runtime_error(
579  "*scat_data* and *scat_species* are inconsistent in size.");
580  if (scat_meta.nelem() != nss)
581  throw runtime_error(
582  "*scat_data* and *scat_meta* are inconsistent in size.");
583  if (pnd_agenda_array.nelem() != nss)
584  throw runtime_error(
585  "*scat_data* and *pnd_agenda_array* are inconsistent "
586  "in size.");
587  if (pnd_agenda_array_input_names.nelem() != nss)
588  throw runtime_error(
589  "*scat_data* and *pnd_agenda_array_input_names* are "
590  "inconsistent in size.");
591  // Further checks of scat_data vs. scat_meta below
592 
593  // Effective lengths of cloudbox
594  const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
595  const Index ip_offset = cloudbox_limits[0];
596  Index nlat = 1;
597  Index ilat_offset = 0;
598  if (atmosphere_dim > 1) {
599  nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
600  ilat_offset = cloudbox_limits[2];
601  }
602  Index nlon = 1;
603  Index ilon_offset = 0;
604  if (atmosphere_dim > 2) {
605  nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
606  ilon_offset = cloudbox_limits[4];
607  }
608 
609  // Check that *particle_bulkprop_field* contains zeros outside and at
610  // cloudbox boundaries
611  const String estring =
612  "*particle_bulkprop_field* allowed to contain"
613  " non-zero values only inside the cloudbox.";
614  // Pressure end ranges
615  Index dp_start, dp_end;
616  if (cloudbox_limits[0] != 0) {
618  joker, Range(0, ip_offset + 1), joker, joker)) > 0 ||
620  joker, Range(0, ip_offset + 1), joker, joker)) < 0)
621  throw runtime_error(estring);
622  dp_start = 1;
623  } else
624  dp_start = 0;
625  if (cloudbox_limits[1] != p_grid.nelem() - 1) {
626  const Index np_above = p_grid.nelem() + 1 - (np + ip_offset);
628  joker, Range(cloudbox_limits[1], np_above), joker, joker)) > 0 ||
630  joker, Range(cloudbox_limits[1], np_above), joker, joker)) < 0)
631  throw runtime_error(estring);
632  dp_end = np - 1;
633  } else
634  dp_end = np;
635 
636  // Cumulative number of scattering elements
637  ArrayOfIndex ncumse(nss + 1);
638  ncumse[0] = 0;
639  for (Index i = 0; i < nss; i++) {
640  if (scat_data[i].nelem() != scat_meta[i].nelem())
641  throw runtime_error(
642  "*scat_data* and *scat_meta* have inconsistent sizes.");
643  ncumse[i + 1] = ncumse[i] + scat_data[i].nelem();
644  }
645 
646  // Allocate output variables
647  //
648  pnd_field.resize(ncumse[nss], np, nlat, nlon);
649  pnd_field = 0.0; // To set all end values to zero
650  //
651  // Help variables for partial derivatives
652  Index nq = 0;
653  ArrayOfArrayOfIndex scatspecies_to_jq;
654  //
655  if (!jacobian_do) {
656  dpnd_field_dx.resize(0);
657  } else {
658  nq = jacobian_quantities.nelem();
659  dpnd_field_dx.resize(nq);
660  scatspecies_to_jq.resize(nss);
661  //
662  for (Index iq = 0; iq < nq; iq++) {
663  if (jacobian_quantities[iq].MainTag() == SCATSPECIES_MAINTAG) {
664  const Index ihit =
666  if (ihit < 0) {
667  ostringstream os;
668  os << "Jacobian quantity with index " << iq << " refers to\n"
669  << " " << jacobian_quantities[iq].Subtag()
670  << "\nbut this species could not be found in *scat_species*.";
671  throw runtime_error(os.str());
672  }
673  scatspecies_to_jq[ihit].push_back(iq);
674  dpnd_field_dx[iq].resize(ncumse[nss], np, nlat, nlon);
675  dpnd_field_dx[iq] = 0.0; // To set all end values to zero
676  }
677  }
678  }
679 
680  // Extract data from pnd-agenda array
681  for (Index is = 0; is < nss; is++) {
682  // Index range with respect to pnd_field
683  Range se_range(ncumse[is], ncumse[is + 1] - ncumse[is]);
684 
685  // Determine how pnd_agenda_array_input_names are related to input fields
686  //
687  const Index nin = pnd_agenda_array_input_names[is].nelem();
688  ArrayOfIndex i_pbulkprop(nin);
689  //
690  for (Index i = 0; i < nin; i++) {
691  i_pbulkprop[i] = find_first(particle_bulkprop_names,
693  if (i_pbulkprop[i] < 0) {
694  ostringstream os;
695  os << "Pnd-agenda with index " << is << " is set to require \""
696  << pnd_agenda_array_input_names[is][i] << "\",\nbut this quantity "
697  << "could not found in *particle_bulkprop_names*.\n"
698  << "(Note that temperature must be written as \"Temperature\")";
699  throw runtime_error(os.str());
700  }
701  }
702 
703  // Set *dpnd_data_dx_names*
704  //
705  Index ndx = 0;
707  //
708  if (jacobian_do) {
709  ndx = scatspecies_to_jq[is].nelem();
710  dpnd_data_dx_names.resize(ndx);
711  for (Index ix = 0; ix < ndx; ix++) {
712  dpnd_data_dx_names[ix] =
713  jacobian_quantities[scatspecies_to_jq[is][ix]].SubSubtag();
714  }
715  }
716 
717  // Loop lat/lon positions and call *pnd_agenda*
718  for (Index ilon = 0; ilon < nlon; ilon++) {
719  for (Index ilat = 0; ilat < nlat; ilat++) {
720  // Note that we don't need any calculations for end points
721 
722  // Here we consider this for lat and lon
723  if ((nlat > 1 && (ilat == 0 || ilat == nlat - 1)) ||
724  (nlon > 1 && (ilon == 0 || ilon == nlon - 1))) {
725  continue;
726  }
727 
728  // Pressure handled here, by not including end points in loops
729 
730  Matrix pnd_agenda_input(np, nin);
731  //
732  for (Index i = 0; i < nin; i++) {
733  for (Index ip = 0; ip < np; ip++) {
734  pnd_agenda_input(ip, i) =
735  particle_bulkprop_field(i_pbulkprop[i],
736  ip_offset + ip,
737  ilat_offset + ilat,
738  ilon_offset + ilon);
739  }
740  }
741 
743  //
744  for (Index ip = 0; ip < np; ip++) {
745  pnd_agenda_input_t[ip] =
746  t_field(ip_offset + ip, ilat_offset + ilat, ilon_offset + ilon);
747  }
748 
749  // Call pnd-agenda array
752  //
754  pnd_data,
755  dpnd_data_dx,
756  is,
762 
763  // Copy to output variables
764  for (Index ip = 0; ip < np; ip++) {
765  pnd_field(se_range, ip, ilat, ilon) = pnd_data(ip, joker);
766  }
767  for (Index ix = 0; ix < ndx; ix++) {
768  for (Index ip = dp_start; ip < dp_end; ip++) {
769  dpnd_field_dx[scatspecies_to_jq[is][ix]](se_range, ip, ilat, ilon) =
770  dpnd_data_dx(ix, ip, joker);
771  }
772  }
773  }
774  }
775  }
776 }
777 
778 /* Workspace method: Doxygen documentation will be auto-generated */
783  const Index& species_index,
784  const String& x_unit,
785  const Numeric& x_fit_start,
786  const Numeric& x_fit_end,
787  const Index& do_only_x,
788  const Verbosity&) {
789  // Checks
790  const Index nss = scat_meta.nelem();
791  if (nss == 0) throw runtime_error("*scat_meta* is empty!");
792  if (nss < species_index + 1) {
793  ostringstream os;
794  os << "Selected scattering species index is " << species_index
795  << " but this "
796  << "is not allowed since *scat_meta* has only " << nss << " elements.";
797  throw runtime_error(os.str());
798  }
799  //
800  const Index nse = scat_meta[species_index].nelem();
801  if (nse < 2)
802  throw runtime_error(
803  "The scattering species must have at least two "
804  "elements to use this method.");
805 
806  // Extract particle masses
807  //
808  Vector mass(nse);
809  //
810  for (Index i = 0; i < nse; i++) {
811  mass[i] = scat_meta[species_index][i].mass;
812  }
813 
814  // Create size grid
815  //
816  scat_species_x.resize(nse);
817  //
818  if (x_unit == "dveq") {
819  for (Index i = 0; i < nse; i++) {
820  scat_species_x[i] = scat_meta[species_index][i].diameter_volume_equ;
821  }
822  //
823  if (do_only_x) {
824  scat_species_a = -1;
825  scat_species_b = -1;
826  } else
830  mass,
831  x_fit_start,
832  x_fit_end);
833  }
834 
835  else if (x_unit == "dmax") {
836  for (Index i = 0; i < nse; i++) {
837  scat_species_x[i] = scat_meta[species_index][i].diameter_max;
838  }
839  //
840  if (do_only_x) {
841  scat_species_a = -1;
842  scat_species_b = -1;
843  } else
847  mass,
848  x_fit_start,
849  x_fit_end);
850  }
851 
852  else if (x_unit == "area") {
853  for (Index i = 0; i < nse; i++) {
854  scat_species_x[i] =
855  scat_meta[species_index][i].diameter_area_equ_aerodynamical;
856  }
857  //
858  if (do_only_x) {
859  scat_species_a = -1;
860  scat_species_b = -1;
861  } else
865  mass,
866  x_fit_start,
867  x_fit_end);
868  }
869 
870  else if (x_unit == "mass") {
871  scat_species_x = mass;
872  //
873  scat_species_a = 1;
874  scat_species_b = 1;
875  }
876 
877  else {
878  ostringstream os;
879  os << "You have selected the x_unit: " << x_unit
880  << "while accepted choices are: \"dveq\", \"dmax\", \"mass\" and \"area\"";
881  throw runtime_error(os.str());
882  }
883 }
884 
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::scat_meta
ArrayOfArrayOfScatteringMetaData scat_meta(Workspace &ws) noexcept
Definition: autoarts.h:5983
ARTS::Var::scat_species_a
Numeric scat_species_a(Workspace &ws) noexcept
Definition: autoarts.h:6103
auto_md.h
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
pndFromPsdBasic
void pndFromPsdBasic(Matrix &pnd_data, Tensor3 &dpnd_data_dx, const Vector &pnd_size_grid, const Matrix &psd_data, const Vector &psd_size_grid, const Tensor3 &dpsd_data_dx, const Index &quad_order, const Verbosity &)
WORKSPACE METHOD: pndFromPsdBasic.
Definition: m_microphysics.cc:130
SCATSPECIES_MAINTAG
const String SCATSPECIES_MAINTAG
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
interpolation.h
Header file for interpolation.cc.
ARTS::Var::pnd_field
Tensor4 pnd_field(Workspace &ws) noexcept
Definition: autoarts.h:5081
abs
#define abs(x)
Definition: legacy_continua.cc:20626
particle_massesFromMetaData
void particle_massesFromMetaData(Matrix &particle_masses, const ArrayOfArrayOfScatteringMetaData &scat_meta, const Verbosity &)
WORKSPACE METHOD: particle_massesFromMetaData.
Definition: m_microphysics.cc:109
ARTS::Var::psd_size_grid
Vector psd_size_grid(Workspace &ws) noexcept
Definition: autoarts.h:5460
bin_quadweights
void bin_quadweights(Vector &w, const Vector &x, const Index &order)
Definition: cloudbox.cc:611
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
ARTS::Var::jacobian_do
Index jacobian_do(Workspace &ws) noexcept
Definition: autoarts.h:3912
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
array.h
This file contains the definition of Array.
derive_scat_species_a_and_b
void derive_scat_species_a_and_b(Numeric &a, Numeric &b, const Vector &x, const Vector &mass, const Numeric &x_fit_start, const Numeric &x_fit_end)
Definition: microphysics.cc:77
ARTS::Var::dpnd_data_dx_names
ArrayOfString dpnd_data_dx_names(Workspace &ws) noexcept
Definition: autoarts.h:3235
ARTS::Var::particle_bulkprop_names
ArrayOfString particle_bulkprop_names(Workspace &ws) noexcept
Definition: autoarts.h:4823
microphysics.h
Internal functions for microphysics calculations (size distributions etc.)
Array
This can be used to make arrays out of anything.
Definition: array.h:108
ARTS::Var::dpsd_data_dx
Tensor3 dpsd_data_dx(Workspace &ws) noexcept
Definition: autoarts.h:3281
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
get_sorted_indexes
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:73
min
#define min(a, b)
Definition: legacy_continua.cc:20628
ARTS::Var::pnd_data
Matrix pnd_data(Workspace &ws) noexcept
Definition: autoarts.h:5051
particle_massesFromMetaDataSingleCategory
void particle_massesFromMetaDataSingleCategory(Matrix &particle_masses, const ArrayOfArrayOfScatteringMetaData &scat_meta, const Verbosity &)
WORKSPACE METHOD: particle_massesFromMetaDataSingleCategory.
Definition: m_microphysics.cc:72
messages.h
Declarations having to do with the four output streams.
ARTS::Var::pnd_agenda_array
ArrayOfAgenda pnd_agenda_array(Workspace &ws) noexcept
Definition: autoarts.h:4979
ARTS::Var::particle_masses
Matrix particle_masses(Workspace &ws) noexcept
Definition: autoarts.h:4850
my_basic_string< char >
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
optproperties.h
Scattering database structure and functions.
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:49
math_funcs.h
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
ARTS::Var::scat_species_x
Vector scat_species_x(Workspace &ws) noexcept
Definition: autoarts.h:6130
max
#define max(a, b)
Definition: legacy_continua.cc:20629
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
pnd_fieldCalcFromParticleBulkProps
void pnd_fieldCalcFromParticleBulkProps(Workspace &ws, Tensor4 &pnd_field, ArrayOfTensor4 &dpnd_field_dx, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const ArrayOfArrayOfScatteringMetaData &scat_meta, const Tensor4 &particle_bulkprop_field, const ArrayOfString &particle_bulkprop_names, const ArrayOfAgenda &pnd_agenda_array, const ArrayOfArrayOfString &pnd_agenda_array_input_names, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: pnd_fieldCalcFromParticleBulkProps.
Definition: m_microphysics.cc:499
parameters.h
This file contains header information for the dealing with command line parameters.
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
ARTS::Var::dpnd_field_dx
ArrayOfTensor4 dpnd_field_dx(Workspace &ws) noexcept
Definition: autoarts.h:3251
cloudbox.h
Internal cloudbox functions.
psd.h
Internal functions associated with size distributions.
ARTS::Var::scat_species_b
Numeric scat_species_b(Workspace &ws) noexcept
Definition: autoarts.h:6112
ARTS::Var::scat_data_checked
Index scat_data_checked(Workspace &ws) noexcept
Definition: autoarts.h:5850
Range
The range class.
Definition: matpackI.h:160
ScatSpeciesSizeMassInfo
void ScatSpeciesSizeMassInfo(Vector &scat_species_x, Numeric &scat_species_a, Numeric &scat_species_b, const ArrayOfArrayOfScatteringMetaData &scat_meta, const Index &species_index, const String &x_unit, const Numeric &x_fit_start, const Numeric &x_fit_end, const Index &do_only_x, const Verbosity &)
WORKSPACE METHOD: ScatSpeciesSizeMassInfo.
Definition: m_microphysics.cc:779
pndFromPsd
void pndFromPsd(Matrix &pnd_data, Tensor3 &dpnd_data_dx, const Vector &pnd_size_grid, const Matrix &psd_data, const Vector &psd_size_grid, const Tensor3 &dpsd_data_dx, const ArrayOfArrayOfSingleScatteringData &scat_data, const Vector &f_grid, const Index &scat_data_checked, const Index &quad_order, const Index &scat_index, const Numeric &threshold_rsec, const Numeric &threshold_bext, const Numeric &threshold_rpnd, const Verbosity &)
WORKSPACE METHOD: pndFromPsd.
Definition: m_microphysics.cc:209
chk_atm_field
void chk_atm_field(const String &x_name, ConstTensor3View x, const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, const bool &chk_lat90)
chk_atm_field (simple fields)
Definition: check_input.cc:1268
is_increasing
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
logic.h
Header file for logic.cc.
chk_atm_grids
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
Definition: check_input.cc:1201
ARTS::Var::particle_bulkprop_field
Tensor4 particle_bulkprop_field(Workspace &ws) noexcept
Definition: autoarts.h:4810
rte.h
Declaration of functions in rte.cc.
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
special_interp.h
Header file for special_interp.cc.
ARTS::Var::dpnd_data_dx
Tensor3 dpnd_data_dx(Workspace &ws) noexcept
Definition: autoarts.h:3223
ARTS::Var::pnd_agenda_input_t
Vector pnd_agenda_input_t(Workspace &ws) noexcept
Definition: autoarts.h:5038
pnd_agenda_arrayExecute
void pnd_agenda_arrayExecute(Workspace &ws, Matrix &pnd_data, Tensor3 &dpnd_data_dx, const Index agenda_array_index, const Vector &pnd_agenda_input_t, const Matrix &pnd_agenda_input, const ArrayOfString &pnd_agenda_input_names, const ArrayOfString &dpnd_data_dx_names, const ArrayOfAgenda &input_agenda_array)
Definition: auto_md.cc:24714
file.h
This file contains basic functions to handle ASCII files.
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
ARTS::Var::scat_species
ArrayOfString scat_species(Workspace &ws) noexcept
Definition: autoarts.h:6086
find_first
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:292
check_input.h
Vector
The Vector class.
Definition: matpackI.h:860
sorting.h
Contains sorting routines.
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Var::pnd_agenda_input
Matrix pnd_agenda_input(Workspace &ws) noexcept
Definition: autoarts.h:5011
ARTS::Var::pnd_size_grid
Vector pnd_size_grid(Workspace &ws) noexcept
Definition: autoarts.h:5119
ARTS::Var::pnd_agenda_array_input_names
ArrayOfArrayOfString pnd_agenda_array_input_names(Workspace &ws) noexcept
Definition: autoarts.h:4991
ARTS::Var::psd_data
Matrix psd_data(Workspace &ws) noexcept
Definition: autoarts.h:5446
arts.h
The global header file for ARTS.
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
xml_io.h
This file contains basic functions to handle XML data files.