m_abs.cc

Go to the documentation of this file.

/* Copyright (C) 2000, 2001 Stefan Buehler   <sbuehler@uni-bremen.de>
                            Patrick Eriksson <patrick@rss.chalmers.se>
                            Axel von Engeln  <engeln@uni-bremen.de>
                            Thomas Kuhn      <tkuhn@uni-bremen.de>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
// 
 
#include <math.h>
#include <algorithm>
#include "arts.h"
#include "matpackI.h"
#include "array.h"
#include "messages.h"
#include "file.h"
#include "absorption.h"
#include "auto_wsv.h"
#include "auto_md.h"
#include "math_funcs.h"
#include "make_array.h"
#include "atm_funcs.h"
#include "continua.h"
#include "make_vector.h"
 
void lines_per_tgSetEmpty(// WS Output:
                          ArrayOfArrayOfLineRecord& lines_per_tg,
                          // WS Input:
                          const TagGroups& tgs)
{
  // Make lines_per_tg the right size:
  lines_per_tg.resize( tgs.nelem() );
  
  for (Index i=0; i<tgs.nelem(); ++i)
    {
      lines_per_tg[i].resize(0);
    }
}
 
void linesReadFromHitran(// WS Output:
                         ArrayOfLineRecord& lines,
                          // Control Parameters:
                         const String& filename,
                         const Numeric& fmin,
                         const Numeric& fmax)
{
  ifstream is;
 
  // Reset lines in case it already existed:
  lines.resize(0);
 
  out2 << "  Reading file: " << filename << '\n';
  open_input_file(is, filename);
 
  bool go_on = true;
  while ( go_on )
    {
      LineRecord lr;
      if ( lr.ReadFromHitranStream(is) )
        {
          // If we are here the read function has reached eof and has
          // returned no data.
          go_on = false;
        }
      else
        {
          if ( fmin <= lr.F() )
            {
              if ( lr.F() <= fmax )
                lines.push_back(lr);
              else
                go_on = false;
            }
        }
    }
  out2 << "  Read " << lines.nelem() << " lines.\n";
}
 
 
void linesReadFromHitran2004(// WS Output:
                         ArrayOfLineRecord& lines,
                          // Control Parameters:
                         const String& filename,
                         const Numeric& fmin,
                         const Numeric& fmax)
{
  ifstream is;
 
  // Reset lines in case it already existed:
  lines.resize(0);
 
  out2 << "  Reading file: " << filename << '\n';
  open_input_file(is, filename);
 
  bool go_on = true;
  while ( go_on )
    {
      LineRecord lr;
      if ( lr.ReadFromHitran2004Stream(is) )
        {
          // If we are here the read function has reached eof and has
          // returned no data.
          go_on = false;
        }
      else
        {
          if ( fmin <= lr.F() )
            {
              if ( lr.F() <= fmax )
                lines.push_back(lr);
              else
                go_on = false;
            }
        }
    }
  out2 << "  Read " << lines.nelem() << " lines.\n";
}
 
 
void linesReadFromMytran2(// WS Output:
                          ArrayOfLineRecord& lines,
                          // Control Parameters:
                          const String& filename,
                          const Numeric& fmin,
                          const Numeric& fmax)
{
  ifstream is;
 
  // Reset lines in case it already existed:
  lines.resize(0);
 
  out2 << "  Reading file: " << filename << '\n';
  open_input_file(is, filename);
 
  bool go_on = true;
  while ( go_on )
    {
      LineRecord lr;
      if ( lr.ReadFromMytran2Stream(is) )
        {
          // If we are here the read function has reached eof and has
          // returned no data.
          go_on = false;
        }
      else
        {
          // lines are not necessarily frequency sorted 
          if ( fmin <= lr.F() )
            if ( lr.F() <= fmax )
              lines.push_back(lr);
        }
    }
  out2 << "  Read " << lines.nelem() << " lines.\n";
}
 
 
void linesReadFromJpl(// WS Output:
                      ArrayOfLineRecord& lines,
                      // Control Parameters:
                      const String& filename,
                      const Numeric& fmin,
                      const Numeric& fmax)
{
  ifstream is;
 
  // Reset lines in case it already existed:
  lines.resize(0);
 
  out2 << "  Reading file: " << filename << '\n';
  open_input_file(is, filename);
 
  bool go_on = true;
  while ( go_on )
    {
      LineRecord lr;
      if ( lr.ReadFromJplStream(is) )
        {
          // If we are here the read function has reached eof and has
          // returned no data.
          go_on = false;
        }
      else
        {
          // we expect lines to be sorted
          if ( fmin <= lr.F() )
            {
              if ( lr.F() <= fmax )
                lines.push_back(lr);
              else
                go_on = false;
            }
        }
    }
  out2 << "  Read " << lines.nelem() << " lines.\n";
}
 
 
void linesReadFromArts(// WS Output:
                       ArrayOfLineRecord& lines,
                       // Control Parameters:
                       const String& filename,
                       const Numeric& fmin,
                       const Numeric& fmax)
{
  // The input stream:
  ifstream is;
 
  // We will use this line record to read in the line records in the
  // file one after another:
  LineRecord lr;
 
  // Reset lines in case it already existed:
  lines.resize(0);
 
  out2 << "  Reading file: " << filename << '\n';
  open_input_file(is, filename);
 
  // Get version tag and check that it corresponds to the current version.
  {
    String v;
    is >> v;
    if ( v!=lr.Version() )
      {
        ostringstream os;
        
        // If what we read is the version String, it should have at elast 9 characters.
        if ( 9 <= v.nelem() )
          {
            if ( "ARTSCAT" == v.substr(0,7) )
            {
              os << "The ARTS line file you are trying contains a version tag "
                 << "different from the current version.\n"
                 << "Tag in file:     " << v << "\n"
                 << "Current version: " << lr.Version();
              throw runtime_error(os.str());
            }
           }
 
            os << "The ARTS line file you are trying to read does not contain a valid version tag.\n"
            << "Probably it was created with an older version of ARTS that used different units.";
            throw runtime_error(os.str());
      }
  }
 
  bool go_on = true;
  while ( go_on )
    {
      if ( lr.ReadFromArtsStream(is) )
        {
          // If we are here the read function has reached eof and has
          // returned no data.
          go_on = false;
        }
      else
        {
          // lines are not necessarily frequency sorted 
          if ( fmin <= lr.F() && lr.F() <= fmax )
            {
              lines.push_back(lr);
            }
        }
    }
  out2 << "  Read " << lines.nelem() << " lines.\n";
}
 
void linesElowToJoule(// WS Output:
                      ArrayOfLineRecord& lines )
{
  for ( Index i=0; i<lines.nelem(); ++i )
    lines[i].melow = wavenumber_to_joule(lines[i].melow); 
}
 
void lines_per_tgReadFromCatalogues(// WS Output:
                                    ArrayOfArrayOfLineRecord& lines_per_tg,
                                    // WS Input:
                                    const TagGroups& tgs,
                                    // Control Parameters:
                                    const ArrayOfString& filenames,
                                    const ArrayOfString& formats,
                                    const Vector& fmin,
                                    const Vector& fmax)
{
  const Index n_tg   = tgs.nelem();     // # tag groups
  const Index n_cat = filenames.nelem();        // # tag Catalogues
 
  // Check that dimensions of the keyword parameters are consistent
  // (must all be the same). 
 
  if ( n_cat != formats.nelem() ||
       n_cat != fmin.nelem() ||
       n_cat != fmax.nelem() )
    {
      ostringstream os;
      os << "lines_per_tgReadFromCatalogues: All keyword\n"
         << "parameters must get the same number of arguments.\n"
         << "You specified:\n"
         << "filenames: " << n_cat         << "\n"
         << "formats:   " << formats.nelem() << "\n"
         << "fmin:      " << fmin.nelem()    << "\n"
         << "fmax:      " << fmax.nelem();
      throw runtime_error(os.str());
    }
  
  // Furthermore, the dimension must be
  // smaller than or equal to the number of tag groups.
 
  if ( n_cat > n_tg )
    {
      ostringstream os;
      os << "lines_per_tgReadFromCatalogues: You specified more\n"
         << "catalugues than you have tag groups.\n"
         << "You specified:\n"
         << "Catalogues: " << n_cat << "\n"
         << "tgs: " << n_tg;
      throw runtime_error(os.str());
    }
 
  // There must be at least one tag group and at least one catalogue:
 
  if ( n_cat < 1 ||
       n_tg   < 1 )
    {
      ostringstream os;
      os << "lines_per_tgReadFromCatalogues: You must have at\n"
         << "least one catalogue and at least one tag group.\n"
         << "You specified:\n"
         << "Catalogues: " << n_cat << "\n"
         << "tgs: " << n_tg;
      throw runtime_error(os.str());
    }
 
  // There can be repetitions in the keyword paramters. We want to read
  // and process each catalogue only once, so we'll compile a set of
  // real catalogues, along with a data structure that tells us which
  // tag groups should use this catalogue.
 
  MakeArray<String>      real_filenames ( filenames[0] );
  MakeArray<String>      real_formats   ( formats[0]   );
  MakeArray<Numeric>     real_fmin      ( fmin[0]      );
  MakeArray<Numeric>     real_fmax      ( fmax[0]      );
 
  Array< ArrayOfIndex > real_tgs(1);
  real_tgs[0].resize(1);
  real_tgs[0][0] = 0;
 
  // The last specified catalogue, to which we should assign all
  // remaining lines. Index of this one in real_ arrays.
  Index last_cat = 0;           
 
  for ( Index i=1; i<n_tg; ++i )
    {
      // Is there a catalogue specified?
      if ( n_cat > i )
        {
          // Yes, there is a catalogue.
 
          // Has this been specified before?
          // We use the STL find algorithm to look for the catalogue
          // name in the real_catalogues. Find returns an iterator, so
          // to get an index we have to take the difference to
          // .begin(). 
          const Index that_cat =
            (Index)(find( real_filenames.begin(),
                          real_filenames.end(),
                          filenames[i] ) - real_filenames.begin());
          if ( that_cat < real_filenames.nelem() )
            {
              // Yes, it has been specified before
              // ==> Assign to that catalogue
              real_tgs[that_cat].push_back(i);
 
              // Verify, that format, fmin, and fmax are consistent:
              if ( formats[i] != real_formats[that_cat] ||
                   fmin[i]    != real_fmin[that_cat]    ||
                   fmax[i]    != real_fmax[that_cat]       )
                {
                  ostringstream os;
                  os << "lines_per_tgReadFromCatalogues: If you specify the\n"
                     << "same catalogue repeatedly, format, fmin, and fmax must be\n"
                     << "consistent. There is an inconsistency between\n"
                     << "catalogue " << that_cat << " and " << i << ".";
                  throw runtime_error(os.str());
                }
            }
          else
            {
              // No, it has not been specified before.
              // ==> Add an entry to real_tgs and the other real_ variables:
              real_tgs.push_back( MakeArray<Index>(i) );
 
              real_filenames.push_back( filenames[i] );
              real_formats.push_back  ( formats[i]   );  
              real_fmin.push_back     ( fmin[i]      );     
              real_fmax.push_back     ( fmax[i]      );
 
              last_cat = i;     // assign remainder of lines to this
                                // catalogue, if there is no other catalogue.
            }
        }
      else
        {
          // No, there is no catalogue.
          // ==> Assign to the last catalogue
          real_tgs[last_cat].push_back(i);
        }
    }
 
  Index n_real_cat = real_filenames.nelem(); // # real catalogues to read
 
  // Some output to low priority stream:
  out3 << "  Catalogues to read and tgs for which these will be used:\n";
  for ( Index i=0; i<n_real_cat; ++i )
    {
      out3 << "  " << real_filenames[i] << ":";
      for ( Index s=0; s<real_tgs[i].nelem(); ++s )
        out3 << " " << real_tgs[i][s];
      out3 << "\n";
    }
 
  // Make lines_per_tg the right size:
  lines_per_tg.resize( tgs.nelem() );
 
  // Loop through the catalogues to read:
  for ( Index i=0; i<n_real_cat; ++i )
    {
      ArrayOfLineRecord   lines;
 
      // Read catalogue:
 
      if ( "HITRAN96"==real_formats[i] )
        {
          linesReadFromHitran( lines, real_filenames[i], real_fmin[i], real_fmax[i] );
        }
      else if ( "HITRAN04"==real_formats[i] )
        {
          linesReadFromHitran2004( lines, real_filenames[i], real_fmin[i], real_fmax[i] );
        }
      else if ( "MYTRAN2"==real_formats[i] )
        {
          linesReadFromMytran2( lines, real_filenames[i], real_fmin[i], real_fmax[i] );
        }
      else if ( "JPL"==real_formats[i] )
        {
          linesReadFromJpl( lines, real_filenames[i], real_fmin[i], real_fmax[i] );
        }
      else if ( "ARTS"==real_formats[i] )
        {
          linesReadFromArts( lines, real_filenames[i], real_fmin[i], real_fmax[i] );
        }
      else
        {
          ostringstream os;
          os << "lines_per_tgReadFromCatalogues: You specified the\n"
             << "format `" << real_formats[i] << "', which is unknown.\n"
             << "Allowd formats are: HITRAN96, HITRAN04, MYTRAN2, JPL, ARTS.";
          throw runtime_error(os.str());
        }
 
      // We need to make subset tgs for the groups that should
      // be read from this catalogue.
      TagGroups  these_tgs(real_tgs[i].nelem());
      for ( Index s=0; s<real_tgs[i].nelem(); ++s )
        {
          these_tgs[s] = tgs[real_tgs[i][s]];
        }
 
      // Create these_lines_per_tg:
      ArrayOfArrayOfLineRecord these_lines_per_tg;
      lines_per_tgCreateFromLines( these_lines_per_tg, lines, these_tgs );
 
      // Put these lines in the right place in lines_per_tg:
      for ( Index s=0; s<real_tgs[i].nelem(); ++s )
        {
          lines_per_tg[real_tgs[i][s]] = these_lines_per_tg[s];
        }
    }
}
 
 
void lines_per_tgCreateFromLines(// WS Output:
                                  ArrayOfArrayOfLineRecord& lines_per_tg,
                                  // WS Input:
                                  const ArrayOfLineRecord&   lines,
                                  const TagGroups&           tgs)
{
  // The species lookup data:
  extern const Array<SpeciesRecord> species_data;
 
  // As a safety feature, we will watch out for the case that a
  // species is included in the calculation, but not all lines are
  // used. For this we need an array to flag the used species:
  
  // For some weird reason, Arrays of bool do not work, although all
  // other types seem to work fine. So in this single case, I'll use
  // the stl vector directly. The other place where this is done is in
  // the function executor in main.cc.
  // FIXME: Fix this when Array<bool> works.
  std::vector<bool> species_used (species_data.nelem(),false);
      
  // Make lines_per_tg the right size:
  lines_per_tg = ArrayOfArrayOfLineRecord(tgs.nelem());
 
  // Unfortunately, MTL conatains a bug that leads to all elements of
  // the outer Array of an Array<Array>> pointing to the same data
  // after creation. So we need to fix this explicitly:
  for ( Index i=0; i<lines_per_tg.nelem(); ++i )
    lines_per_tg[i] = ArrayOfLineRecord();
 
  // Loop all lines in the input line list:
  for ( Index i=0; i<lines.nelem(); ++i )
    {
      // Get a convenient reference to the current line:
      const LineRecord& this_line = lines[i];
 
      // We have to test each tag group in turn (in the order in which
      // they appear in the controlfile). The line is assigned to the
      // first tag group that fits.
 
      // The flag found is used to break the for loops when the right
      bool found = false;
 
      // We need to define j here, since we need the value outside the
      // for loop:
      Index j;
 
      // Loop the tag groups:
      for ( j=0; j<tgs.nelem() && !found ; ++j ) 
        {
          // A tag group can contain several tags:
          for ( Index k=0; k<tgs[j].nelem() && !found; ++k )
            {
              // Get a reference to the current tag (not really
              // necessary, but makes for nicer notation below):
              const OneTag& this_tag = tgs[j][k];
 
              // Now we will test different attributes of the line
              // against matching attributes of the tag. If any
              // attribute does not match, we continue with the next tag
              // in the tag group. (Exception: Species, see just below.)
 
              // Test species. If this attribute does not match we don`t
              // have to test the other tags in this group, since all
              // tags must belong to the same species.
              if ( this_tag.Species() != this_line.Species() ) break;
 
              // Test isotope. The isotope can either match directly, or
              // the Isotope of the tag can be one larger than the
              // number of isotopes, which means `all'. Test the second
              // condition first, since this will probably be more often
              // used.
              if ( this_tag.Isotope() != this_line.SpeciesData().Isotope().nelem() )
                if ( this_tag.Isotope() != this_line.Isotope() )
                  continue;
 
            // Test frequncy range we take both the lower (Lf) and the
            // upper (Uf) border frequency to include the `equal' case.
            // Both Lf and Uf can also be negative, which means `no limit'
 
            // Take the lower limit first:
              if ( this_tag.Lf() >= 0 )
                if ( this_tag.Lf() > this_line.F() )
                  continue;
 
            // Then the upper limit:
              if ( this_tag.Uf() >= 0 )
                if ( this_tag.Uf() < this_line.F() )
                  continue;
 
            // When we get here, this_tag has survived all tests. That
            // means it matches the line perfectly!
              found = true;
            }
        }
 
      // If a matching tag was found, this line can be used in the
      // calculation. Add it to the line list for this tag group.
      if (found)
        {
          // We have to use j-1 here, since j was still increased by
          // one after the matching tag has been found.
          lines_per_tg[j-1].push_back(this_line);
 
          // Flag this species as used, if not already done:
          if ( !species_used[this_line.Species()] )
            species_used[this_line.Species()] = true;
        }
      else
        {
          // Safety feature: Issue a warning messages if the lines for a
          // species are only partly covered by tags.
          if ( species_used[this_line.Species()] )
            {
              out2 << "Your tags include other lines of species "
                   << this_line.SpeciesData().Name()
                   << ",\n"
                   << "why do you not include line "
                   << i
                   << " (at "
                   << this_line.F()
                   << " Hz)?\n";
            }
        }
   }
 
  // Write some information to the lowest priority output stream.
  for (Index i=0; i<tgs.nelem(); ++i)
    {
        out3 << "  " << i << ":";
 
        for (Index s=0; s<tgs[i].nelem(); ++s)
          out3 << " " << tgs[i][s].Name();
 
        out3 << ": " << lines_per_tg[i].nelem() << " lines\n";
    }
 
}
 
void lines_per_tgAddMirrorLines(// WS Output:
                                ArrayOfArrayOfLineRecord& lines_per_tg)
{
  // We will simply append the mirror lines after the original
  // lines. This way we don't have to make a backup copy of the
  // original lines. 
 
  for ( Index i=0; i<lines_per_tg.nelem(); ++i )
    {
      // Get a reference to the current list of lines to save typing:
      ArrayOfLineRecord& ll = lines_per_tg[i];
      
      // For increased efficiency, reserve the necessary space:
      ll.reserve(2*ll.nelem());
 
      // Loop through all lines of this tag group:
      {
        // It is important that we determine the size of ll *before*
        // we start the loop. After all, we are adding elements. And
        // we cerainly don't want to continue looping the newly added
        // elements, we want to loop only the original elements.
        Index n=ll.nelem();
        for ( Index j=0; j<n; ++j )
          {
            LineRecord dummy = ll[j];
            dummy.setF( -dummy.F() );
            //      cout << "Adding ML at f = " << dummy.F() << "\n";
            ll.push_back(dummy);
          }
      }
    }
 
}
 
void lines_per_tgCompact(// WS Output:
                         ArrayOfArrayOfLineRecord& lines_per_tg,
                         // WS Input:
                         const ArrayOfLineshapeSpec& lineshape,
                         const Vector& f_mono)
{
 
  // Make sure lines_per_tg and lineshape have the same dimension:
  if ( lines_per_tg.nelem() != lineshape.nelem() ) 
    {
      ostringstream os;
      os << "Dimension of lines_per_tg does\n"
         << "not match that of lineshape.";
      throw runtime_error(os.str());
    }
  
  // Make sure that the frequency grid is properly sorted:
  for ( Index s=0; s<f_mono.nelem()-1; ++s )
    {
      if ( f_mono[s+1] <= f_mono[s] )
        {
          ostringstream os;
          os << "The frequency grid f_mono is not properly sorted.\n"
             << "f_mono[" << s << "] = " << f_mono[s] << "\n"
             << "f_mono[" << s+1 << "] = " << f_mono[s+1];
          throw runtime_error(os.str());
        }
    }
 
  // Cycle through all tag groups:
  for ( Index i=0; i<lines_per_tg.nelem(); ++i )
    {
      // Get cutoff frequency of this tag group:
      Numeric cutoff = lineshape[i].Cutoff();
 
      // Check whether cutoff is defined:
      if ( cutoff != -1)
        {
          // Get a reference to the current list of lines to save typing:
          ArrayOfLineRecord& ll = lines_per_tg[i];
 
          // Calculate the borders:
          Numeric upp = f_mono[f_mono.nelem()-1] + cutoff;
          Numeric low = f_mono[0] - cutoff;
 
          // Cycle through all lines within this tag group. 
          Array<ArrayOfLineRecord::iterator> keep;
          for ( ArrayOfLineRecord::iterator j=ll.begin(); j<ll.end(); ++j )
            {
              // Center frequency:
              const Numeric F0 = j->F();
 
              if ( ( F0 >= low) && ( F0 <= upp) )
                {
                  // Build list of elements which should be kept
                  keep.push_back (j);
                  // The original implementation just erased the non-wanted
                  // elements from the array. Problem: On every erase the
                  // whole array will be copied which actually kills
                  // performance.
                }
            }
 
          // If next comparison is false, some elements have to be removed
          if (keep.nelem () != ll.nelem ())
            {
              ArrayOfLineRecord nll;
              // Copy all elements that should be kept to a new Array
              for (Array<ArrayOfLineRecord::iterator>::iterator j
                   = keep.begin(); j != keep.end(); j++)
                {
                  nll.push_back (**j);
                }
              // Overwrite the old array with the new one
              ll.resize (nll.nelem ());
              ll = nll;
            }
        }
    }
}
 
 
void linesWriteAscii(// WS Input:
                      const ArrayOfLineRecord& lines,
                      // Control Parameters:
                      const String& f)
{
  String filename = f;
  
  // If filename is empty then set the default filename:
  if ( "" == filename )
    {
      extern const String out_basename;                       
      filename = out_basename+".lines.aa";
    }
 
  ofstream os;
 
  out2 << "  Writing file: " << filename << '\n';
  open_output_file(os, filename);
 
  write_lines_to_stream(os,lines);
}
 
 
void lines_per_tgWriteAscii(// WS Input:
                              const ArrayOfArrayOfLineRecord& lines_per_tg,
                              // Control Parameters:
                              const String& f)
{
  String filename = f;
  
  // If filename is empty then set the default filename:
  if ( "" == filename )
    {
      extern const String out_basename;                       
      filename = out_basename+".lines_per_tg.aa";
    }
 
  ofstream os;
 
  out2 << "  Writing file: " << filename << '\n';
  open_output_file(os, filename);
 
  os << lines_per_tg.nelem() << '\n';
 
  for ( Index i=0; i<lines_per_tg.nelem(); ++i )
    {
      const ArrayOfLineRecord& lines = lines_per_tg[i];
      write_lines_to_stream( os, lines );
    }
}
 
 
void tgsDefine(// WS Output:
               TagGroups& tgs,
               // Control Parameters:
               const ArrayOfString& tags)
{
  tgs.resize(tags.nelem());
 
  //cout << "Tags: " << tags << "\n";
 
  // Each element of the array of Strings tags defines one tag
  // group. Let's work through them one by one.
  for ( Index i=0; i<tags.nelem(); ++i )
    {
      // There can be a comma separated list of tag definitions, so we
      // need to break the String apart at the commas.
      ArrayOfString tag_def;
 
      bool go_on = true;
      String these_tags = tags[i];
      while (go_on)
        {
          Index n = (Index)these_tags.find(',');
          if ( n == these_tags.npos ) // npos indicates `not found'
            {
              // There are no more commas.
              //              cout << "these_tags: (" << these_tags << ")\n";
              tag_def.push_back(these_tags);
              go_on = false;
            }
          else
            {
              tag_def.push_back(these_tags.substr(0,n));
              these_tags.erase(0,n+1);
            }
        }
 
      // tag_def now holds the different tag Strings for this group.
      //      cout << "tag_def =\n" << tag_def << endl;
 
      // Set size to zero, in case the method has been called before.
      tgs[i].resize(0);
 
      for ( Index s=0; s<tag_def.nelem(); ++s )
        {
          // Remove leading whitespace, if there is any:
          while ( ' '  == tag_def[s][0] ||
                  '\t' == tag_def[s][0]    )    tag_def[s].erase(0,1);
 
          OneTag this_tag(tag_def[s]);
 
          // Safety check: For s>0 check that the tags belong to the same species.
          if (s>0)
            if ( tgs[i][0].Species() != this_tag.Species() )
              throw runtime_error("Tags in a tag group must belong to the same species.");
 
          tgs[i].push_back(this_tag);
        }
    }
 
  // Print list of tag groups to the most verbose output stream:
  out3 << "  Defined tag groups:";
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      out3 << "\n  " << i << ":";
      for ( Index s=0; s<tgs[i].nelem(); ++s )
        {
          out3 << " " << tgs[i][s].Name();
        }
    }
  out3 << '\n';
 
//  cout << endl << endl << tgs << endl;
}
 
void tgsDefineAllInScenario(// WS Output:
                            TagGroups& tgs,
                            // Control Parameters:
                            const String& basename)
{
  // Species lookup data:
  extern const Array<SpeciesRecord> species_data;
 
  // We want to make lists of included and excluded species:
  ArrayOfString included(0), excluded(0);
 
  tgs.resize(0);
 
  for ( Index i=0; i<species_data.nelem(); ++i )
    {
      const String specname = species_data[i].Name();
      const String filename = basename + "." + specname + ".aa";
 
      // Try to open VMR file:
      try
        {
          ifstream file;
          open_input_file(file, filename);
 
          // Ok, if we get here the file was found.
 
          // Add to included list:
          included.push_back(specname);
 
          // Convert name of species to a OneTag object:
          OneTag this_tag(specname);
 
          // Create Array of OneTags with length 1 (our tag group has only one tag):
          Array<OneTag> this_group(1);
          this_group[0] = this_tag;
 
          // Add this tag group to tgs:
          tgs.push_back(this_group);
        }
      catch (runtime_error x)
        {
          // Ok, the file for the species could not be found.
          excluded.push_back(specname);
        }
    }
  
  // Some nice output:
  out2 << "  Included Species (" << included.nelem() << "):\n";
  for ( Index i=0; i<included.nelem(); ++i )
    out2 << "     " << included[i] << "\n";
 
  out2 << "  Excluded Species (" << excluded.nelem() << "):\n";
  for ( Index i=0; i<excluded.nelem(); ++i )
    out2 << "     " << excluded[i] << "\n";
}
 
void lineshapeDefine(// WS Output:
                     ArrayOfLineshapeSpec&    lineshape,
                     // WS Input:
                     const TagGroups&         tgs,
                     const String&            shape,
                     const String&            normalizationfactor,
                     const Numeric&           cutoff)
{
  // Make lineshape and normalization factor data visible:
  extern const Array<LineshapeRecord> lineshape_data;
  extern const Array<LineshapeNormRecord> lineshape_norm_data;
 
 
  // generate the right number of elements
  Index tag_sz = tgs.nelem();
  lineshape.resize(tag_sz);
 
  // Is this lineshape available?
  Index found0=-1;
  for ( Index i=0; i<lineshape_data.nelem() && (found0 == -1) ; ++i )
    {
      const String& str = lineshape_data[i].Name();
      if (str == shape) 
        {
          out2 << "  Selected lineshape: " << str << "\n";
          found0=i;
        }
    }
 
  // Is this normalization to the lineshape available?
  Index found1=-1;
  for ( Index i=0; i<lineshape_norm_data.nelem() && (found1 == -1); ++i )
    {
      const String& str = lineshape_norm_data[i].Name();
      if (str == normalizationfactor) 
        {
          out2 << "  Selected normalization factor  : " << normalizationfactor << "\n";
 
          if ( (cutoff != -1) && (cutoff < 0.0) )
            throw runtime_error("  Cutoff must be -1 or positive.");
          out2 << "  Selected cutoff frequency [Hz] : " << cutoff << "\n";
          found1=i;
        }
    }
 
 
  // did we find the lineshape and normalization factor?
  if (found0 == -1)
    throw runtime_error("Selected lineshape not available.");
  if (found1 == -1)
    throw runtime_error("Selected normalization to lineshape not available.");
 
 
  // now set the lineshape  
  for (Index i=0; i<tag_sz; i++)
    {
      lineshape[i].SetInd_ls( found0 );
      lineshape[i].SetInd_lsn( found1 );
      lineshape[i].SetCutoff( cutoff );
    }
}
 
void lineshape_per_tgDefine(// WS Output:
                            ArrayOfLineshapeSpec& lineshape,
                            // WS Input:
                            const TagGroups&      tgs,
                            const ArrayOfString&  shape,
                            const ArrayOfString&  normalizationfactor,
                            const Vector&         cutoff )
{
  // Make lineshape and normalization factor data visible:
  extern const Array<LineshapeRecord> lineshape_data;
  extern const Array<LineshapeNormRecord> lineshape_norm_data;
 
  // check that the number of elements are equal
  Index tg_sz = tgs.nelem();
  if ( (tg_sz != shape.nelem()) ||
       (tg_sz != normalizationfactor.nelem()) || 
       (tg_sz != cutoff.nelem()) )
    {
      ostringstream os;
      os << "lineshape_per_tgDefine: number of elements does\n"
         << "not match the number of tag groups defined.";
      throw runtime_error(os.str());
    }
      
 
  // generate the right number of elements
  lineshape.resize(tg_sz);
 
  // Is this lineshape available?
  for (Index k=0; k<tg_sz; ++k)
    {
      Index found0=-1;
      for ( Index i=0; i<lineshape_data.nelem() && (found0 == -1); ++i )
        {
          const String& str = lineshape_data[i].Name();
          if (str == shape[k]) 
            {
              out2 << "  Tag Group: [";
              for (Index s=0; s<tgs[k].nelem()-1; ++s)
                out2 << tgs[k][s].Name() << ", "; 
              out2 << tgs[k][tgs[k].nelem()-1].Name() << "]\n";
              out2 << "  Selected lineshape: " << str << "\n";
              found0=i;
            }
        }
 
      // Is this normalization to the lineshape available?
      Index found1=-1;
      for ( Index i=0; i<lineshape_norm_data.nelem() && (found1 == -1); ++i )
        {
          const String& str = lineshape_norm_data[i].Name();
          if (str == normalizationfactor[k]) 
            {
              out2 << "  Selected normalization factor: " << normalizationfactor[k] << "\n";
              if ( (cutoff[k] != -1) && (cutoff[k] < 0.0) )
                throw runtime_error("  Cutoff must be -1 or positive.");
              out2 << "  Selected cutoff frequency    : " << cutoff[k] << "\n";
              found1=i;
            }
        }
 
 
      // did we find the lineshape and normalization factor?
      if (found0 == -1)
        {
          ostringstream os;
          os << "Selected lineshape not available: "<< shape[k] <<"\n";
          throw runtime_error(os.str());
        }
      if (found1 == -1)
        {
          ostringstream os;
          os << "Selected normalization to lineshape not available: "<< 
            normalizationfactor[k] <<"\n";
          throw runtime_error(os.str());
        }
 
      // now set the lineshape variables 
      lineshape[k].SetInd_ls( found0 );
      lineshape[k].SetInd_lsn( found1 );
      lineshape[k].SetCutoff( cutoff[k] );
    }
}
 
 
 
void raw_vmrsReadFromScenario(// WS Output:
                              ArrayOfMatrix&   raw_vmrs,
                              // WS Input:
                              const TagGroups&     tgs,
                              // Control Parameters:
                               const String&    basename)
{
  // The species lookup data:
  extern const Array<SpeciesRecord> species_data;
 
  // We need to read one profile for each tag group.
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      // Determine the name.
      String name =
        basename + "." +
        species_data[tgs[i][0].Species()].Name() + ".aa";
      
      // Add an element for this tag group to the vmr profiles:
      raw_vmrs.push_back(Matrix());
      
      // Read the VMR:
      // (We use the workspace method MatrixReadAscii for this.)
      MatrixReadAscii(raw_vmrs[i],"",name);
      
      // state the source of profile.
      out3 << "  " << species_data[tgs[i][0].Species()].Name()
           << " profile read from file: " << name << "\n";
    }
}
 
void raw_vmrsReadFromFiles(// WS Output:
                           ArrayOfMatrix&   raw_vmrs,
                           // WS Input:
                           const TagGroups& tgs,
                           // Control Parameters:
                           const ArrayOfString&  seltags,
                           const ArrayOfString&  filenames,
                           const String&         basename)
{
  // The species lookup data:
  extern const Array<SpeciesRecord> species_data;
  ArrayOfIndex i_seltags;
 
  // Get indices of seltags in tgs. The function will throw an error
  // unless each tg in seltags corresponds exactly to a tg
  // in tgs. 
  get_tagindex_for_Strings( i_seltags, tgs, seltags );
 
  // Now we have to build an Array of the filenames to use for all the
  // tag groups. Initialize it with basename...
  ArrayOfString true_filenames(tgs.nelem());
  true_filenames = basename + '.';
 
  // Initialize to the standard scenario (see function
  // raw_vmrsReadFromScenario): 
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
     true_filenames[i] +=
        species_data[tgs[i][0].Species()].Name() + ".aa";
     // Should be identical to how the filenames are constructed in
     // raw_vmrsReadFromScenario! 
    }
 
  // Replace the names of the tgs given in seltags with the names in
  // filenames: 
  for ( Index i=0; i<seltags.nelem(); ++i )
    {
      true_filenames[i_seltags[i]] = filenames[i];
    }
    
  // Read the files:
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      // Add an element for this tag group to the vmr profiles:
      raw_vmrs.push_back(Matrix());
      
      // Read the VMR:
      // (We use the workspace method MatrixReadAscii for this.)
      MatrixReadAscii(raw_vmrs[i],"",true_filenames[i]);
      
      // state the source of profile.
      out3 << "  " << species_data[tgs[i][0].Species()].Name()
           << " profile read from file: " << true_filenames[i] << "\n";
    }
}
 
void WaterVaporSaturationInClouds( // WS Input/Output
                                  Matrix&           vmrs,  // manipulates this WS
                                  Vector&           p_abs, // manipulates this WS
                                  // WS Input
                                  const Vector&     t_abs, // constant
                                  const TagGroups&  tgs  ) // constant
{
 
  // make sure that the VMR and pressure grid are the same
  assert ( vmrs.ncols() == p_abs.nelem() );
 
  // The species lookup data
  extern const Array<SpeciesRecord> species_data;
  // cloud tag numbers:
  Index liquid_index = 1+tgs.nelem();
  Index ice_index    = 1+tgs.nelem();
  Index h2o_index[tgs.nelem()];
 
  // check size of input vectors.
  assert ( t_abs.nelem() == p_abs.nelem() ); 
  assert ( vmrs.ncols()  == p_abs.nelem() ); 
 
  // find tags for clouds and for water vapor
  Index u = 0;
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      String tag_name = species_data[tgs[i][0].Species()].Name();
      if (tag_name == "liquidcloud") // <=== liquid water clouds tag name
        {
          liquid_index = i;
        }
      if (tag_name == "icecloud") // <====== ice water clouds tag name
        {
          ice_index = i;
        }
      if (tag_name == "H2O") // <=========== water vapor tags to change VMR
        {
          h2o_index[u++] = i;
          //cout << "tag_name=" << tag_name << ",  tag=" << i << ",  u=" << u 
          //     << ",  h2o_index[u]=" << h2o_index[u] << "\n";
        }
    }
 
 
  // if no water vapor profile is in use do not go further
  if (u < 1)
    {
      out2 << "WaterVaporSaturationInClouds: no H2O profile found to adjust for clouds.\n"
           << "Therefore no saturation calculation is performed\n";
      return;
    }
 
  // ------------------------< saturation over liquid water >------------------------
  // modify the water vapor profiles for saturation over liquid water in a liquid water cloud
  if ( (liquid_index >= 0) && (liquid_index < tgs.nelem()) )
    {
      // sauration over liquid water 
      for (Index uu=0; uu<u; ++uu)              // --- loop over all H2O tags
        {
          for (Index i=0; i<vmrs.ncols() ; ++i) // --- loop over altitude grid
            {
              if (vmrs(liquid_index,i) > 0.000) // --- cloud present or not?
                {
                  // dry air part of the total pressure
                  Numeric p_dry         =  p_abs[i] * ( 1.000e0 - vmrs(h2o_index[uu],i) );
                  // water vapor saturation over liquid water 
                  Numeric e_s           =  WVSatPressureLiquidWater( t_abs[i] );
                  // new total pressure = dry air pressure +  water vapor saturation pressure
                  p_abs[i]              =  e_s + p_dry;
                  // new water vapor volume mixing ratio
                  vmrs(h2o_index[uu],i) =  e_s / p_abs[i];
                  // check if VMR has strange value
                  if ( (vmrs(h2o_index[uu],i) < 0.000e0) || (vmrs(h2o_index[uu],i) > 1.000e0) ) 
                    {
                      ostringstream os;
                      os << "WaterVaporSaturationInClouds: The water vapor VMR value " 
                         << vmrs(h2o_index[uu],i) << "\n"
                         << " looks strange after setting it to the saturation pressure over liquid water.";
                      throw runtime_error(os.str());
                      return;
                    }
                }
            }
        }
    }
 
 
  // ------------------------< saturation over ice water >------------------------
  // modify the water vapor profiles for saturation over ice water in a ice water cloud
  if ( (ice_index >= 0) && (ice_index < tgs.nelem()) )
    {
      for (Index uu=0; uu<u; ++uu)              // --- loop over all H2O tags
        {
          for (Index i=0; i<vmrs.ncols() ; ++i) // --- loop over altitude grid
            {
              if (vmrs(ice_index,i) > 0.000)    // --- cloud present or not?
                {
                  // dry air part of the total pressure
                  Numeric p_dry         =  p_abs[i] * ( 1.000e0 - vmrs(h2o_index[uu],i) );
                  // water vapor saturation over ice
                  Numeric e_s           =  WVSatPressureIce( t_abs[i] );
                  // new total pressure = dry air pressure +  water vapor saturation pressure
                  p_abs[i]              =  e_s + p_dry;
                  // new water vapor volume mixing ratio
                  vmrs(h2o_index[uu],i) =  e_s / p_abs[i];
                  // check if VMR has strange value
                  if ( (vmrs(h2o_index[uu],i) < 0.000e0)  || (vmrs(h2o_index[uu],i) > 1.000e0) )
                    {
                      ostringstream os;
                      os << "WaterVaporSaturationInClouds: The water vapor VMR value " 
                         << vmrs(h2o_index[uu],i) << "\n"
                         << " looks strange after setting it to the saturation pressure over ice.";
                      throw runtime_error(os.str());
                      return;
                    }
                }
            }
        }
    }
  
  return;
}
 
 
 
void AtmFromRaw(// WS Output:
                  Vector&        t_abs,
                  Vector&        z_abs,
                  Matrix&        vmrs,
                  // WS Input:      
                  const TagGroups&       tgs,
                  const Vector&          p_abs,
                  const Matrix&          raw_ptz,
                  const ArrayOfMatrix&   raw_vmrs)
{
  
  //---------------< 1. Interpolation of temperature and altitude >---------------
  {  
    // Safety check: Make sure that raw_ptz really is a [x,3] matrix:
    if ( 3 != raw_ptz.ncols() )
      {
        ostringstream os;
        os << "The variable raw_ptz does not have the right dimensions,\n"
           << "ncols() should be 3, but is actually "<< raw_ptz.ncols();
        throw runtime_error(os.str());
      }
 
 
    // Contents of raw_ptz: p_raw is column 1,  tz_raw is column 2-3.
 
    // Interpolate tz_raw to p_abs grid.
    // The reason why we take tz_raw as a matrix is that the
    // interpolation can then be done simultaneously, hence slightly
    // more efficient.
 
    // For the interpolated profiles:
    Matrix tz_intp( 2, p_abs.nelem() );
 
    interpp( tz_intp,
             raw_ptz(Range(joker),0),
             transpose(raw_ptz(Range(joker),Range(1,joker))),
             p_abs );
 
    // The first Matpack expression selects the first column of
    // raw_ptz as a vector. The second Matpack expression gives the
    // transpose of the last two columns of raw_ptz. The function
    // interpp can be called with these selections directly.
 
    // Extract t_abs:
    t_abs.resize( tz_intp.ncols() );
    t_abs = tz_intp(0,Range(joker));    // Matpack can copy the first row of
                                        // tz_intp to t_abs like this. But
                                        // t_abs has to have the right size!
 
    // Extract z_abs:
    z_abs.resize( tz_intp.ncols() );
    z_abs = tz_intp(1,Range(joker));    // Matpack can copy the second row of
                                        // tz_intp to t_abs like this. But
                                        // t_abs has to have the right size!
  }
 
  //---------------< 2. Interpolation of VMR profiles >---------------
  {
    // check size of input String vectors.
    assert ( tgs.nelem() == raw_vmrs.nelem() ); 
 
    // Make vmrs the right size:
    vmrs.resize( raw_vmrs.nelem(), p_abs.nelem() );
    
    // For sure, we need to loop through all VMR profiles:
    for ( Index j=0; j<raw_vmrs.nelem(); ++j )
      {
 
        // Get a reference to the profile we are concerned with:
        const Matrix& raw = raw_vmrs[j];
 
        // Raw should be a matrix with dimension [x,2], the first column
        // is the raw pressure grid, the second column the VMR values.
      
        // Safety check to ensure this:
        if ( 2 != raw.ncols() )
          {
            ostringstream os;
            os << "The variable raw_vmrs("
               << j
               << ") does not have the right dimensions,\n"
               << "ncols() should be 2, but is actually "<< raw.ncols();
            throw runtime_error(os.str());
          }
 
        // Interpolate the profile to the predefined pressure grid:
        //String tag_name = species_data[tgs[j][0].Species()].Name(); // name of the tag
        // if ( (tag_name == "liquidcloud") || (tag_name == "icecloud") )
        //   {
        //     // Interpolate linearly the cloud profiles
        //     interpp_cloud( vmrs(j,Range(joker)),
        //                 raw(Range(joker),0),
        //                 raw(Range(joker),1),
        //                 p_abs );
            //  out3 << "This VMR: " << vmrs(j,Range(joker)) << "\n";
        //   }
        // else
        //   {
            // Interpolate VMRs:
            interpp( vmrs(j,Range(joker)),
                     raw(Range(joker),0),
                     raw(Range(joker),1),
                     p_abs );
            // out3 << "This VMR: " << vmrs(j,Range(joker)) << "\n";
        //   }
        // The calls to interpp_cloud and inerpp contain some nice
        // Matpack  features:
        // 1. vmrs(j,Range(joker)) selects the jth row of vmrs.
        // 2. raw(Range(joker),0) and raw(Range(joker),1) select the
        // first and second column of raw. We don't need transpose
        // here, since the selected objects are vectors. 
        //
        // Note that you can call the interpolation functions directly
        // with these selections. 
      }
  }
}
 
 
void hseSet(
          Vector&    hse,
    const Numeric&   pref,
    const Numeric&   zref,
    const Numeric&   g0,
    const Index&     niter )
{
  hse.resize( 5 );
  
  hse[0] = 1;
  hse[1] = pref;
  hse[2] = zref;
  hse[3] = g0;
  hse[4] = Numeric( niter );
}
 
void hseSetFromLatitude(
                        Vector&    hse,
                const Numeric&   pref,
                const Numeric&   zref,
                const Numeric&   latitude,
                const Index&     niter )
{
 
  hse.resize( 5 );
  
  hse[0] = 1;
  hse[1] = pref;
  hse[2] = zref;
  hse[3] = g_of_lat(latitude);
  hse[4] = Numeric( niter );
}
 
void hseSetFromLatitudeIndex(
                              Vector&    hse,
                      const Vector&    p_abs,
                      const Vector&    z_abs,
                      const Numeric&   latitude,
                      const Index&     index,
                      const Index&     niter )
{
 
  /* check index range */
  check_if_in_range( 0, (Numeric)p_abs.nelem()-1, (Numeric)index, "index" );
 
  hse.resize( 5 );
  
  hse[0] = 1;
  hse[1] = p_abs[index];
  hse[2] = z_abs[index];
  hse[3] = g_of_lat(latitude);
  hse[4] = Numeric( niter );
}
 
 
void hseFromBottom(
          Vector&    hse,
    const Vector&    p_abs,
    const Vector&    z_abs,
    const Numeric&   g0,
    const Index&       niter )
{
  hseSet( hse, p_abs[0], z_abs[0], g0, niter );
}
 
 
 
void hseOff(
          Vector&    hse )
{
  hse.resize( 1 );
  hse[0] = 0;
}
 
 
 
// Algorithm based on equations from my (PE) thesis (page 274) and the book
// Meteorology today for scientists and engineers by R.B. Stull (pages 9-10). 
//
void hseCalc(
          Vector&    z_abs,
    const Vector&    p_abs,
    const Vector&    t_abs,
    const Vector&    h2o_abs,
    const Numeric&   r_geoid,   
    const Vector&    hse )
{
  check_if_bool( static_cast<Index>(hse[0]), 
                                        "the HSE flag (first element of hse)");
  
  if ( hse[0] )
  {
    if ( hse.nelem() != 5 )
      throw runtime_error("The length of the *hse* vector must be 5.");
 
    const Index   np = p_abs.nelem();
          Index   i;                     // altitude index
          Numeric  g;                     // gravitational acceleration
          Numeric  r;                     // water mixing ratio in gram/gram
          Numeric  tv;                    // virtual temperature
          Numeric  dz;                    // step geometrical altitude
          Vector   ztmp(np);              // temporary storage for z_abs
  
    // Pick out values from hse
    const Numeric   pref  = hse[1];
    const Numeric   zref  = hse[2];
    const Numeric   g0    = hse[3];
    const Index     niter = Index( hse[4] );
  
    check_lengths( z_abs, "z_abs", t_abs, "t_abs" );  
    check_lengths( z_abs, "z_abs", h2o_abs, "h2o_abs" );  
    if ( niter < 1 )
      throw runtime_error("The number of iterations must be > 0.");
  
    for ( Index iter=0; iter<niter; iter++ )
    {
      // Init ztmp
      ztmp[0] = z_abs[0];
  
      // Calculate new altitudes (relative z_abs(1)) and store in ztmp
      for ( i=0; i<np-1; i++ )
      {
        // Calculate g 
        g  = ( g_of_z(r_geoid,g0,z_abs[i]) + 
               g_of_z(r_geoid,g0,z_abs[i+1]) ) / 2.0;
  
        // Calculate weight mixing ratio for water assuming constant average
        // molecular weight of the air
        r  = 18/28.96 * (h2o_abs[i]+h2o_abs[i+1])/2;
  
        // The virtual temperature (no liquid water)
        tv = (1+0.61*r) * (t_abs[i]+t_abs[i+1])/2;
  
        // The change in vertical altitude from i to i+1 
        dz = 287.053*tv/g * log( p_abs[i]/p_abs[i+1] );
        ztmp[i+1] = ztmp[i] + dz;
      }
  
      // Match the altitude of the reference point
      dz = interpp( p_abs, ztmp, pref ) - zref;
 
      //  z_abs = ztmp - dz;
      z_abs = ztmp;
      z_abs -= dz;              
    }
  }
}
 
 
void h2o_absSet(
                Vector&          h2o_abs,
                const TagGroups& tgs,
                const Matrix&    vmrs )
{
  const Index   n = tgs.nelem();
  Index   found = -1;
  String  s;
 
  for( Index i=0; i<n && found<0; i++ ) 
  {
    s = tgs[i][0].Name();
 
    if ( s.substr(0,4) == "H2O-" )
      found = i;
  }
 
  /* ----- original version -------------------------------------------
 
  if ( found < 0 )
    throw runtime_error("h2o_absSet: No tag group contains water!");
  
  h2o_abs.resize( vmrs.ncols() );
  h2o_abs = vmrs(found,Range(joker));   
  // Matpack can copy the contents of vectors like this. The
  // dimensions must be the same! The expression
  // vmrs(found,Range(joker)) selects the row with index corresponding
  // to found.
 
  --------------------------------------------------------------------- 
  */
 
  h2o_abs.resize( vmrs.ncols() );
  if ( found >= 0 )
    {
      h2o_abs = vmrs(found,Range(joker));       
    } else {
      out2 << "  WARNING in h2o_absSet: could not find any H2O tag. So H2O vmr is set to zero!\n";
      for( Index i=0; i<h2o_abs.nelem(); i++ ) h2o_abs[i] = 0.0000000000e0;
    }
 
 
 
 
}
 
 
 
 
void vmrsScale(
               Matrix&                vmrs,
               const TagGroups&       tgs,
               const ArrayOfString&   scaltgs,
               const Vector&          scalfac)
{
  Index                            itag;
  ArrayOfIndex                     tagindex;      
 
  if ( scalfac.nelem() != scaltgs.nelem()  )
    throw runtime_error("vmrScale: Number of tgs and fac are different!");
  
  get_tagindex_for_Strings( tagindex, tgs, scaltgs );
 
  const Index   n = tagindex.nelem();
 
  for ( itag=0; itag<n; itag++ )
    {
      //out2 << scalfac[itag] << ".\n";
      vmrs(tagindex[itag],Range(joker)) *= scalfac[itag]; 
      // Matpack can multiply all elements of a vector with a constant
      // factor like this. In this case the vector is the selected row
      // of Matrix vmrs.
  
      //out2 << vmrs(tagindex[itag],Range(joker)) << ".\n";
    }
}
 
 
 
void n2_absSet(
               Vector&            n2_abs,
               const   TagGroups& tgs,
               const   Matrix&    vmrs )
{
  const Index   n = tgs.nelem();
  Index     found = -1;
  String  s;
 
  for( Index i=0; i<n && found<0; i++ ) 
  {
    s = tgs[i][0].Name();
    
    if ( s.substr(0,3) == "N2-" )
      found = i;
  }
 
  /* ----- original version -------------------------------------------
 
  if ( found < 0 )
    throw runtime_error("n2_absSet: No tag group contains nitrogen!");
  
  n2_abs.resize( vmrs.ncols() );
  n2_abs = vmrs(found,Range(joker));
  // Matpack can copy the contents of vectors like this. The
  // dimensions must be the same! The expression
  // vmrs(found,Range(joker)) selects the row with index corresponding
  // to found.
 
  --------------------------------------------------------------------- 
  */
 
  n2_abs.resize( vmrs.ncols() );
    if ( found >= 0 )
    {
      n2_abs = vmrs(found,Range(joker));        
    } else {
      out2 << "  WARNING in n2_absSet: could not find any N2 tag. So N2 vmr is set to zero!\n";
      for( Index i=0; i<n2_abs.nelem(); i++ ) n2_abs[i] = 0.0000000000e0;
    }
}
 
 
void absCalc(// WS Output:
             Matrix&                         abs,
             ArrayOfMatrix&                  abs_per_tg,
             // WS Input:                 
             const TagGroups&                tgs,
             const Vector&                   f_mono,
             const Vector&                   p_abs,
             const Vector&                   t_abs,
             const Vector&                   n2_abs,
             const Vector&                   h2o_abs,
             const Matrix&                   vmrs,
             const ArrayOfArrayOfLineRecord& lines_per_tg,
             const ArrayOfLineshapeSpec&     lineshape,
             const ArrayOfString&            cont_description_names,
             const ArrayOfString&            cont_description_models,
             const ArrayOfVector&            cont_description_parameters)
{
  // Dimension checks are performed in the executed functions
 
  // allocate local variable to hold the cross sections per tag group
  ArrayOfMatrix xsec_per_tg;
  
  xsec_per_tgInit( xsec_per_tg, tgs, f_mono, p_abs );
 
  xsec_per_tgAddLines( xsec_per_tg,
                       tgs,
                       f_mono,
                       p_abs,
                       t_abs,
                       h2o_abs,
                       vmrs,
                       lines_per_tg,
                       lineshape );
 
  xsec_per_tgAddConts( xsec_per_tg,
                       tgs,
                       f_mono,
                       p_abs,
                       t_abs,
                       n2_abs,
                       h2o_abs,
                       vmrs,
                       cont_description_names,
                       cont_description_parameters,
                       cont_description_models);
 
  absCalcFromXsec(abs,
                  abs_per_tg,
                  xsec_per_tg,
                  vmrs );
 
}
 
void absCalcSaveMemory(// WS Output:
             Matrix&                         abs,
             // WS Input:                 
             const TagGroups&                tgs,
             const Vector&                   f_mono,
             const Vector&                   p_abs,
             const Vector&                   t_abs,
             const Vector&                   n2_abs,
             const Vector&                   h2o_abs,
             const Matrix&                   vmrs,
             const ArrayOfArrayOfLineRecord& lines_per_tg,
             const ArrayOfLineshapeSpec&     lineshape,
             const ArrayOfString&            cont_description_names,
             const ArrayOfString&            cont_description_models,
             const ArrayOfVector&            cont_description_parameters)
{
  // Dimension checks are performed in the executed functions
 
  // Allocate local variable to hold the cross sections per tag group:
  ArrayOfMatrix xsec_per_tg;
 
  // Allocate local variable to hold the absorption for each tag group:
  Matrix this_abs;
 
  // Allocate local variable to hold abs_per_tg for each tag
  // group. This is just a dummy, we need it, since it is a formal
  // argument of absCalcFromXsec.
  ArrayOfMatrix this_abs_per_tg;
 
  // Define variable to hold a dummy list of tag groups with only 1 element:
  TagGroups this_tg(1);
 
  // Local list of VMRs, only 1 element:
  Matrix this_vmr(1,p_abs.nelem());
 
  // Local lines_per_tg, only 1 element:
  ArrayOfArrayOfLineRecord these_lines(1);
 
  // Local lineshape list, only 1 element:
  ArrayOfLineshapeSpec     this_lineshape(1);
 
  // Initialize the output variable abs:
  abs.resize( f_mono.nelem(), p_abs.nelem() );
  abs = 0;                      // Matpack can set all elements like this.
 
  out2 << "  Number of tag groups to do: " << tgs.nelem() << "\n";
 
  // Loop over all species:
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      out2 << "  Doing tag group " << i << ".\n";
 
      // Get a dummy list of tag groups with only the current element:
      this_tg[0].resize(tgs[i].nelem());
      this_tg[0] = tgs[i];
 
      // VMR for this species:
      this_vmr(0,joker) = vmrs(i,joker);
 
      // List of lines:
      these_lines[0].resize(lines_per_tg[i].nelem());
      these_lines[0] = lines_per_tg[i];
 
      // List of lineshapes:
      this_lineshape[0] = lineshape[i];
 
      xsec_per_tgInit( xsec_per_tg, this_tg, f_mono, p_abs );
 
      xsec_per_tgAddLines( xsec_per_tg,
                           this_tg,
                           f_mono,
                           p_abs,
                           t_abs,
                           h2o_abs,
                           this_vmr,
                           these_lines,
                           this_lineshape );
 
      xsec_per_tgAddConts( xsec_per_tg,
                           this_tg,
                           f_mono,
                           p_abs,
                           t_abs,
                           n2_abs,
                           h2o_abs,
                           this_vmr,
                           cont_description_names,
                           cont_description_parameters,
                           cont_description_models);
 
      absCalcFromXsec(this_abs,
                      this_abs_per_tg,
                      xsec_per_tg,
                      this_vmr );
 
      // Add absorption of this species to total absorption:
      assert(abs.nrows()==this_abs.nrows());
      assert(abs.ncols()==this_abs.ncols());
      abs += this_abs;
    }
}
 
 
void absCalcFromXsec(// WS Output:
                     Matrix&                         abs,
                     ArrayOfMatrix&                  abs_per_tg,
                     // WS Input:                 
                     const ArrayOfMatrix&            xsec_per_tg,
                     const Matrix&                   vmrs)
{
  // Check that vmrs and xsec_per_tg really have compatible
  // dimensions. In vmrs there should be one row for each tg:
  if ( vmrs.nrows() != xsec_per_tg.nelem() )
    {
      ostringstream os;
      os << "Variable vmrs must have compatible dimension to xsec_per_tg.\n"
         << "vmrs.nrows() = " << vmrs.nrows() << '\n'
         << "xsec_per_tg.nelem() = " << xsec_per_tg.nelem();
      throw runtime_error(os.str());
    }
 
  // Check that number of altitudes are compatible. We only check the
  // first element, this is possilble because within arts all elements
  // are on the same altitude grid.
  if ( vmrs.ncols() != xsec_per_tg[0].ncols() )
    {
      ostringstream os;
      os << "Variable vmrs must have same numbers of altitudes as xsec_per_tg.\n"
         << "vmrs.ncols() = " << vmrs.ncols() << '\n'
         << "xsec_per_tg[0].ncols() = " << xsec_per_tg[0].ncols();
      throw runtime_error(os.str());
    }  
  
  // Initialize abs and abs_per_tg. The array dimension of abs_per_tg
  // is the same as that of xsec_per_tg. The dimension of abs should
  // be equal to one of the xsec_per_tg enries.
  abs.resize( xsec_per_tg[0].nrows(), xsec_per_tg[0].ncols() );
  abs = 0;                      // Matpack can set all elements like this.
 
  abs_per_tg.resize( xsec_per_tg.nelem() );
 
  out2 << "  Computing abs and abs_per_tg from xsec_per_tg.\n";
 
  // Loop through all tag groups
  for ( Index i=0; i<xsec_per_tg.nelem(); ++i )
    {
      out2 << "  Tag group " << i << '\n';
 
      // Make this element of xsec_per_tg the right size:
      abs_per_tg[i].resize( xsec_per_tg[i].nrows(), xsec_per_tg[i].ncols() );
      abs_per_tg[i] = 0;        // Initialize all elements to 0.
 
      // Loop through all altitudes
      for ( Index j=0; j<xsec_per_tg[i].ncols(); j++)
        {
          // Loop through all frequencies
          for ( Index k=0; k<xsec_per_tg[i].nrows(); k++)
            {
              abs_per_tg[i](k,j) = xsec_per_tg[i](k,j) * vmrs(i,j);
            }
        }
 
      // Add up to the total absorption:
      abs += abs_per_tg[i];     // In Matpack you can use the +=
                                // operator to do elementwise addition.
    }
}
 
 
void xsec_per_tgInit(// WS Output:
                     ArrayOfMatrix&   xsec_per_tg,
                     // WS Input:
                     const TagGroups& tgs,
                     const Vector&    f_mono,
                     const Vector&    p_abs
                     )
{
  // Initialize xsec_per_tg. The array dimension of xsec_per_tg
  // is the same as that of lines_per_tg.
  xsec_per_tg.resize( tgs.nelem() );
 
  // Loop xsec_per_tg and make each matrix the right size,
  // initializing to zero:
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      // Make this element of abs_per_tg the right size:
      xsec_per_tg[i].resize( f_mono.nelem(), p_abs.nelem() );
      xsec_per_tg[i] = 0;       // Matpack can set all elements like this.
    }
 
  out2 << "  Initialized xsec_per_tg.\n"
       << "  Number of frequencies        : " << f_mono.nelem() << "\n"
       << "  Number of pressure levels    : " << p_abs.nelem() << "\n";
}
 
void xsec_per_tgAddLines(// WS Output:
                         ArrayOfMatrix&                   xsec_per_tg,
                         // WS Input:             
                         const TagGroups&                 tgs,
                         const Vector&                    f_mono,
                         const Vector&                    p_abs,
                         const Vector&                    t_abs,
                         const Vector&                    h2o_abs,
                         const Matrix&                    vmrs,
                         const ArrayOfArrayOfLineRecord&  lines_per_tg,
                         const ArrayOfLineshapeSpec&      lineshape)
{
  // Check that all paramters that should have the number of tag
  // groups as a dimension are consistent.
  {
    const Index n_tgs    = tgs.nelem();
    const Index n_xsec   = xsec_per_tg.nelem();
    const Index n_vmrs   = vmrs.nrows();
    const Index n_lines  = lines_per_tg.nelem();
    const Index n_shapes = lineshape.nelem();
 
    if ( n_tgs != n_xsec  ||
         n_tgs != n_vmrs  ||
         n_tgs != n_lines ||
         n_tgs != n_shapes   )
      {
        ostringstream os;
        os << "The following variables must all have the same dimension:\n"
           << "tgs:          " << tgs.nelem() << '\n'
           << "xsec_per_tg:  " << xsec_per_tg.nelem() << '\n'
           << "vmrs:         " << vmrs.nrows() << '\n'
           << "lines_per_tg: " << lines_per_tg.nelem() << '\n'
           << "lineshape:    " << lineshape.nelem();
        throw runtime_error(os.str());
      }
  }  
 
  // Print information:
  //
  out2 << "  Calculating line spectra.\n";
  //
  // Uncomment the part below if you temporarily want detailed info about 
  // transitions to be done
 
      // The variables defined here (in particular the frequency
      // conversion) are just to make the output nice. They are not used
      // in subsequent calculations.
  //    cout << "  Transitions to do: \n";
  //    Index nlines = 0;
  //    String funit;
  //    Numeric ffac;
  //    if ( f_mono[0] < 3e12 )
  //      {
  //        funit = "GHz"; ffac = 1e9;
  //      }
  //    else
  //      {
  //        extern const Numeric SPEED_OF_LIGHT;
  //        funit = "cm-1"; ffac = SPEED_OF_LIGHT*100;
  //      }
  //    for ( Index i=0; i<lines_per_tg.nelem(); ++i )
  //      {
  //        for ( Index l=0; l<lines_per_tg[i].nelem(); ++l )
  //          {
  //          cout << "    " << lines_per_tg[i][l].Name() << " @ " 
  //               << lines_per_tg[i][l].F()/ffac  << " " << funit << " ("
  //               << lines_per_tg[i][l].I0() << "  "
  //               << lines_per_tg[i][l].Agam() << ")\n"; 
  //          nlines++;
  //        }
  //    }
  //  out2 << "  Total number of transistions : " << nlines << "\n";
 
  // Call xsec_species for each tag group.
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      out2 << "  Tag group " << i
           << " (" << get_tag_group_name(tgs[i]) << "): ";
      
      // Get a pointer to the line list for the current species. This
      // is just so that we don't have to type lines_per_tg[i] over
      // and over again.
      const ArrayOfLineRecord& ll = lines_per_tg[i];
 
      // Also get a pointer to the lineshape specification:
      const LineshapeSpec& ls = lineshape[i];
      
      // Skip the call to xsec_per_tg if the line list is empty.
      if ( 0 < ll.nelem() )
        {
 
          // Get the name of the species. The member function name of a
          // LineRecord returns the full name (species + isotope). So
          // it is for us more convenient to get the species index
          // from the LineRecord member function Species(), and then
          // use this to look up the species name in species_data.
          extern const Array<SpeciesRecord> species_data;
          String species_name = species_data[ll[0].Species()].Name();
 
          // Get the name of the lineshape. For that we use the member
          // function Ind_ls() to the lineshape data ls, which returns
          // an index. With that index we can go into lineshape_data
          // to get the name.
          // We will need this for safety checks later on.
          extern const Array<LineshapeRecord> lineshape_data;
          String lineshape_name = lineshape_data[ls.Ind_ls()].Name();
 
 
          // The use of overlap parameters for oxygen lines only makes
          // sense if the special Rosenkranz lineshape is used
          // (Rosenkranz_Voigt_Drayson or Rosenkranz_Voigt_Kuntz6). 
          // Likewise, the use of these lineshapes only makes sense if
          // overlap parameters are available. We will test both these
          // conditions.
 
          // First of all, let's find out if the species we are
          // dealing with is oxygen. 
          if ( "O2" == species_name )
            {
              // Do we have overlap parameters in the aux fields of
              // the LineRecord?
              if ( 0 != ll[0].Naux() )
                {
                  // Yes. So let's make sure that the lineshape is one
                  // that can use these parameters. 
                  if ( "Rosenkranz_Voigt_Drayson" != lineshape_name &&
                       "Rosenkranz_Voigt_Kuntz6"  != lineshape_name    )
                    {
                      ostringstream os;
                      os 
                        << "You are using a line catalogue that contains auxiliary parameters to\n"
                        << "take care of overlap for oxygen lines. But you are not using a\n"
                        << "lineshape that uses these parameters. Use Rosenkranz_Voigt_Drayson or\n"
                        << "Rosenkranz_Voigt_Kuntz6.";
                      throw runtime_error(os.str());                  
                    }
                }               
            }
 
          // Now we go the opposite way. Let's see if the Rosenkranz
          // lineshapes are used.
          if ( "Rosenkranz_Voigt_Drayson" == lineshape_name ||
               "Rosenkranz_Voigt_Kuntz6"  == lineshape_name    )
            {
              // Is the species oxygen, as it should be?
              if ( "O2" != species_name )
                {
                  ostringstream os;
                  os 
                    << "You are trying to use one of the special Rosenkranz lineshapes with\n"
                    << "overlap correction for a species other than oxygen. Your species is\n"
                    << species_name << ".\n"
                    << "Select another lineshape for this species.";
                    throw runtime_error(os.str());                    
                }
              else
                {
                  // Do we have overlap parameters in the aux fields of
                  // the LineRecord?
                  if ( 0 == ll[0].Naux() )
                    {
                      ostringstream os;
                      os 
                        << "You are trying to use one of the special Rosenkranz lineshapes with\n"
                        << "overlap correction. But your line file does not contain aux\n"
                        << "parameters. (I've checked only the first LineRecord). Use a line file\n"
                        << "with overlap parameters.";
                        throw runtime_error(os.str());                
                    }
                }
            }
 
          out2 << ll.nelem() << " transitions\n";
          xsec_species( xsec_per_tg[i],
                        f_mono,
                        p_abs,
                        t_abs,
                        h2o_abs,
                        vmrs(i,Range(joker)),
                        ll,
                        ls.Ind_ls(),
                        ls.Ind_lsn(),
                        ls.Cutoff());
          // Note that we call xsec_species with a row of vmrs,
          // selected by the above Matpack expression. This is
          // possible, because xsec_species is using Views.
        }
      else
        {
          out2 << ll.nelem() << " transitions, skipping\n";
        }
    }
}
 
void xsec_per_tgAddConts(// WS Output:
                         ArrayOfMatrix&                   xsec_per_tg,
                         // WS Input:             
                         const TagGroups&                 tgs,
                         const Vector&                    f_mono,
                         const Vector&                    p_abs,
                         const Vector&                    t_abs,
                         const Vector&                    n2_abs,
                         const Vector&                    h2o_abs,
                         const Matrix&                    vmrs,
                         const ArrayOfString&             cont_description_names,
                         const ArrayOfVector&             cont_description_parameters,
                         const ArrayOfString&             cont_description_models )
{
  // Check that all paramters that should have the number of tag
  // groups as a dimension are consistent.
  {
    const Index n_tgs    = tgs.nelem();
    const Index n_xsec   = xsec_per_tg.nelem();
    const Index n_vmrs   = vmrs.nrows();
 
    if ( n_tgs != n_xsec || n_tgs != n_vmrs )
      {
        ostringstream os;
        os << "The following variables must all have the same dimension:\n"
           << "tgs:          " << tgs.nelem() << '\n'
           << "xsec_per_tg:  " << xsec_per_tg.nelem() << '\n'
           << "vmrs:         " << vmrs.nrows();
        throw runtime_error(os.str());
      }
  }
 
  // Check, that dimensions of cont_description_names and
  // cont_description_parameters are consistent...
  if ( cont_description_names.nelem() !=
       cont_description_parameters.nelem() )
    {
      for (Index i=0; i < cont_description_names.nelem(); ++i) 
        {
          cout << "xsec_per_tgAddConts: " << i << " name : " << cont_description_names[i] << "\n";
        }
      for (Index i=0; i < cont_description_parameters.nelem(); ++i) 
        {
          cout << "xsec_per_tgAddConts: " << i << " param: " << cont_description_parameters[i] << "\n";
        }
      for (Index i=0; i < cont_description_models.nelem(); ++i) 
        {
          cout << "xsec_per_tgAddConts: " << i << " option: " << cont_description_models[i] << "\n";
        }
      ostringstream os;
        os << "The following variables must have the same dimension:\n"
           << "cont_description_names:      " << cont_description_names.nelem() << '\n'
           << "cont_description_parameters: " << cont_description_parameters.nelem();
        throw runtime_error(os.str());
    }
 
  // ...and that indeed the names match valid continuum models:
  for ( Index i=0; i<cont_description_names.nelem(); ++i )
    {
      check_continuum_model(cont_description_names[i]);
    }
 
  out2 << "  Calculating continuum spectra.\n";
 
  // Loop tag groups:
  for ( Index i=0; i<tgs.nelem(); ++i )
    {
      extern const Array<SpeciesRecord> species_data; 
 
      // Go through the tags in the current tag group to see if they
      // are continuum tags:  
      for ( Index s=0; s<tgs[i].nelem(); ++s )
        {
          // First of all, we have to make sure that this is not a
          // tag that means `all isotopes', because this should not
          // include continuum. For such tags, tag.Isotope() will
          // return the number of isotopes (i.e., one more than the
          // allowed index range).
          if ( tgs[i][s].Isotope() <
               species_data[tgs[i][s].Species()].Isotope().nelem() )
            {
              // If we get here, it means that the tag describes a
              // specific isotope. Could be a continuum tag!
                
              // The if clause below checks whether the abundance of this tag
              // is negative. (Negative abundance marks continuum tags.)
              // It does the following:
              //
              // 1. species_data contains the lookup table of species
              //          specific data. We need the right entry in this
              //          table. The index of this is obtained by calling member function
              //          Species() on the tag. Thus we have:
              //          species_data[tgs[i][s].Species()].
              //
              // 2. Member function Isotope() on the above biest gives
              //    us the array of isotope specific data. This we have
              //    to subscribe with the right isotope index, which we
              //    get by calling member function Isotope on the
              //    tag. Thus we have:
              //    Isotope()[tgs[i][s].Isotope()]
              //
              // 3. Finally, from the isotope data we need to get the mixing
              //    ratio by calling the member function Abundance().
              if ( 0 >
                   species_data[tgs[i][s].Species()].Isotope()[tgs[i][s].Isotope()].Abundance() )
                {
                  // We have identified a continuum tag!
 
                  // Get only the continuum name. The full tag name is something like:
                  // H2O-HITRAN96Self-*-*. We want only the `H2O-HITRAN96Self' part:
                  const String name =
                    species_data[tgs[i][s].Species()].Name() + "-"
                    + species_data[tgs[i][s].Species()].Isotope()[tgs[i][s].Isotope()].Name();
  
                  // Check that this corresponds to a valid continuum
                  // model:
                  check_continuum_model(name);
 
                  // Check, if we have parameters for this model. For
                  // this, the model name must be listed in
                  // cont_description_names.
                  const Index n =
                    (Index)(find( cont_description_names.begin(),
                                  cont_description_names.end(),
                                  name ) - cont_description_names.begin());
 
                  // n==cont_description_names.nelem() indicates that
                  // the name was not found.
                  if ( n==cont_description_names.nelem() )
                    {
                      ostringstream os;
                      os << "Cannot find model " << name
                         << " in cont_description_names.";
                      throw runtime_error(os.str());                  
                    }
 
                  // Ok, the tag specifies a valid continuum model and
                  // we have continuum parameters.
                  
                  out2 << "  Adding " << name
                       << " to tag group " << i << ".\n";
 
                  // find the options for this continuum tag from the input array
                  // of options. The actual field of the array is n:
                  const String ContOption = cont_description_models[n];
 
                  // Add the continuum for this tag. The parameters in
                  // this call should be clear. The vmr is in
                  // vmrs(i,Range(joker)). The other vmr variable, `h2o_abs'
                  // contains the real H2O vmr, which is needed for
                  // the oxygen continuum.
                  xsec_continuum_tag( xsec_per_tg[i],
                                      name,
                                      cont_description_parameters[n],
                                      cont_description_models[n], 
                                      f_mono,
                                      p_abs,
                                      t_abs,
                                      n2_abs,
                                      h2o_abs,
                                      vmrs(i,Range(joker)) );
                  // Calling this function with a row of Matrix vmrs
                  // is possible because it uses Views.
                }
            }
        }
    }
 
}
 
void abs_per_tgReduce(// WS Output:
                      ArrayOfMatrix&         abs_per_tg,
                      // WS Input:
                      const TagGroups&       tgs,
                      const TagGroups&       wfs_tgs)
{
 
  // Make a safety check that the dimensions of tgs and
  // abs_per_tg are the same (could happen that we call this workspace
  // method twice by accident).
  if ( abs_per_tg.nelem()!=tgs.nelem() )
    throw(runtime_error("The variables abs_per_tg and tgs must\n"
                        "have the same dimension."));
 
  // Erase could be very inefficient in this case, since elements
  // behind the erased one are copied in order to fill the
  // gap. Therefore, we will construct a new abs_per_tg, and finally
  // use it to replace the old one.
  ArrayOfMatrix abs_per_tg_out( wfs_tgs.nelem() );
 
  // Go through the weighting function tag groups:
  for ( Index i=0; i<wfs_tgs.nelem(); ++i )
    {
      // Index to the elements of wfs_tgs in tgs:
      Index n;
      get_tag_group_index_for_tag_group( n, tgs, wfs_tgs[i] );
 
      abs_per_tg_out[i].resize( abs_per_tg[n].nrows(), abs_per_tg[n].ncols() );
      abs_per_tg_out[i] = abs_per_tg[n]; // Matpack can copy the contents of
                                         // matrices like this. The dimensions
                                         // must be the same! 
    }  
 
  // Copy the generated matrices back to abs_per_tg
  abs_per_tg.resize( wfs_tgs.nelem() );
  abs_per_tg = abs_per_tg_out;  // FIXME: It should be checked whether
                                // this works correctly. Does my Array
                                // implementation work as it should? 
}
 
 
 
//======================================================================
//             Methods related to refraction
//======================================================================
 
void refrSet( 
              Index&    refr,
              Index&    refr_lfac,
              String&   refr_model,
        const Index&    on,
        const String&   model,
        const Index&    lfac )
{
  check_if_bool( on, "on" );
  if ( lfac < 1 )
    throw runtime_error("The length factor cannot be smaller than 1.");      
 
  refr       = on;
  refr_lfac  = lfac;
  refr_model = model;
}
 
 
 
void refrOff( 
              Index&    refr,
              Index&    refr_lfac,
              String&   refr_model )
{
  refrSet( refr, refr_lfac, refr_model, 0, "", 1 );
  refr_lfac = 0;
}
 
 
 
void refrCalc (
                    Vector&   refr_index,
              const Vector&   p_abs,
              const Vector&   t_abs,
              const Vector&   h2o_abs,
              const Index&    refr,
              const String&   refr_model )
{
  check_if_bool( refr, "refr" );
 
  if ( refr == 0 )
    refr_index.resize( 0 );
 
  else
  {
    if ( refr_model == "Unity" )
    {
      const Index n = p_abs.nelem();
      refr_index.resize( n );
      refr_index = 1.0;
    }
    
    else if ( refr_model == "Boudouris" )
      refr_index_Boudouris( refr_index, p_abs, t_abs, h2o_abs );
 
    else if ( refr_model == "BoudourisDryAir" )
      refr_index_BoudourisDryAir( refr_index, p_abs, t_abs );
 
    else
    {
      ostringstream os;
      os << "Unknown refraction parameterization: " << refr_model << "\n"
         << "Existing parameterizations are: \n"
         << "Unity\n"
         << "Boudouris\n"
         << "BoudourisDryAir";
      throw runtime_error(os.str());
    }
  }
}
 
 
 
 
//======================================================================
//             Methods related to continua
//======================================================================
 
void cont_descriptionInit(// WS Output:
                          ArrayOfString& names,
                          ArrayOfString& options, 
                          ArrayOfVector& parameters)
{
  names.resize(0);
  options.resize(0);
  parameters.resize(0);
  out2 << "  Initialized cont_description_names \n"
          "              cont_description_models\n"
          "              cont_description_parameters.\n";
}
 
 
void cont_descriptionAppend(// WS Output:
                       ArrayOfString& cont_description_names,
                       ArrayOfString& cont_description_models,
                       ArrayOfVector& cont_description_parameters,
                       // Control Parameters:
                       const String& tagname,
                       const String& model,
                       const Vector& userparameters)
{
  // First we have to check that name matches a continuum species tag.
  check_continuum_model(tagname);
 
  //cout << "   + tagname:    " << tagname << "\n";
  //cout << "   + model:      " << model << "\n";
  //cout << "   + parameters: " << userparameters << "\n";
 
  // Add name and parameters to the apropriate variables:
  cont_description_names.push_back(tagname);
  cont_description_models.push_back(model);
  cont_description_parameters.push_back(userparameters);
}