ARTS  2.2.66
m_batch.cc
Go to the documentation of this file.
1 /* Copyright (C) 2004-2012
2  Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3  Stefan Buehler <sbuehler@ltu.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
20 
21 
22 /*===========================================================================
23  === File description
24  ===========================================================================*/
25 
39 /*===========================================================================
40  === External declarations
41  ===========================================================================*/
42 
43 #include <cmath>
44 using namespace std;
45 
46 #include "arts.h"
47 #include "arts_omp.h"
48 #include "auto_md.h"
49 #include "math_funcs.h"
50 #include "physics_funcs.h"
51 #include "rte.h"
52 #include "xml_io.h"
53 
54 extern const Numeric PI;
55 extern const Numeric DEG2RAD;
56 extern const Numeric RAD2DEG;
57 extern const Index GFIELD3_P_GRID;
58 
59 
60 /*===========================================================================
61  === The functions (in alphabetical order)
62  ===========================================================================*/
63 
64 
65 /* Workspace method: Doxygen documentation will be auto-generated
66 
67  2008-07-21 Stefan Buehler */
68 void ForLoop(Workspace& ws,
69  // WS Input:
70  const Agenda& forloop_agenda,
71  // Control Parameters:
72  const Index& start,
73  const Index& stop,
74  const Index& step,
75  const Verbosity& verbosity)
76 {
78 
79  for (Index i=start; i<=stop; i+=step)
80  {
81  out1 << " Executing for loop body, index: " << i << "\n";
82  forloop_agendaExecute(ws, i, forloop_agenda);
83  }
84 }
85 
86 
87 /* Workspace method: Doxygen documentation will be auto-generated
88 
89  2004-09-15 Patrick Eriksson
90 
91  Added keyword to control if row or column is extracted.
92 
93  2007-07-24 Stefan Buehler */
95  // WS Generic Output:
96  Vector& v,
97  // WS Input:
98  // WS Generic Input:
99  const Matrix& m,
100  const Index& index,
101  // Control Parameters:
102  const String& direction,
103  const Verbosity&)
104 {
105  if (direction=="row")
106  {
107  if( index >= m.nrows() )
108  {
109  ostringstream os;
110  os << "The index " << index
111  << " is outside the row range of the Matrix.";
112  throw runtime_error( os.str() );
113 
114  }
115 
116  v.resize( m.ncols() );
117  v = m( index, joker );
118  }
119  else if (direction=="column")
120  {
121  if( index >= m.ncols() )
122  {
123  ostringstream os;
124  os << "The index " << index
125  << " is outside the column range of the Matrix.";
126  throw runtime_error( os.str() );
127 
128  }
129 
130  v.resize( m.nrows() );
131  v = m( joker, index );
132  }
133  else
134  {
135  ostringstream os;
136  os << "Keyword *direction* must be either *row* or *column*,"
137  << "but you gave: " << direction << ".";
138  throw runtime_error( os.str() );
139  }
140 }
141 
142 
143 /* Workspace method: Doxygen documentation will be auto-generated */
145  // WS Generic Output:
146  Matrix& m,
147  // WS Input:
148  // WS Generic Input:
149  const Tensor3& t3,
150  const Index& index,
151  // Control Parameters:
152  const String& direction,
153  const Verbosity&)
154 {
155  if (direction=="page")
156  {
157  if( index >= t3.npages() )
158  {
159  ostringstream os;
160  os << "The index " << index
161  << " is outside the page range of the Matrix.";
162  throw runtime_error( os.str() );
163 
164  }
165 
166  m.resize( t3.nrows(), t3.ncols() );
167  m = t3( index, joker, joker );
168  }
169  else if (direction=="row")
170  {
171  if( index >= t3.nrows() )
172  {
173  ostringstream os;
174  os << "The index " << index
175  << " is outside the row range of the Matrix.";
176  throw runtime_error( os.str() );
177 
178  }
179 
180  m.resize( t3.npages(), t3.ncols() );
181  m = t3( joker, index, joker );
182  }
183  else if (direction=="column")
184  {
185  if( index >= t3.ncols() )
186  {
187  ostringstream os;
188  os << "The index " << index
189  << " is outside the column range of the Matrix.";
190  throw runtime_error( os.str() );
191 
192  }
193 
194  m.resize( t3.npages(), t3.nrows() );
195  m = t3( joker, joker, index );
196  }
197  else
198  {
199  ostringstream os;
200  os << "Keyword *direction* must be either *page* or *row* or *column*,"
201  << "but you gave: " << direction << ".";
202  throw runtime_error( os.str() );
203  }
204 }
205 
206 
207 /* Workspace method: Doxygen documentation will be auto-generated */
209  // WS Output:
210  ArrayOfVector& ybatch,
211  ArrayOfArrayOfVector& ybatch_aux,
212  ArrayOfMatrix& ybatch_jacobians,
213  // WS Input:
214  const Index& ybatch_start,
215  const Index& ybatch_n,
216  const Agenda& ybatch_calc_agenda,
217  // Control Parameters:
218  const Index& robust,
219  const Verbosity& verbosity)
220 {
221  CREATE_OUTS;
222 
223  Index first_ybatch_index = 0;
224 
225  ArrayOfString fail_msg;
226  bool do_abort = false;
227 
228  // We allow a start index ybatch_start that is different from 0. We
229  // will calculate ybatch_n jobs starting at the start
230  // index. Internally, we count from zero, which is the right
231  // counting for the output array ybatch. When we call
232  // ybatch_calc_agenda, we add ybatch_start to the internal index
233  // count.
234 
235  // We create a counter, so that we can generate nice output about
236  // how many jobs are already done. (All parallel threads later will
237  // increment this, so that we really get an accurate total count!)
238  Index job_counter = 0;
239 
240  // Resize the output arrays:
241  ybatch.resize(ybatch_n);
242  ybatch_aux.resize(ybatch_n);
243  ybatch_jacobians.resize(ybatch_n);
244 
245  // We have to make a local copy of the Workspace and the agendas because
246  // only non-reference types can be declared firstprivate in OpenMP
247  Workspace l_ws(ws);
248  Agenda l_ybatch_calc_agenda(ybatch_calc_agenda);
249 
250  // Go through the batch:
251 
252  if (ybatch_n)
253 #pragma omp parallel for \
254  schedule(dynamic) \
255  if (!arts_omp_in_parallel() \
256  && ybatch_n > 1) \
257  firstprivate(l_ws, l_ybatch_calc_agenda)
258  for(Index ybatch_index = first_ybatch_index;
259  ybatch_index<ybatch_n;
260  ybatch_index++ )
261  {
262  Index l_job_counter; // Thread-local copy of job counter.
263 
264  if (do_abort) continue;
265 #pragma omp critical (ybatchCalc_job_counter)
266  {
267  l_job_counter = ++job_counter;
268  }
269 
270  {
271  ostringstream os;
272  os << " Job " << l_job_counter << " of " << ybatch_n
273  << ", Index " << ybatch_start+ybatch_index << ", Thread-Id "
274  << arts_omp_get_thread_num() << "\n";
275  out2 << os.str();
276  }
277 
278  try
279  {
280  Vector y;
281  ArrayOfVector y_aux;
282  Matrix jacobian;
283 
284  ybatch_calc_agendaExecute(l_ws, y, y_aux, jacobian,
285  ybatch_start+ybatch_index,
286  l_ybatch_calc_agenda);
287 
288  if (y.nelem())
289  {
290 #pragma omp critical (ybatchCalc_assign_y)
291  ybatch[ybatch_index] = y;
292 #pragma omp critical (ybatchCalc_assign_y_aux)
293  ybatch_aux[ybatch_index] = y_aux;
294 
295  // Dimensions of Jacobian:
296  const Index Knr = jacobian.nrows();
297  const Index Knc = jacobian.ncols();
298 
299  if (Knr != 0 || Knc != 0)
300  {
301  if (Knr != y.nelem())
302  {
303  ostringstream os;
304  os << "First dimension of Jacobian must have same length as the measurement *y*.\n"
305  << "Length of *y*: " << y.nelem() << "\n"
306  << "Dimensions of *jacobian*: (" << Knr << ", " << Knc << ")\n";
307  // A mismatch of the Jacobian dimension is a fatal error
308  // and should result in program termination. By setting abort
309  // to true, this will result in a runtime error in the catch
310  // block even if robust == 1
311 #pragma omp critical (ybatchCalc_setabort)
312  do_abort = true;
313 
314  throw runtime_error(os.str());
315  }
316 
317  ybatch_jacobians[ybatch_index] = jacobian;
318 
319  // After creation, all individual Jacobi matrices in the array will be
320  // empty (size zero). No need for explicit initialization.
321  }
322  }
323  }
324  catch (runtime_error e)
325  {
326  if (robust && !do_abort)
327  {
328  ostringstream os;
329  os << "WARNING! Job at ybatch_index " << ybatch_start+ybatch_index << " failed.\n"
330  << "y Vector in output variable ybatch will be empty for this job.\n"
331  << "The runtime error produced was:\n"
332  << e.what() << "\n";
333  out0 << os.str();
334  }
335  else
336  {
337  // The user wants the batch job to fail if one of the
338  // jobs goes wrong.
339 #pragma omp critical (ybatchCalc_setabort)
340  do_abort = true;
341 
342  ostringstream os;
343  os << " Job at ybatch_index "
344  << ybatch_start+ybatch_index << " failed. Aborting...\n";
345  out1 << os.str();
346  }
347  ostringstream os;
348  os << "Run-time error at ybatch_index "
349  << ybatch_start+ybatch_index << ": \n" << e.what();
350 #pragma omp critical (ybatchCalc_push_fail_msg)
351  fail_msg.push_back(os.str());
352  }
353  }
354 
355  if (fail_msg.nelem())
356  {
357  ostringstream os;
358 
359  if (!do_abort) os << "\nError messages from failed batch cases:\n";
360  for (ArrayOfString::const_iterator it = fail_msg.begin(); it != fail_msg.end(); it++)
361  os << *it << '\n';
362 
363  if (do_abort)
364  throw runtime_error(os.str());
365  else
366  out0 << os.str();
367  }
368 
369 }
370 
371 
372 /* Workspace method: Doxygen documentation will be auto-generated */
374  Workspace& ws,
375  //Output
376  ArrayOfVector& ybatch,
377  //Input
378  const ArrayOfArrayOfSpeciesTag& abs_species,
379  const Agenda& met_profile_calc_agenda,
380  const Vector& f_grid,
381  const Matrix& met_amsu_data,
382  const Matrix& sensor_pos,
383  const Vector& refellipsoid,
384  const Vector& lat_grid,
385  const Vector& lon_grid,
386  const Index& atmosphere_dim,
387  const ArrayOfSingleScatteringData& scat_data_array,
388  //Keyword
389  const Index& nelem_p_grid,
390  const String& met_profile_path,
391  const String& met_profile_pnd_path,
392  const Verbosity& verbosity)
393 {
394  GriddedField3 t_field_raw;
395  GriddedField3 z_field_raw;
396  ArrayOfGriddedField3 vmr_field_raw;
397  ArrayOfGriddedField3 pnd_field_raw;
398  Vector p_grid;
399  Matrix sensor_los;
400  Index cloudbox_on;
401  ArrayOfIndex cloudbox_limits;
402  Matrix z_surface;
403  Vector y;
404  Index no_profiles = met_amsu_data.nrows();
405 
406  // *vmr_field_raw* is an ArrayOfArrayOfTensor3 where the first array
407  //holds the gaseous species.
408  //Resize *vmr_field_raw* according to the number of gaseous species
409  //elements
410  vmr_field_raw.resize(abs_species.nelem());
411 
412  //pnd_field_raw is an ArrayOfArrayOfTensor3 where the first array
413  //holds particle species.
414  // Number of particle types:
415  const Index N_pt = scat_data_array.nelem();
416 
417  pnd_field_raw.resize(N_pt);
418 
419  // The satellite zenith angle is read in from the amsu data
420  // and converted to arts sensor_los
421  ConstVectorView sat_za_from_data = met_amsu_data(Range(joker),3);
422 
423  sensor_los.resize(1,1);
424 
425  // The lat and lon are extracted to get the proper file names of
426  // profiles
427  ConstVectorView lat = met_amsu_data(Range(joker),0);
428  ConstVectorView lon = met_amsu_data(Range(joker),1);
429 
430  z_surface.resize(1,1);
431 
432  // The spectra .
433  y.resize(f_grid.nelem());
434 
435  // The batch spectra.
436  ybatch.resize(no_profiles);
437 
438  // Loop over the number of profiles.
439  for (Index i = 0; i < no_profiles; ++ i)
440  {
441  ostringstream lat_os, lon_os;
442 
443  Index lat_prec = 3;
444  if(lat[i] < 0) lat_prec--;
445  if(abs(lat[i])>=10 )
446  {
447  lat_prec--;
448  if(abs(lat[i])>=100 ) lat_prec--;
449  }
450 
451  lat_os.setf (ios::showpoint | ios::fixed);
452  lat_os << setprecision((int)lat_prec) << lat[i];
453 
454  Index lon_prec = 4;
455  if(lon[i] < 0) lon_prec--;
456  if(abs(lon[i])>=10 )
457  {
458  lon_prec--;
459  if(abs(lon[i])>=100 ) lon_prec--;
460  }
461  lon_os.setf (ios::showpoint | ios::fixed);
462  lon_os << setprecision((int)lon_prec) << lon[i];
463 
464  sensor_los(0,0) =
465  180.0 - (asin(refellipsoid[0] * sin(sat_za_from_data[i] * DEG2RAD) /sensor_pos(0,0)))* RAD2DEG;
466 
467  //Reads the t_field_raw from file
468  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".t.xml",
469  t_field_raw, verbosity);
470 
471  //Reads the z_field_raw from file
472  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".z.xml",
473  z_field_raw, verbosity);
474 
475  //Reads the humidity from file - it is only an ArrayofTensor3
476  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".H2O.xml",
477  vmr_field_raw[0], verbosity);
478 
479  //Reads the pnd_field_raw for one particle
480  //xml_read_from_file("/rinax/storage/users/rekha/uk_data/profiles/new_obs/newest_forecastfields/reff100/profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".pnd100.xml", pnd_field_raw[0]);
481 
482  //xml_read_from_file(met_profile_pnd_path +"reff100_newformat/profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".pnd100.xml", pnd_field_raw[0]);
483 
484  xml_read_from_file(met_profile_pnd_path +"lwc_reff15/profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".pnd15.xml",
485  pnd_field_raw[0], verbosity);
486  //Write the profile number into a file.
487  // xml_write_to_file("profile_number.xml", i);
488 
489  // Set z_surface from lowest level of z_field
490  z_surface(0,0) = z_field_raw.data(0,0,0);
491 
492  /* The vmr_field_raw is an ArrayofArrayofTensor3 where the outer
493  array is for species.
494 
495  The oxygen and nitrogen VMRs are set to constant values of 0.209
496  and 0.782, respectively and are used along with humidity field
497  to generate *vmr_field_raw*.*/
498 
499  /*The second element of the species. The first 3 Tensors in the
500  array are the same . They are pressure grid, latitude grid and
501  longitude grid. The third tensor which is the vmr is set to a
502  constant value of 0.782, corresponding to N2.*/
503 
504  vmr_field_raw[1].resize(vmr_field_raw[0]);
505  vmr_field_raw[1].copy_grids(vmr_field_raw[0]);
506  vmr_field_raw[1] = 0.782;//vmr of N2
507 
508  /*the third element of the species. the first 3 Tensors in the
509  array are the same . They are pressure grid, latitude grid and
510  longitude grid. The third tensor which is the vmr is set to a
511  constant value of 0.209, corresponding to O2.*/
512  vmr_field_raw[2].resize(vmr_field_raw[0]);
513  vmr_field_raw[2].copy_grids(vmr_field_raw[0]);
514  vmr_field_raw[2] = 0.209;//vmr of O2
515 
516  const ConstVectorView tfr_p_grid = t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
517  // N_p is the number of elements in the pressure grid
518  Index N_p = tfr_p_grid.nelem();
519 
520  //Making a p_grid with the first and last element taken from the profile.
521  VectorNLogSpace(p_grid,
522  nelem_p_grid,
523  tfr_p_grid[0],
524  tfr_p_grid[N_p -1],
525  verbosity);
526 
527  /*To set the cloudbox limits, the lower and upper cloudbox limits
528  are to be set. The lower cloudbox limit is set to the lowest
529  pressure level. The upper level is the highest level where the
530  ice water content is non-zero.*/
531  Numeric cl_grid_min, cl_grid_max;
532 
533  //Lower limit = lowest pressure level of the original grid.
534  //Could it be the interpolated p_grid? FIXME STR
535  cl_grid_min = tfr_p_grid[0];
536 
537  // A counter for non-zero ice content
538  Index level_counter = 0;
539 
540  // Loop over all pressure levels
541  for (Index ip = 0; ip< N_p; ++ip)
542  {
543  //Checking for non-zero ice content. 0.001 is a threshold for
544  //ice water content.
545  // if((pnd_field_raw[0].data()(ip, 0, 0) > 0.001) || (pnd_field_raw[1](ip, 0, 0) > 0.001))
546  if(pnd_field_raw[0].data(ip, 0, 0) > 0.001)
547  {
548  ++level_counter;
549  //if non-zero ice content is found, it is set to upper
550  //cloudbox limit. Moreover, we take one level higher
551  // than the upper limit because we want the upper limit
552  //to have 0 pnd.
553  cl_grid_max = tfr_p_grid[ip +1];
554  }
555  }
556 
557  //cloudbox limits have dimensions 2*atmosphere_dim
558  cloudbox_limits.resize( atmosphere_dim*2 );
559 
560  //if there is no cloud in the considered profile, still we
561  //need to set the upper limit. I here set the first level
562  //for the upper cloudbox limit.
563  if(level_counter == 0)
564  {
565  cl_grid_max = p_grid[1];
566  }
567 
568  //Cloudbox is set.
569  cloudboxSetManually(cloudbox_on,
570  cloudbox_limits,
571  atmosphere_dim,
572  p_grid,
573  lat_grid,
574  lon_grid,
575  cl_grid_min,
576  cl_grid_max,
577  0,0,0,0,
578  verbosity);
579 
580  /*executing the met_profile_calc_agenda
581  Agenda communication variables are
582  Output of met_profile_calc_agenda : y
583  Input to met_profile_calc_agenda : t_field_raw,
584  z_field_raw, vmr_field_raw, pnd_field_raw, p_grid,
585  sensor_los, cloudbox_on, cloudbox_limits, z_surface, */
586 
587  met_profile_calc_agendaExecute (ws, y, t_field_raw, vmr_field_raw,
588  z_field_raw, pnd_field_raw, p_grid,
589  sensor_los, cloudbox_on,
590  cloudbox_limits, z_surface,
591  met_profile_calc_agenda );
592 
593  //putting in the spectra *y* for each profile, thus assigning y
594  //to the ith row of ybatch
595  ybatch[i] = y;
596 
597  }// closing the loop over profile basenames
598 }
599 
600 
601 /* Workspace method: Doxygen documentation will be auto-generated */
603  Workspace& ws,
604  //Output
605  ArrayOfVector& ybatch,
606  //Input
607  const ArrayOfArrayOfSpeciesTag& abs_species,
608  const Agenda& met_profile_calc_agenda,
609  const Vector& f_grid,
610  const Matrix& met_amsu_data,
611  const Matrix& sensor_pos,
612  const Vector& refellipsoid,
613  //Keyword
614  const Index& nelem_p_grid,
615  const String& met_profile_path,
616  const Verbosity& verbosity)
617 {
618  GriddedField3 t_field_raw;
619  GriddedField3 z_field_raw;
620  ArrayOfGriddedField3 vmr_field_raw;
621  ArrayOfGriddedField3 pnd_field_raw;
622  Vector p_grid;
623  Matrix sensor_los;
624  Index cloudbox_on = 0;
625  ArrayOfIndex cloudbox_limits;
626  Matrix z_surface;
627  Vector y;
628  Index no_profiles = met_amsu_data.nrows();
629  //Index no_profiles = met_profile_basenames.nelem();
630  // The humidity data is stored as an ArrayOfTensor3 whereas
631  // vmr_field_raw is an ArrayOfArrayOfTensor3
632  GriddedField3 vmr_field_raw_h2o;
633 
634  vmr_field_raw.resize(abs_species.nelem());
635 
636  y.resize(f_grid.nelem());
637  ybatch.resize(no_profiles);
638 
639  Vector sat_za_from_profile;
640  sat_za_from_profile = met_amsu_data(Range(joker),3);
641  Numeric sat_za;
642 
643  sensor_los.resize(1,1);
644 
645  Vector lat, lon;
646  lat = met_amsu_data(Range(joker),0);
647  lon = met_amsu_data(Range(joker),1);
648 
649  Vector oro_height;
650  oro_height = met_amsu_data(Range(joker),5);
651 
652  z_surface.resize(1,1);
653  for (Index i = 0; i < no_profiles; ++ i)
654  {
655  ostringstream lat_os, lon_os;
656 
657  Index lat_prec = 3;
658  if(lat[i] < 0) lat_prec--;
659  if(abs(lat[i])>=10 )
660  {
661  lat_prec--;
662  if(abs(lat[i])>=100 ) lat_prec--;
663  }
664 
665  lat_os.setf (ios::showpoint | ios::fixed);
666  lat_os << setprecision((int)lat_prec) << lat[i];
667 
668  Index lon_prec = 4;
669  if(lon[i] < 0) lon_prec--;
670  if(abs(lon[i])>=10 )
671  {
672  lon_prec--;
673  if(abs(lon[i])>=100 ) lon_prec--;
674  }
675  lon_os.setf (ios::showpoint | ios::fixed);
676  lon_os << setprecision((int)lon_prec) << lon[i];
677  cout<<lat_os.str()<<endl;
678  cout<<lon_os.str()<<endl;
679 
680 
681  sat_za = sat_za_from_profile[i];
682 
683  sensor_los(Range(joker),0) =
684  180.0 - (asin(refellipsoid[0] * sin(sat_za * PI/180.) /sensor_pos(0,0)))*180./PI;
685  cout<<"sensor_los"<<sat_za_from_profile[i]<<endl;
686  cout<<"sensor_los"<<sat_za<<endl;
687  cout<<"sensor_los"<<sensor_los<<endl;
688  //Reads the t_field_raw from file
689 
690  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".t.xml",
691  t_field_raw, verbosity);
692  //Reads the z_field_raw from file
693  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".z.xml",
694  z_field_raw, verbosity);
695 
696  //Reads the humidity from file - it is only an ArrayofTensor3
697  // The vmr_field_raw is an ArrayofArrayofTensor3 where the outer
698  // array is for species
699  xml_read_from_file(met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str() + ".H2O.xml",
700  vmr_field_raw_h2o, verbosity);
701  //xml_read_from_file("/home/home01/rekha/uk/profiles/sat_vmr/profile.lat_"+lat_os.str()//+".lon_"+lon_os.str() + ".H2O_es.xml",
702  // vmr_field_raw_h2o, verbosity);
703 
704  cout << "--------------------------------------------------------------------------"<<endl;
705  cout << "The file" << met_profile_path +"profile.lat_"+lat_os.str()+".lon_"+lon_os.str()<< "is executed now"<<endl;
706  cout << "--------------------------------------------------------------------------"<<endl;
707  xml_write_to_file("profile_number.xml", i, FILE_TYPE_ASCII, 0, verbosity);
708  // the first element of the species is water vapour.
709 
710  // N_p is the number of elements in the pressure grid
711  //z_surface(0,0) = oro_height[i]+ 0.01;
712  z_surface(0,0) = z_field_raw.data(0,0,0);
713  cout<<"z_surface"<<z_surface<<endl;
714  const ConstVectorView tfr_p_grid = t_field_raw.get_numeric_grid(GFIELD3_P_GRID);
715  Index N_p = tfr_p_grid.nelem();
716 
717  vmr_field_raw[0] = vmr_field_raw_h2o;
718 
719  // the second element of the species. the first 3 Tensors in the
720  //array are the same . They are pressure grid, latitude grid and
721  // longitude grid. The third tensor which is the vmr is set to a
722  // constant value of 0.782.
723  vmr_field_raw[1].resize(vmr_field_raw[0]);
724  vmr_field_raw[1].copy_grids(vmr_field_raw[0]);
725  vmr_field_raw[1].data(joker, joker, joker) = 0.782;
726 
727  // the second element of the species. the first 3 Tensors in the
728  //array are the same . They are pressure grid, latitude grid and
729  // longitude grid. The third tensor which is the vmr is set to a
730  // constant value of 0.209.
731  vmr_field_raw[2].resize(vmr_field_raw[0]);
732  vmr_field_raw[2].copy_grids(vmr_field_raw[0]);
733  vmr_field_raw[2].data(joker, joker, joker) = 0.209;
734 
735  //xml_write_to_file(met_profile_basenames[i]+ ".N2.xml", vmr_field_raw[1]);
736  //xml_write_to_file(met_profile_basenames[i]+ ".O2.xml", vmr_field_raw[2]);
737 
738  //Making a p_grid with the first and last element taken from the profile.
739  // this is because of the extrapolation problem.
740 
741  VectorNLogSpace(p_grid,
742  nelem_p_grid,
743  tfr_p_grid[0],
744  tfr_p_grid[N_p -1],
745  verbosity);
746  cout<<"t_field_raw[0](0,0,0)"<<tfr_p_grid[0]<<endl;
747  cout<<"t_field_raw[0](N_p -1,0,0)"<<tfr_p_grid[N_p -1] <<endl;
748  xml_write_to_file("p_grid.xml", p_grid, FILE_TYPE_ASCII, 0, verbosity);
749 
750  // executing the met_profile_calc_agenda
751  met_profile_calc_agendaExecute (ws, y, t_field_raw, vmr_field_raw,
752  z_field_raw, pnd_field_raw, p_grid,
753  sensor_los, cloudbox_on,
754  cloudbox_limits, z_surface,
755  met_profile_calc_agenda );
756 
757  //putting in the spectra *y* for each profile
758  ybatch[i] = y;
759 
760  }// closing the loop over profile basenames
761 }
762 
763 
764 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void cloudboxSetManually(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Numeric &p1, const Numeric &p2, const Numeric &lat1, const Numeric &lat2, const Numeric &lon1, const Numeric &lon2, const Verbosity &)
WORKSPACE METHOD: cloudboxSetManually.
Definition: m_cloudbox.cc:365
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
ConstVectorView get_numeric_grid(Index i) const
Get a numeric grid.
The Vector class.
Definition: matpackI.h:556
int arts_omp_get_thread_num()
Wrapper for omp_get_thread_num.
Definition: arts_omp.cc:83
The range class.
Definition: matpackI.h:148
void ybatchMetProfiles(Workspace &ws, ArrayOfVector &ybatch, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &met_profile_calc_agenda, const Vector &f_grid, const Matrix &met_amsu_data, const Matrix &sensor_pos, const Vector &refellipsoid, const Vector &lat_grid, const Vector &lon_grid, const Index &atmosphere_dim, const ArrayOfSingleScatteringData &scat_data_array, const Index &nelem_p_grid, const String &met_profile_path, const String &met_profile_pnd_path, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMetProfiles.
Definition: m_batch.cc:373
const Numeric RAD2DEG
const Index GFIELD3_P_GRID
This file contains basic functions to handle XML data files.
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:146
const Numeric DEG2RAD
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
const Numeric PI
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
void met_profile_calc_agendaExecute(Workspace &ws, Vector &y, const GriddedField3 &t_field_raw, const ArrayOfGriddedField3 &vmr_field_raw, const GriddedField3 &z_field_raw, const ArrayOfGriddedField3 &pnd_field_raw, const Vector &p_grid, const Matrix &sensor_los, const Index cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Matrix &z_surface, const Agenda &input_agenda)
Definition: auto_md.cc:14367
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
#define CREATE_OUTS
Definition: messages.h:216
#define CREATE_OUT1
Definition: messages.h:212
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:831
#define abs(x)
Definition: continua.cc:20458
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
void ybatchMetProfilesClear(Workspace &ws, ArrayOfVector &ybatch, const ArrayOfArrayOfSpeciesTag &abs_species, const Agenda &met_profile_calc_agenda, const Vector &f_grid, const Matrix &met_amsu_data, const Matrix &sensor_pos, const Vector &refellipsoid, const Index &nelem_p_grid, const String &met_profile_path, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchMetProfilesClear.
Definition: m_batch.cc:602
void VectorNLogSpace(Vector &x, const Index &n, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
void ForLoop(Workspace &ws, const Agenda &forloop_agenda, const Index &start, const Index &stop, const Index &step, const Verbosity &verbosity)
WORKSPACE METHOD: ForLoop.
Definition: m_batch.cc:68
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
A constant view of a Vector.
Definition: matpackI.h:292
void MatrixExtractFromTensor3(Matrix &m, const Tensor3 &t3, const Index &index, const String &direction, const Verbosity &)
WORKSPACE METHOD: MatrixExtractFromTensor3.
Definition: m_batch.cc:144
Workspace class.
Definition: workspace_ng.h:47
void ybatch_calc_agendaExecute(Workspace &ws, Vector &y, ArrayOfVector &y_aux, Matrix &jacobian, const Index ybatch_index, const Agenda &input_agenda)
Definition: auto_md.cc:15174
void VectorExtractFromMatrix(Vector &v, const Matrix &m, const Index &index, const String &direction, const Verbosity &)
WORKSPACE METHOD: VectorExtractFromMatrix.
Definition: m_batch.cc:94
Header file for helper functions for OpenMP.
void forloop_agendaExecute(Workspace &ws, const Index forloop_index, const Agenda &input_agenda)
Definition: auto_md.cc:13569
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
Definition: xml_io.cc:982
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
void ybatchCalc(Workspace &ws, ArrayOfVector &ybatch, ArrayOfArrayOfVector &ybatch_aux, ArrayOfMatrix &ybatch_jacobians, const Index &ybatch_start, const Index &ybatch_n, const Agenda &ybatch_calc_agenda, const Index &robust, const Verbosity &verbosity)
WORKSPACE METHOD: ybatchCalc.
Definition: m_batch.cc:208
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580