ARTS  2.2.66
disort.cc
Go to the documentation of this file.
1 /* Copyright (C) 2006-2012 Claudia Emde <claudia.emde@dlr.de>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA.
17 */
18 
28 #include <stdexcept>
29 #include <iostream>
30 #include <cstdlib>
31 #include <cmath>
32 #include "array.h"
33 #include "agenda_class.h"
34 #include "messages.h"
35 #include "xml_io.h"
36 #include "logic.h"
37 #include "check_input.h"
38 #include "auto_md.h"
39 #include "interpolation.h"
40 
41 extern const Numeric PI;
42 extern const Numeric PLANCK_CONST;
43 extern const Numeric SPEED_OF_LIGHT;
44 extern const Numeric BOLTZMAN_CONST;
45 
47 
69  VectorView dtauc,
70  VectorView ssalb,
71  const Agenda& opt_prop_part_agenda,
72  const Agenda& propmat_clearsky_agenda,
73  const Agenda& spt_calc_agenda,
74  ConstTensor4View pnd_field,
75  ConstTensor3View t_field,
76  ConstTensor3View z_field,
77  ConstVectorView p_grid,
78  ConstTensor4View vmr_field,
79  ConstVectorView f_mono
80  )
81 {
82 
83  const Index N_pt = pnd_field.nbooks();
84  // In DISORT the "cloudbox" must cover the whole atmosphere
85  const Index Np_cloud = pnd_field.npages();
86  const Index stokes_dim = 1;
87 
88  assert( dtauc.nelem() == Np_cloud-1);
89  assert( ssalb.nelem() == Np_cloud-1);
90 
91  // Local variables to be used in agendas
92  Matrix abs_vec_spt_local(N_pt, stokes_dim, 0.);
93  Tensor3 ext_mat_spt_local(N_pt, stokes_dim, stokes_dim, 0.);
94  Matrix abs_vec_local;
95  Tensor3 ext_mat_local;
96  Numeric rtp_temperature_local;
97  Numeric rtp_pressure_local;
98  Tensor4 propmat_clearsky_local;
99  Vector ext_vector(Np_cloud);
100  Vector abs_vector(Np_cloud);
101  Vector rtp_vmr_local(vmr_field.nbooks());
102  // Calculate ext_mat, abs_vec and sca_vec for all pressure points.
103 
104  propmat_clearsky_local = 0.;
105 
106 
107  for(Index scat_p_index_local = 0; scat_p_index_local < Np_cloud;
108  scat_p_index_local ++)
109  {
110  rtp_temperature_local =
111  t_field(scat_p_index_local, 0, 0);
112 
113  //Calculate optical properties for single particle types:
115  ext_mat_spt_local,
116  abs_vec_spt_local,
117  scat_p_index_local, 0, 0, //position
118  rtp_temperature_local,
119  0, 0, // angles, only needed for za=0
120  spt_calc_agenda);
121 
123  ext_mat_local, abs_vec_local,
124  ext_mat_spt_local,
125  abs_vec_spt_local,
126  scat_p_index_local, 0, 0,
127  opt_prop_part_agenda);
128 
129  ext_vector[scat_p_index_local] = ext_mat_local(0,0,0);
130  abs_vector[scat_p_index_local] = abs_vec_local(0,0);
131  }
132 
133 
134 
135  // Calculate layer averaged single scattering albedo and optical depth
136  for (Index i = 0; i < Np_cloud-1; i++)
137  {
138  Numeric ext = 0.;
139  Numeric abs = 0.;
140 
141  ext=.5*(ext_vector[i]+ext_vector[i+1]);
142  abs=.5*(abs_vector[i]+abs_vector[i+1]);
143 
144  if (ext!=0)
145  ssalb[Np_cloud-2-i]=(ext-abs)/ext;
146 
147  rtp_pressure_local = 0.5 * (p_grid[i] + p_grid[i+1]);
148  rtp_temperature_local = 0.5 * (t_field(i,0,0) + t_field(i+1,0,0));
149 
150  // Average vmrs
151  for (Index j = 0; j < vmr_field.nbooks(); j++)
152  rtp_vmr_local[j] = 0.5 * (vmr_field(j, i, 0, 0) +
153  vmr_field(j, i+1, 0, 0));
154 
155  const Vector rtp_mag_dummy(3,0);
156  const Vector ppath_los_dummy;
157 
159  propmat_clearsky_local,
160  f_mono, // monochromatic calculation
161  rtp_mag_dummy,ppath_los_dummy,
162  rtp_pressure_local,
163  rtp_temperature_local,
164  rtp_vmr_local,
165  propmat_clearsky_agenda);
166 
167  Numeric abs_total = propmat_clearsky_local(joker,0,0,0).sum(); //Assuming non-polarized light and only one frequency
168 
169  dtauc[Np_cloud-2-i]=(ext+abs+abs_total)*
170  (z_field(i+1, 0, 0)-z_field(i, 0, 0));
171  }
172 }
173 
175 
189 void phase_functionCalc(//Output
190  MatrixView phase_function,
191  //Input
192  const ArrayOfSingleScatteringData& scat_data_array_mono,
193  ConstTensor4View pnd_field)
194 {
195  Matrix phase_function_level(pnd_field.npages(),
196  scat_data_array_mono[0].za_grid.nelem(), 0.);
197 
198  //Loop over pressure levels
199  for (Index i_p = 0; i_p < pnd_field.npages(); i_p++)
200  {
201  // Loop over scattering angles
202  for (Index i_t = 0; i_t < scat_data_array_mono[0].za_grid.nelem(); i_t++)
203  {
204  // Calculate ensemble averaged extinction coefficient
205  Numeric sca_coeff=0.;
206 
207  for (Index j = 0; j < scat_data_array_mono.nelem(); j++)
208  sca_coeff += pnd_field(j, i_p, 0, 0) *
209  (scat_data_array_mono[j].ext_mat_data(0, 0, 0, 0, 0)-
210  scat_data_array_mono[j].abs_vec_data(0, 0, 0, 0, 0));
211 
212  // Phase function
213  for (Index j = 0; j < scat_data_array_mono.nelem(); j++)
214  {
215  if (sca_coeff != 0)
216  phase_function_level(i_p, i_t) +=
217  pnd_field(j, i_p, 0, 0) *
218  scat_data_array_mono[j].pha_mat_data(0, 0, i_t, 0, 0, 0, 0)
219  *4*PI/sca_coeff;// Normalization
220  }
221 
222  }
223  }
224 
225 
226  // Calculate average phase function for the layers
227  for (Index i_l = 0; i_l < pnd_field.npages()-1; i_l++)
228  {
229  for (Index i_t=0; i_t < phase_function_level.ncols(); i_t++)
230  {
231  if (phase_function_level(i_l, i_t) !=0 &&
232  phase_function_level(i_l+1, i_t) !=0)
233  phase_function(i_l, i_t) = .5*
234  (phase_function_level(i_l, i_t)+
235  phase_function_level(i_l+1, i_t));
236  }
237  }
238 
239 }
240 
242 
255 void pmomCalc(//Output
256  MatrixView pmom,
257  //Input
258  ConstMatrixView phase_function,
259  ConstVectorView scat_angle_grid,
260  const Index n_legendre,
261  const Verbosity& verbosity)
262 {
263  Numeric pint; //integrated phase function
264  Numeric p0_1, p0_2, p1_1, p1_2, p2_1, p2_2;
265 
266  Vector za_grid(181);
267  Vector u(181);
268 
269  for (Index i = 0; i< 181; i++)
270  za_grid[i] = double(i);
271 
272  ArrayOfGridPos gp(181);
273  gridpos(gp, scat_angle_grid, za_grid);
274 
275  Matrix itw(gp.nelem(),2);
276  interpweights(itw,gp);
277 
278  Matrix phase_int(phase_function.nrows(),181);
279  for (Index i_l=0; i_l < phase_function.nrows(); i_l++)
280  interp(phase_int(i_l, joker), itw, phase_function(i_l, joker), gp);
281 
282  for (Index i = 0; i<za_grid.nelem(); i++)
283  u[i] = cos(za_grid[i] *PI/180.);
284 
285  for (Index i_l=0; i_l < phase_function.nrows(); i_l++)
286  {
287  pint = 0.;
288  // Check if phase function is normalized
289  for (Index i = 0; i<za_grid.nelem()-1; i++)
290  pint+=0.5*(phase_int(i_l, i)+phase_int(i_l, i+1))*
291  abs(u[i+1] - u[i]);
292 
293  if (pint != 0){
294  if (abs(2.-pint) > 1e-4)
295  {
296  CREATE_OUT1;
297  out1 << "Warning: The phase function is not normalized to 2\n"
298  << "The value is:" << pint << "\n";
299  }
300 
301  pmom(i_l, joker)= 0.;
302 
303  for (Index i = 0; i<za_grid.nelem()-1; i++)
304  {
305  p0_1=1.;
306  p0_2=1.;
307 
308  pmom(phase_function.nrows()-1-i_l,0)=1.;
309 
310  //pmom(phase_function.nrows()-1-i_l,0)+=0.5*0.5*(phase_int(i_l, i)+
311  // phase_int(i_l, i+1))
312  //*abs(u[i+1]-u[i]);
313 
314  p1_1=u[i];
315  p1_2=u[i+1];
316 
317  pmom(phase_function.nrows()-1-i_l,1)+=0.5*0.5*
318  (p1_1*phase_int(i_l, i)+
319  p1_2*phase_int(i_l, i+1))
320  *abs(u[i+1]-u[i]);
321 
322  for (Index l=2; l<n_legendre; l++)
323  {
324  p2_1=(2*(double)l-1)/(double)l*u[i]*p1_1-((double)l-1)/
325  (double)l*p0_1;
326  p2_2=(2*(double)l-1)/(double)l*u[i+1]*p1_2-((double)l-1)/
327  (double)l*p0_2;
328 
329  pmom(phase_function.nrows()-1-i_l, l)+=0.5*0.5*
330  (p2_1*phase_int(i_l, i)+
331  p2_2*phase_int(i_l, i+1))
332  *abs(u[i+1]-u[i]);
333 
334  p0_1=p1_1;
335  p0_2=p1_2;
336  p1_1=p2_1;
337  p1_2=p2_2;
338  }
339  }
340  // cout << "pmom : " << pmom(phase_function.nrows()-1-i_l, joker) << endl;
341 
342  }
343  }
344 }
345 
347 
365  const Numeric& f,
366  const Numeric& t )
367 {
368  assert( f > 0 );
369  assert( t >= 0 );
370 
371  // Double must be used here (if not, a becomes 0 when using float)
372  static const double a = 2 * PLANCK_CONST / (SPEED_OF_LIGHT*SPEED_OF_LIGHT);
373  static const double b = PLANCK_CONST / BOLTZMAN_CONST;
374 
375  return a * f*f*f / ( exp( b*f/t ) - 1 );
376 }
377 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void opt_prop_part_agendaExecute(Workspace &ws, Tensor3 &ext_mat, Matrix &abs_vec, const Tensor3 &ext_mat_spt, const Matrix &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Agenda &input_agenda)
Definition: auto_md.cc:14459
The VectorView class.
Definition: matpackI.h:372
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
Declarations having to do with the four output streams.
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
The Tensor4 class.
Definition: matpackIV.h:383
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
This file contains basic functions to handle XML data files.
Header file for interpolation.cc.
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
const Numeric BOLTZMAN_CONST
A constant view of a Tensor4.
Definition: matpackIV.h:141
This file contains the definition of Array.
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
void pmomCalc(MatrixView pmom, ConstMatrixView phase_function, ConstVectorView scat_angle_grid, const Index n_legendre, const Verbosity &verbosity)
pmomCalc
Definition: disort.cc:255
Declarations for agendas.
Numeric planck2(const Numeric &f, const Numeric &t)
planck
Definition: disort.cc:364
void phase_functionCalc(MatrixView phase_function, const ArrayOfSingleScatteringData &scat_data_array_mono, ConstTensor4View pnd_field)
phase_functionCalc
Definition: disort.cc:189
#define CREATE_OUT1
Definition: messages.h:212
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void propmat_clearsky_agendaExecute(Workspace &ws, Tensor4 &propmat_clearsky, const Vector &f_grid, const Vector &rtp_mag, const Vector &rtp_los, const Numeric rtp_pressure, const Numeric rtp_temperature, const Vector &rtp_vmr, const Agenda &input_agenda)
Definition: auto_md.cc:13039
#define abs(x)
Definition: continua.cc:20458
const Joker joker
const Numeric SPEED_OF_LIGHT
const Numeric PI
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Header file for logic.cc.
This can be used to make arrays out of anything.
Definition: array.h:40
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
A constant view of a Matrix.
Definition: matpackI.h:596
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
void spt_calc_agendaExecute(Workspace &ws, Tensor3 &ext_mat_spt, Matrix &abs_vec_spt, const Index scat_p_index, const Index scat_lat_index, const Index scat_lon_index, const Numeric rtp_temperature, const Index scat_za_index, const Index scat_aa_index, const Agenda &input_agenda)
Definition: auto_md.cc:14960
const Numeric PLANCK_CONST
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void dtauc_ssalbCalc(Workspace &ws, VectorView dtauc, VectorView ssalb, const Agenda &opt_prop_part_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &spt_calc_agenda, ConstTensor4View pnd_field, ConstTensor3View t_field, ConstTensor3View z_field, ConstVectorView p_grid, ConstTensor4View vmr_field, ConstVectorView f_mono)
dtauc_ssalbCalc
Definition: disort.cc:68
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832