diff --git a/ChangeLog b/ChangeLog
index dd89f84799d3a8f73d7859c0986ba34fb64f3578..a61ca2bcaff27d72db22c4670c8abd7c2ab44fa2 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,12 @@
+2015-11-17  Oliver Lemke  <olemke@core-dump.info>
+
+	* ARTS-XML-DATA-2-3-19
+	
+	* planets/Earth/ECMWF/IFS/Eresmaa_137L:
+
+	Added Eresmaa_137L dataset. It is a enhancement of the Chevallier
+	dataset. The ARTS XML files were created by John Mrziglod.
+
 2015-10-13  Richard Larsson  <ric.larsson@gmail.com>
 
 	* ARTS-XML-DATA-2-3-18
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/ConvertCompactAtm_ArrayOfMatrix-to-ArrayOfGriddedField4.arts b/planets/Earth/ECMWF/IFS/Eresmaa_137L/ConvertCompactAtm_ArrayOfMatrix-to-ArrayOfGriddedField4.arts
new file mode 100644
index 0000000000000000000000000000000000000000..64374cd2b5a5f1fca5edb0209544b2d9766c1e8f
--- /dev/null
+++ b/planets/Earth/ECMWF/IFS/Eresmaa_137L/ConvertCompactAtm_ArrayOfMatrix-to-ArrayOfGriddedField4.arts
@@ -0,0 +1,117 @@
+#DEFINITIONS:  -*-sh-*-
+
+Arts2 {
+    ## This script converts the xml output of 'extract_arts_l137.f90' to ARTS compatible data format ArrayOfGriddedField4
+
+    ## Call this script by running 
+    ## ./arts/build/src/arts arts-xml-data/planets/Earth/ECMWF/IFS/Eresmaa_137L/ConvertCompactAtm_ArrayOfMatrix-to-ArrayOfGriddedField4.arts
+
+    StringCreate( datapath )
+    StringSet( datapath, "arts-xml-data/planets/Earth/ECMWF/IFS/Eresmaa_137L/" )
+    ArrayOfStringCreate( filenames )
+    ArrayOfStringCreate( fieldnames )
+    ArrayOfStringCreate( filenames_extra )
+    ArrayOfStringCreate( fieldnames_extra )
+
+    # ECMWF data (Eresmaa and similar)
+    ArrayOfStringSet( filenames, [
+        "eresmaal137_all_ccol.xml",
+        "eresmaal137_all_oz.xml",
+        "eresmaal137_all_q.xml",
+        "eresmaal137_all_rcol.xml",
+        "eresmaal137_all_t.xml"
+    ] )
+    ArrayOfStringSet( fieldnames, [
+                     "T", "z",
+                     "scat_species-LWC-mass_density", "scat_species-IWC-mass_density",
+                     "scat_species-RR-mass_flux", "scat_species-SR-mass_flux",
+                     "abs_species-H2O", "abs_species-O3"
+    ] )
+    
+    ArrayOfStringSet( filenames_extra, [
+        "eresmaal137_all_ccol_extra.xml",
+        "eresmaal137_all_oz_extra.xml",
+        "eresmaal137_all_q_extra.xml",
+        "eresmaal137_all_rcol_extra.xml",
+        "eresmaal137_all_t_extra.xml"
+    ] )
+    ArrayOfStringSet( fieldnames_extra, [
+                     "cloud cover", "vertical velocity"
+    ] )
+
+
+
+
+    ###########################################
+
+    StringCreate( fileloc )
+    StringCreate( filename_output )
+    StringCreate( filename_input )
+    Copy( filename_output, datapath )
+    IndexCreate( ncases )
+    nelemGet( ncases, filenames )
+
+    # Some prep settings
+    #-------------------
+    output_file_formatSetZippedAscii
+
+    #output_file_formatSetAscii
+    IndexSet (atmosphere_dim, 1)
+    ArrayOfMatrixCreate( arrayofmatrix_1 )
+
+
+    ## atmospheric fields (required by ARTS)
+    AgendaSet( forloop_agenda ){
+      Extract( fileloc, filenames, forloop_index )
+      Append( filename_output, fileloc )
+      
+      # the temporary files of the fortran program are named 'FILENAME.xml.temp'
+      StringCompose( filename_input, filename_output, ".temp" )
+
+      # Read ArrayOfMatrix batch data
+      #------------------------------
+      ReadXML( arrayofmatrix_1, filename_input )
+
+
+      # Convert to batch_atm_fields_compact
+      # -----------------------------------
+      batch_atm_fields_compactFromArrayOfMatrix(
+        atmospheres_fields=arrayofmatrix_1,
+        field_names=fieldnames )
+    
+      
+      WriteXML( in=batch_atm_fields_compact, filename=filename_output )
+    }
+
+    IndexStepDown( ncases, ncases )
+    ForLoop( forloop_agenda, 0, ncases, 1 )
+    
+    ## additional atmospheric fields (cannot processed by ARTS directly)
+    nelemGet( ncases, filenames_extra )
+    
+    AgendaSet( forloop_agenda){
+      Extract( fileloc, filenames_extra, forloop_index )
+      Append( filename_output, fileloc )
+      
+      # the temporary files of the fortran program are named 'FILENAME.xml.temp'
+      StringCompose( filename_input, filename_output, ".temp" )
+
+      # Read ArrayOfMatrix batch data
+      #------------------------------
+      ReadXML( arrayofmatrix_1, filename_input )
+
+
+      # Convert to batch_atm_fields_compact
+      # -----------------------------------
+      batch_atm_fields_compactFromArrayOfMatrix(
+        atmospheres_fields=arrayofmatrix_1,
+        field_names=fieldnames_extra )
+    
+      
+      WriteXML( in=batch_atm_fields_compact, filename=filename_output )
+    }
+
+    IndexStepDown( ncases, ncases )
+    ForLoop( forloop_agenda, 0, ncases, 1 )
+
+} # End of Main
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/README b/planets/Earth/ECMWF/IFS/Eresmaa_137L/README
new file mode 100644
index 0000000000000000000000000000000000000000..f088336831fbbcd9bbf08f188a7494aa7ef49a5b
--- /dev/null
+++ b/planets/Earth/ECMWF/IFS/Eresmaa_137L/README
@@ -0,0 +1,73 @@
+Usage
+-----
+The ready-to-use arts xml files (formatted in ArrayOfGriddedField4) are stored under this path. Please use only the "*.xml.gz" files and ignore potential "*.temp" files.
+
+OR
+
+If you want to convert the ECMWF profile datasets by yourself to the ArrayOfGriddedField4 format, please follow these steps:
+1) Download the original files from the NWP SAF wep page
+2) Extract the downloaded file 'profiles137.tar.gz' to this path.
+3) Compile "extract_arts_l137.f90" with your favourite FORTRAN compiler (e.g. f95)
+4) Run the compiled output file
+5) Run "ConvertCompactAtm_ArrayOfMatrix-to-ArrayOfGriddedField4.arts" with your arts distribution
+
+
+Files
+-----
+*.xml.gz [ARTS (gridded field) format] - Contains a matrix of the atmospheric profiles (row: pressure level, col: atmospheric variables):
+    1) pressure [Pa] 
+    2) Temperature [K] 
+    3) Height [m]
+    4) CLW [kg/m^3]
+    5) CIW[kg/m^3]
+    6) Rain[kg/(m2*s)]
+    7) Snow[kg/(m2*s)]
+    8) H2O_VMR[1]
+    9) O3_VMR[1]
+    
+*_surface.xml.gz - contains the profile index, surface variables and geographic location:
+    1) Profile index [1-25000]
+    2) Ln of Surf Pressure in [Pa]
+    3) Surface geopotential [m2/s2] 
+    4) Surface Skin Temperature [K]
+    5) 2m Temperature [K]
+    6) 2m Dew point temperature [K]
+    7) 10m wind speed [m/s]
+    8) 10m wind direction [0-360 DEG]
+    9) Stratiform precipitation at surface [m]
+    10) Convective precipitation at surface [m]
+    11) Snowfall at surface [m] (water equival.)
+    12) Land/sea Mask [0-1]
+    13) Latitude [deg]
+    14) Longitude [deg]
+    15) Year
+    16) Month    
+    17) Day    
+    18) Albedo [0-1]  
+    19) Roughness [m]  
+    20) Low  vegetation type [0-20] 
+    21) High vegetation type [0-19] 
+    22) Low  vegetation fractional cover [0-1]
+    23) High vegetation fractional cover [0-1]
+    24) SeaIce cover [0-1]                    
+    25) Snow albedo [0-1]                  
+    26) Snow density [kg/m3]       
+    27) Snow temperature [K]              
+    28) Snow depth [m]                       
+    29) Soil Layer 1 temperature [K]           
+    30) Soil Layer 2 temperature [K]   
+    31) Soil Layer 3 temperature [K]      
+    32) Soil Layer 4 temperature [K]         
+    33) Volumetric Soil Water Layer 1 [m3/m3] 
+    34) Volumetric Soil Water Layer 2 [m3/m3]  
+    35) Volumetric Soil Water Layer 3 [m3/m3]  
+    36) Volumetric Soil Water Layer 4 [m3/m3]  
+    37) Ice temperature Layer 1 [K]       
+    38) Ice temperature Layer 2 [K]       
+    39) Ice temperature Layer 3 [K]           
+    40) Ice temperature Layer 4 [K]
+    
+*_extra.xml.gz [ARTS (gridded field) format] - Contains additional atmospheric variables on different pressure levels:
+    1) Pressure [Pa]
+    2) Fractional cloud cover [0-1]
+    3) Vertical velocity [m/s]
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..ed9efcebd015cd59f7ca4faaf5543c42b351244f
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_extra.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_extra.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..c788b959a41740da8e261c3be3ac2ec7e871371c
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_extra.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_surface.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_surface.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f36eb0a99ba91c7d2174e96724f37bd38175a743
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_ccol_surface.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..57addc0c59bee1ec66eb31eb00d4381cf0720239
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_extra.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_extra.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..b34a1f442ec4d37468eb6f16abd34cb9456b12d6
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_extra.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_surface.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_surface.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..439fda794d33484854d87a6765b48e3309e42198
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_oz_surface.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..ed4bb630add5b7d48bef91b625d53b1179c3205e
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_extra.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_extra.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..438db444bc73530b6194acb2c40726f995536f11
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_extra.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_surface.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_surface.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..4456835cbd4f08187909483b2316110b83ddc738
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_q_surface.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f89ea4af3fd87de4670acfc08048a82b400ed00b
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_extra.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_extra.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f1e416477d63259f898e37c23581b8b79136dfe4
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_extra.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_surface.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_surface.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..fe83bff733db7bc542cff758c7ea3341cea68da1
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_rcol_surface.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..2de11c3b30e0248f71f59ee502b9a7c831710b8d
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_extra.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_extra.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..e96d856990be76326db12402cadb10ceb8b4a739
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_extra.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_surface.xml.gz b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_surface.xml.gz
new file mode 100644
index 0000000000000000000000000000000000000000..c72b77445ac39c5a92770370552c7c17aab1be99
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/eresmaal137_all_t_surface.xml.gz differ
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/extract_arts_l137.f90 b/planets/Earth/ECMWF/IFS/Eresmaa_137L/extract_arts_l137.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6634137d9b104016265f5bdf980d6f207e258a8d
--- /dev/null
+++ b/planets/Earth/ECMWF/IFS/Eresmaa_137L/extract_arts_l137.f90
@@ -0,0 +1,918 @@
+program extract_arts_l137
+
+! This program dumps the Eresmaa/McNally 137-level profile dataset to ARTS-XML format.
+! The source code is combined by using extract_arts_all.f90 (by SAB and DK) from the provided Chevallier datasets and readsaf137.f90 
+! 
+! John Vinzenz Mrziglod 2015-09-01
+
+!   Modified by Daniel Kreyling 12/12/2010
+!
+!   All parameters will be written to one single file per data set.
+!   The files are denoted by "_all_". E.g.: chevallierl91_all_t.xml
+!
+!   Simple program to dump the variables we want from the Chevallier
+!   data set to an ARTS data file. Create a modified copy of this if
+!   you want additional or different variables.
+!
+!   This is the version that extracts also cloud parameters.
+!
+!   Original header from Frederic below. I had to change real to real*8
+!   to make the reading routine work.
+!
+!   SAB 2007-07-19
+!
+!     Description:
+!
+!     Read 137-level diverse profile dataset from ECMWF
+!     Owner:
+!     Frederic Chevallier
+!
+!     History:
+!     Version      Date        Comment
+!     1            22/11/2006  Frederic Chevallier, LSCE
+!     1.1          02/09/2014  Update to 137 levels, Reima Eresmaa, ECMWF 
+!
+!     Code description:
+!       Language:              Fortran 90.
+!       Software Standards:    "European Standards for Writing and Documenting
+!                              Exchangeable Fortran 90 code".
+
+  implicit none
+
+  call readsaf137("t         ", 0)
+  call readsaf137("q             ", 1)
+  call readsaf137("oz            ", 2)
+  call readsaf137("ccol          ", 3)
+  call readsaf137("rcol          ", 4)
+
+end program extract_arts_l137
+
+!--------------------------------------------------------------------
+
+
+SUBROUTINE dewpoint_temp_to_vmrh2o(Td, p, x)
+    ! Convert the dew point temperature and air pressure to volume mixing ratio x.
+    ! FIXME: I am not sure about this formula and its result (JVM)
+    implicit none
+
+    real*8 :: Td                    ! input
+    real*8 :: p                     ! input
+    real*8 :: x                     ! output
+
+    ! Coefficients for Ew:
+    real*8 :: a = -6096.9385
+    real*8 :: b = 21.2409642
+    real*8 :: c = -2.711193e-2
+    real*8 :: d = 1.673952e-5
+    real*8 :: e = 2.433502
+
+    real*8 :: ew
+    ew = exp( a/Td + b + c*Td + d*Td**2 + e*log(Td) )
+
+    ! I worked out this formula after 3.65 of W&H, second edition:
+    x = ew / p
+    
+    return
+end SUBROUTINE dewpoint_temp_to_vmrh2o
+
+SUBROUTINE specific_to_vmrh2o(q,x)
+  ! Convert specific humidity q [kg/kg] to volume mixing ratio x.
+  implicit none
+
+  ! molecular mass of water vapor relative to dry air, from Wallace&Hobbs:
+  real*8, parameter :: eps=0.622 
+
+  real*8 :: q                   ! input
+  real*8 :: x                   ! output
+
+  ! I worked out this formula after 3.59 of W&H, second edition:
+  x = q / ( q + eps*(1-q) )
+  ! It is tempting to divide out q, but then the formula is 
+  ! numerically unstable for small q (division by zero).
+  !
+  ! There is also a simplified formula which could be used 
+  ! approximately for small q values: 
+  !    x = q / ( eps*(1-q) )
+  ! The difference is that in the correct formula the q 
+  ! shows up in the denominator again.
+
+  return
+end SUBROUTINE specific_to_vmrh2o
+
+!--------------------------------------------------------------------
+
+SUBROUTINE specific_to_vmro3(q,x)
+  ! Convert specific ozone q [kg/kg] to volume mixing ratio x.
+  ! This is a clone from the similar subroutine for h2o.
+  implicit none
+
+  ! atomic mass of oxygen: mo3 = 15.9994 
+  !                       (Ulshöfer und Hornschuh Formelsammlung, 3. Auflage)
+  ! effectiv molecular weight of air: mair = 28.97 (Wallace and Hobbs, 2nd edition)
+  ! 
+  ! Molecular weight of O3 relative to (dry) air:
+  ! eps = 3*mo3 / mair
+  real*8, parameter :: eps=1.657
+
+  real*8 :: q                   ! input
+  real*8 :: x                   ! output
+
+  ! I worked out this formula after 3.59 of W&H, second edition:
+  x = q / ( q + eps*(1-q) )
+  ! It is tempting to divide out q, but then the formula is 
+  ! numerically unstable for small q (division by zero).
+  !
+  ! There is also a simplified formula which could be used 
+  ! approximately for small q values: 
+  !    x = q / ( eps*(1-q) )
+  ! The difference is that in the correct formula the q 
+  ! shows up in the denominator again.
+
+  return
+end SUBROUTINE specific_to_vmro3
+
+!--------------------------------------------------------------------
+
+SUBROUTINE condensate_kg_per_kg_to_kg_per_meter3(in, p, T, out)
+  ! Convert condensate in kg/kg to kg/m^3
+  ! in  : condendaste in kg/kg
+  ! p   : pressure in Pa
+  ! T   : temperature in K
+  ! out : condensate in kg/m^3
+  !
+  ! I calculate the volume of 1 kg of air to do the conversion. For
+  ! this, we assume dry air, i.e., I neglect the different mass of
+  ! the H2O that will be also present. I make this simple choice,
+  ! because I am not sure which definition is used inside the ECMWF
+  ! model.  
+  implicit none
+
+  real*8 :: in                  ! input
+  real*8 :: p                   ! input
+  real*8 :: T                   ! input
+  real*8 :: out                 ! output
+
+  ! Gas constant of dry air, according to W&H, second edition:
+  real*8 :: Rd = 287
+
+  ! Calculate volume of 1kg of dry air, according to ideal gas law:
+  real*8 :: V
+  V = Rd*T/p
+
+  ! We have to divide by that to go from kg/kg to kg/m^3:
+  out = in / V
+
+  return
+end SUBROUTINE condensate_kg_per_kg_to_kg_per_meter3
+
+!--------------------------------------------------------------------
+
+SUBROUTINE delta_z(p1,t1,h2o1,p2,t2,h2o2,dz)
+  ! dz : Vertical thickness (in m) between two (p,t,h2o) gridpoints
+  ! p1,p2 : Pressure [Pa]
+  ! t1,t2 : Temperature [K]
+  ! h2o1,h2o2 : Humidity [VMR]
+
+  real*8 :: p1,t1,h2o1,p2,t2,h2o2   ! input
+  real*8 :: dz                      ! output
+
+  ! molecular mass of water vapor relative to dry air, from Wallace&Hobbs:
+  real*8, parameter :: eps=0.622 
+  ! Gravitational constant (in m/s2) taken from Chevalliers f90 routine
+  real*8, parameter :: grav=9.80665  
+  ! Gas constant of dry air in J/(K kg) from Wallace&Hobbs:
+  real*8, parameter :: Rd=287
+
+
+  ! Virtual temperatures [K]:
+  real*8 :: Tv1, Tv2, Tv
+
+  ! Scale height [m]:
+  real*8 :: H1,H2,H
+
+  ! Calculate virtual mean temperature of the layer.
+  ! Virtual temperature from W&H 2nd edition:
+  !  Tv = T / ( 1 - e/p * (1-eps))
+  ! Where e and p are humidity partial pressure and total pressure,
+  ! respectively. (Thus e/p is the H2O VMR.)
+
+  Tv1 = t1 / (1 - h2o1*(1-eps))
+  Tv2 = t2 / (1 - h2o2*(1-eps))
+  Tv = 0.5 * (Tv1+Tv2)
+
+  ! Calculate scale height.
+  ! W&H define it as H = Rd/grav * T
+
+  H1 = Rd/grav * t1
+  H2 = Rd/grav * t2
+  H = 0.5 * (H1+H2)
+
+  dz = H * log(p1/p2)
+
+  return
+end SUBROUTINE delta_z
+
+!--------------------------------------------------------------------
+
+SUBROUTINE wind_speed_and_direction(vx,vy,wspeed,wdir)
+  ! wspeed : wind speed [m/s]
+  ! wdir   : wind direction [0-360 degrees]
+  !           CAUTION: The meteorological wind direction is the direction where the wind comes from.
+  ! vx, vy : 2d vector components of the wind
+
+  implicit none
+
+  real*8 :: vx, vy               ! input
+  real*8 :: wspeed, wdir         ! output
+
+  wspeed = sqrt(vx**2 + vy**2)
+  
+  ! calculate the wind direction
+  if (vx == 0.0) then
+    ! it is a straight wind either from the north or the south
+    if (vy < 0) then
+        wdir = 0    ! wind comes straight from the north
+    else 
+        wdir = 180. ! wind comes straight from the south
+    endif
+  else
+      ! other wind directions
+      wdir = 90 - atan(vy/vx)
+      if ((vx > 0. .and. vy > 0.) .or. (vx > 0 .and. vy < 0) ) then
+        wdir = wdir + 180.
+      endif
+  endif
+  
+
+  return
+end SUBROUTINE wind_speed_and_direction
+
+!--------------------------------------------------------------------
+
+  SUBROUTINE readsaf137(cvar, findex)
+
+    implicit none
+
+    integer, parameter      :: nlev = 137 ! number of vertical levels
+    
+    ! Gravitational constant (in m/s2) taken from Chevalliers f90 routine
+    real*8, parameter :: grav=9.80665  
+
+    real*8, dimension(nlev)   :: temp, hum, ozo, cc, clw, ciw
+    real*8, dimension(nlev)   :: rain, snow, w, pap
+    real*8, dimension(nlev+1) :: paph
+    real*8                    :: lnpsurf, psurf, z0, tsurf, t2m, td2m, hum2m, u10, v10, wspeed, wdirection
+    real*8                    :: stratrsrf, convrsrf, snowsurf, elev
+    real*8                    :: lsm, lat, long
+    real*8                    :: alb, sr, tvl, tvh, cvl, cvh, seaice
+    real*8                    :: asn, rsn, tsn, sd, stl1, stl2, stl3, stl4
+    real*8                    :: swvl1, swvl2, swvl3, swvl4, istl1, istl2, istl3, istl4
+    integer                 :: year, month, day, step, gpoint, ind, i, j, findex
+    character(len=10)       :: cvar
+    
+    ! VMR parameters of H2O and O3:
+    real*8                    :: vmrh2o ! H2O Volume mixing ratio
+    real*8                    :: vmro3  ! O3 Volume mixing ratio
+
+    ! Condensate parameters per volume:
+    real*8                    :: vol_clw    
+    real*8                    :: vol_ciw    
+
+    ! For vertical thickness calculation:
+    real*8 :: p1,t1,h2o1,p2,t2,h2o2,dz,zprof
+
+    if (trim(cvar) .ne. 't'    .and. trim(cvar) .ne. 'q'   .and. trim(cvar) .ne. 'oz' .and. &
+        trim(cvar) .ne. 'ccol' .and. trim(cvar) .ne. 'rcol') then
+      write(*,*) 'Wrong variable, start again'
+      stop
+    endif
+
+    open(1,file='profiles137/nwp_saf_'//trim(cvar)//'_sampled.atm',form='formatted')
+    open(2,file='profiles137/nwp_saf_'//trim(cvar)//'_sampled.sfc',form='formatted')
+    
+    ! ARTS-1-1 output file atmosphere data
+    open(3,file='eresmaal137_all_'//trim(cvar)//'.xml.temp',form='formatted')
+    write(3,'(A)') '<?xml version="1.0"?>'
+    write(3,'(A)') '<arts format="ascii" version="1">                        '
+    write(3,'(A)') '<comment>                                                '
+    write(3,'(A)') 'This file contains one matrix for each atmospheric profile.'
+    write(3,'(A)') 'Each matrix row corresponds to one pressure level. The   '
+    write(3,'(A)') 'meaning of the columns is:                               '
+!    write(3,'(A)') 'p[Pa] T[K] z[m] H2O_VMR[1] O3_VMR[1] CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)]'
+    write(3,'(A)') 'p[Pa] T[K] z[m] CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)] H2O_VMR[1] O3_VMR[1]'
+!    write(3,'(A)') 'CLW[kg/m^3] CIW[kg/m^3] Rain[kg/(m2*s)] Snow[kg/(m2*s)]'
+    write(3,'(A)') '</comment>                                               '
+    write(3,'(A)') '<Array type="Matrix" nelem="5000">                       '
+    
+    ! ARTS-1-1 output file surface variables and meta data
+    open(4,file='eresmaal137_all_'//trim(cvar)//'_surface.xml',form='formatted')
+    write(4,'(A)') '<?xml version="1.0"?>'
+    write(4,'(A)') '<arts format="ascii" version="1">                           '
+    write(4,'(A)') '<comment>                                                   '
+    write(4,'(A)') 'This file contains one row for each atmospheric profile. '
+    write(4,'(A)') 'The meaning of the columns is:'
+    write(4,'(A)') '1) profile index [1-25000]'//NEW_LINE('A')//'&
+                    &2) Ln of Surf Pressure in [Pa]'//NEW_LINE('A')//'&
+                    &3) Surface geopotential [m2/s2] '//NEW_LINE('A')//'&
+                    &4) Surface Skin Temperature [K]'//NEW_LINE('A')//'&
+                    &5) 2m Temperature [K]'//NEW_LINE('A')//'&
+                    &6) 2m Dew point temperature [K]'//NEW_LINE('A')//'&
+                    &7) 10m wind speed [m/s]'//NEW_LINE('A')//'&
+                    &8) 10m wind direction [0-360 DEG]'//NEW_LINE('A')//'&
+                    &9) Stratiform precipitation at surface [m]'//NEW_LINE('A')//'&
+                    &10) Convective precipitation at surface [m]'//NEW_LINE('A')//'&
+                    &11) Snowfall at surface [m] (water equival.)'//NEW_LINE('A')//'&
+                    &12) Land/sea Mask [0-1]'//NEW_LINE('A')//'&
+                    &13) Latitude [deg]'//NEW_LINE('A')//'&
+                    &14) Longitude [deg]'//NEW_LINE('A')//'&
+                    &15) Year'//NEW_LINE('A')//'&
+                    &16) Month    '//NEW_LINE('A')//'&
+                    &17) Day    '//NEW_LINE('A')//'&
+                    &18) Albedo [0-1]  '//NEW_LINE('A')//'&
+                    &19) Roughness [m]  '//NEW_LINE('A')//'&
+                    &20) Low  vegetation type [0-20] '//NEW_LINE('A')//'&
+                    &21) High vegetation type [0-19] '//NEW_LINE('A')//'&
+                    &22) Low  vegetation fractional cover [0-1]'//NEW_LINE('A')//'&
+                    &23) High vegetation fractional cover [0-1]'//NEW_LINE('A')//'&
+                    &24) SeaIce cover [0-1]                    '//NEW_LINE('A')//'&
+                    &25) Snow albedo [0-1]                  '//NEW_LINE('A')//'&
+                    &26) Snow density [kg/m3]       '//NEW_LINE('A')//'&
+                    &27) Snow temperature [K]              '//NEW_LINE('A')//'&
+                    &28) Snow depth [m]                       '//NEW_LINE('A')//'&
+                    &29) Soil Layer 1 temperature [K]           '//NEW_LINE('A')//'&
+                    &30) Soil Layer 2 temperature [K]   '//NEW_LINE('A')//'&
+                    &31) Soil Layer 3 temperature [K]      '//NEW_LINE('A')//'&
+                    &32) Soil Layer 4 temperature [K]         '//NEW_LINE('A')//'&
+                    &33) Volumetric Soil Water Layer 1 [m3/m3] '//NEW_LINE('A')//'&
+                    &34) Volumetric Soil Water Layer 2 [m3/m3]  '//NEW_LINE('A')//'&
+                    &35) Volumetric Soil Water Layer 3 [m3/m3]  '//NEW_LINE('A')//'&
+                    &36) Volumetric Soil Water Layer 4 [m3/m3]  '//NEW_LINE('A')//'&
+                    &37) Ice temperature Layer 1 [K]       '//NEW_LINE('A')//'&
+                    &38) Ice temperature Layer 2 [K]       '//NEW_LINE('A')//'&
+                    &39) Ice temperature Layer 3 [K]           '//NEW_LINE('A')//'&
+                    &40) Ice temperature Layer 4 [K]'
+
+    write(4,'(A)') '</comment>                                                  '
+    !write(4,'(A)') '<Vector nelem="8">'
+    !write(4,'(A)') 'pindex[1-25000] psurf[Pa] Tsurf[K] wspeed[m/s] wdir[0-360°] land-sea[0-1] sea-ice[0-1] snow-depth[m]'
+    !write(4,'(A)') '</Vector>'
+    write(4,'(A)') '<Matrix nrows="5000" ncols="40">'
+
+    ! ARTS-1-1 output file fractional cloud cover
+    open(5,file='eresmaal137_all_'//trim(cvar)//'_extra.xml.temp',form='formatted')
+    write(5,'(A)') '<?xml version="1.0"?>'
+    write(5,'(A)') '<arts format="ascii" version="1">                           '
+    write(5,'(A)') '<comment>                                                   '
+    write(5,'(A)') 'This file contains one matrix for each atmospheric profile. '
+    write(5,'(A)') 'The meaning of column is:'
+    write(5,'(A)') 'p[Pa] cc[0-1] w[m/s]'
+    write(5,'(A)') '</comment>                                                  '
+    write(5,'(A)') '<Array type="Matrix" nelem="5000">                          '
+
+    i = 0
+    do
+      ! read the atmospheric profile
+      read(1,'(1246(e13.6,x),i5,3i3,i7,i5)',end=1) &
+         temp(:),     &! 1) Temperature [K]                          (1-137)
+         hum(:),      &! 2) Humidity [kg/kg]                         (138-274)
+         ozo(:),      &! 3) Ozone [kg/kg]                            (275-411)
+         cc(:),       &! 4) Cloud Cover [0-1]                        (412-548)
+         clw(:),      &! 5) C Liquid W [kg/kg]                       (549-685)
+         ciw(:),      &! 6) C Ice W [kg/kg]                          (686-822)
+         rain(:),     &! 7) Rain [kg/(m2 *s)]                        (823-959)
+         snow(:),     &! 8) Snow [kg/(m2 *s)]                        (960-1096)
+         w(:),        &! 9) Vertical Velocity [Pa/s]                 (1097-1233)
+         lnpsurf,     &!10) Ln of Surf Pressure in [Pa]              (1234)
+         z0,          &!11) Surface geopotential [m2/s2]             (1235) 
+         tsurf,       &!12) Surface Skin Temperature [K]             (1236)
+         t2m,         &!13) 2m Temperature [K]                       (1237)
+         td2m,        &!14) 2m Dew point temperature [K]             (1238)
+         u10,         &!15) 10m wind speed U component [m/s]         (1239)
+         v10,         &!16) 10m wind speed V component [m/s]         (1240)
+         stratrsrf,   &!17) Stratiform precipitation at surface [m]  (1241)
+         convrsrf,    &!18) Convective precipitation at surface [m]  (1242)
+         snowsurf,    &!19) Snowfall at surface [m] (water equival.) (1243)
+         lsm,         &!20) Land/sea Mask [0-1]                      (1244)
+         lat,         &!21) Latitude [deg]                           (1245)
+         long,        &!22) Longitude [deg]                          (1246)
+         year,        &!23) Year                                     (1247)
+         month,       &!24) Month                                    (1248)
+         day,         &!25) Day                                      (1249)
+         step,        &!26) Step                                     (1250)
+         gpoint,      &!27) Grid point [1-2140702]                   (1251)
+         ind           !28) Index (rank-sorted)                      (1252) 
+
+      read(2,'(26(e13.6,x),i5,3i3,i8,i5)',end=1) &
+         alb,         &! 1) Albedo [0-1]                             (1)
+         sr,          &! 2) Roughness [m]                            (2)
+         tvl,         &! 3) Low  vegetation type [0-20]              (3)
+         tvh,         &! 4) High vegetation type [0-19]              (4)
+         cvl,         &! 5) Low  vegetation fractional cover [0-1]   (5)
+         cvh,         &! 6) High vegetation fractional cover [0-1]   (6)
+         seaice,      &! 7) SeaIce cover [0-1]                       (7)
+         asn,         &! 8) Snow albedo [0-1]                        (8)
+         rsn,         &! 9) Snow density [kg/m3]                     (9)
+         tsn,         &!10) Snow temperature [K]                    (10)
+         sd,          &!11) Snow depth [m]                          (11)
+         stl1,        &!12) Soil Layer 1 temperature [K]            (12)
+         stl2,        &!13) Soil Layer 2 temperature [K]            (13)
+         stl3,        &!14) Soil Layer 3 temperature [K]            (14)
+         stl4,        &!15) Soil Layer 4 temperature [K]            (15)
+         swvl1,       &!16) Volumetric Soil Water Layer 1 [m3/m3]   (16)
+         swvl2,       &!17) Volumetric Soil Water Layer 2 [m3/m3]   (17)
+         swvl3,       &!18) Volumetric Soil Water Layer 3 [m3/m3]   (18)
+         swvl4,       &!19) Volumetric Soil Water Layer 4 [m3/m3]   (19)
+         istl1,       &!20) Ice temperature Layer 1 [K]             (20)
+         istl2,       &!21) Ice temperature Layer 2 [K]             (21)
+         istl3,       &!22) Ice temperature Layer 3 [K]             (22)
+         istl4,       &!23) Ice temperature Layer 4 [K]             (23)
+    ! The following variables also appear in .atm file
+         lsm,         &!%24) Land/sea Mask [0-1]                    (24)
+         lat,         &!%25) Latitude [deg]                         (25)
+         long,        &!%26) Longitude[deg]                         (26)
+         year,        &!%27) Year of start of forecast              (27)
+         month,       &!%28) Month of start of forecast             (28)
+         day,         &!%29) Day of start of forecast               (29)
+         step,        &!%30) Forecast step [hours]                  (30)
+         gpoint,      &!%31) Grid point [1-2140702]                 (31)
+         ind           !%32) Index (rank-sorted)                    (32) 
+
+      ! compute surface pressure [Pa]
+      psurf = exp(lnpsurf)
+      
+      ! compute wind speed and direction
+      call wind_speed_and_direction(u10, v10, wspeed, wdirection)
+      
+      ! Write to ARTS output file (surface variables and meta data):
+      write(4,'(I5,13(e13.6,x),i5,2i3,23(e13.6,x))') &
+        findex*5000+i+1,    &!1) file index
+        lnpsurf,            &!2) Ln of Surf Pressure in [Pa]
+        z0,                 &!3) Surface geopotential [m2/s2] 
+        tsurf,              &!4) Surface Skin Temperature [K]
+        t2m,                &!5) 2m Temperature [K]
+        td2m,               &!6) 2m Dew point temperature [K]
+        wspeed,             &!7) 10m wind speed [m/s]
+        wdirection,         &!8) 10m wind direction [0-360 DEG]
+        stratrsrf,          &!9) Stratiform precipitation at surface [m]
+        convrsrf,           &!10) Convective precipitation at surface [m]
+        snowsurf,           &!11) Snowfall at surface [m] (water equival.)
+        lsm,                &!12) Land/sea Mask [0-1]
+        lat,                &!13) Latitude [deg]
+        long,               &!14) Longitude [deg]
+        year,               &!15) Year
+        month,              &!16) Month    
+        day,                &!17) Day    
+        alb,                &!18) Albedo [0-1]  
+        sr,                 &!19) Roughness [m]  
+        tvl,                &!20) Low  vegetation type [0-20] 
+        tvh,                &!21) High vegetation type [0-19] 
+        cvl,                &!22) Low  vegetation fractional cover [0-1]
+        cvh,                &!23) High vegetation fractional cover [0-1]
+        seaice,             &!24) SeaIce cover [0-1]                    
+        asn,                &!25) Snow albedo [0-1]                  
+        rsn,                &!26) Snow density [kg/m3]       
+        tsn,                &!27) Snow temperature [K]              
+        sd,                 &!28) Snow depth [m]                       
+        stl1,               &!29) Soil Layer 1 temperature [K]           
+        stl2,               &!30) Soil Layer 2 temperature [K]   
+        stl3,               &!31) Soil Layer 3 temperature [K]      
+        stl4,               &!32) Soil Layer 4 temperature [K]         
+        swvl1,              &!33) Volumetric Soil Water Layer 1 [m3/m3] 
+        swvl2,              &!34) Volumetric Soil Water Layer 2 [m3/m3]  
+        swvl3,              &!35) Volumetric Soil Water Layer 3 [m3/m3]  
+        swvl4,              &!36) Volumetric Soil Water Layer 4 [m3/m3]  
+        istl1,              &!37) Ice temperature Layer 1 [K]       
+        istl2,              &!38) Ice temperature Layer 2 [K]       
+        istl3,              &!39) Ice temperature Layer 3 [K]           
+        istl4              !40) Ice temperature Layer 4 [K]     
+
+!      write(4,101) findex*5000+i+1, lnpsurf, tsurf, wspeed, wdirection, lsm, seaice, sd, lat, long
+!101   format(I5, ES11.4, F7.2, F7.2, F7.2, 2F10.6, F7.2, 2F10.3)
+
+      ! compute surface geometric height [m]
+      elev = z0 / 9.80665
+
+      ! compute full-level pressure (pap) and half-level pressures (paph) [Pa]
+      ! note that all profile variables above are given at full levels
+      call ec_p137l(psurf,pap,paph)
+      
+      ! Write to ARTS output file (atmospheric profile):
+
+       ! We first output the quantities for the surface. Profile
+       ! quantities given by Chevallier are layer means, not point
+       ! values!
+       ! So, we add another point at the bottom, then simply
+       ! interpret the profile values as if they were point values.
+       ! For the surface, we take the 2m temperature and humidity. We
+       ! do not have a 2m O3 value, so we simply copy the value from
+       ! the lowest layer. 
+
+       ! Calculate vmrh2o from the dew point temperature and air pressure
+       call dewpoint_temp_to_vmrh2o(td2m,psurf,vmrh2o)
+
+       ! Calculate vmro3 from specific ozone:
+       call specific_to_vmro3(ozo(nlev),vmro3)
+
+       ! The 4 condensate parameters are also copied from the lowest
+       ! layer. However, we use psurf and t2m for the conversion from
+       ! kg/kg to kg/m^3.
+
+       ! Calculate cloud liquid water content:
+       call condensate_kg_per_kg_to_kg_per_meter3(clw(nlev), psurf, t2m, vol_clw)
+
+       ! Calculate cloud ice water content:
+       call condensate_kg_per_kg_to_kg_per_meter3(ciw(nlev), psurf, t2m, vol_ciw)
+
+       write(3,'(A)') '<Matrix nrows="138" ncols="9">'
+!       write(3,100) psurf, t2m, elev, vmrh2o, vmro3, &
+!            & vol_clw, vol_ciw, rain(nlev), snow(nlev)
+!100    format(ES11.4, F7.2, F9.2, 2ES9.2, 2ES9.2, 2ES10.2)
+       write(3,100) psurf, t2m, elev, vol_clw, vol_ciw, rain(nlev), snow(nlev), vmrh2o, vmro3
+100    format(ES11.4, F7.2, F9.2, 2ES9.2, 2ES10.2, 2ES9.2)
+
+       !  Write to ARTS output file (fractional cloud cover):
+       write(5,'(A)') '<Matrix nrows="137" ncols="3">'
+
+       ! Initialize quantities for hydrostatic z profile:
+       zprof = elev
+       p1 = psurf
+       t1 = t2m
+       h2o1 = vmrh2o
+
+
+       do j=nlev,1,-1
+          ! Calculate vmrh2o from specific humidity:
+          call specific_to_vmrh2o(hum(j),vmrh2o)
+       
+          ! Calculate vmro3 from specific ozone:
+          call specific_to_vmro3(ozo(j),vmro3)
+
+          ! Quantities for hydrostatic z profile:
+          p2 = pap(j)
+          t2 = temp(j)
+          h2o2 = vmrh2o
+
+          ! Calculate thickness between this point and the last:
+          call delta_z(p1,t1,h2o1,p2,t2,h2o2,dz)
+          
+          zprof = zprof + dz
+
+          ! Calculate cloud liquid water content:
+          call condensate_kg_per_kg_to_kg_per_meter3(clw(j), pap(j), temp(j), vol_clw)
+
+          ! Calculate cloud ice water content:
+          call condensate_kg_per_kg_to_kg_per_meter3(ciw(j), pap(j), temp(j), vol_ciw)
+
+!          write(3,100) pap(j), temp(j), zprof, vmrh2o, vmro3, &
+!            & vol_clw, vol_ciw, rain(j), snow(j)
+          write(3,100) pap(j), temp(j), zprof, vol_clw, vol_ciw, rain(j), snow(j), vmrh2o, vmro3
+!          write(3,100) vol_clw, vol_ciw, rain(j), snow(j)
+
+          ! output fractional cloud cover
+          write(5,'(ES11.4, F10.6, F10.6)') pap(j), cc(j), w(j)
+
+          ! Update hydrostatic quantities for next iteration
+          p1 = p2
+          t1 = t2
+          h2o1 = h2o2
+
+       end do
+
+       write(3,'(A)') '</Matrix>'
+       write(5,'(A)') '</Matrix>'
+      
+      i = i + 1
+    enddo
+    1 write(*,*) 'Number of profiles found in the files:',i
+    ! Close files:
+    close(1)
+    close(2)
+
+    write(3,'(A)') '</Array>'
+    write(3,'(A)') '</arts> '
+    close(3)
+    
+    write(4,'(A)') '</Matrix>'
+    write(4,'(A)') '</arts> '
+    close(4)
+    
+    write(5,'(A)') '</Array>'
+    write(5,'(A)') '</arts> '
+    close(5)
+  end SUBROUTINE readsaf137
+
+! ----------------------------------------
+!
+  SUBROUTINE EC_P137l(spres,pap,paph)
+!
+!    This software was developed within the context of
+!    the EUMETSAT Satellite Application Facility on
+!    Numerical Weather Prediction (NWP SAF), under the
+!    Cooperation Agreement dated 25 November 1998, between
+!    EUMETSAT and the Met Office, UK, by one or more partners
+!    within the NWP SAF. The partners in the NWP SAF are
+!    the Met Office, ECMWF, KNMI and MeteoFrance.
+!
+!    Copyright 2014, EUMETSAT, All Rights Reserved.
+!
+!     Description:
+!     Computes the 137-level vertical pressure grid
+!       associated to the input surface pressure
+!     All pressures are in Pa
+
+      implicit none
+      integer, parameter    :: nlev=137
+      integer               :: jk
+      real*8      :: spres
+      real*8      :: aam(nlev+1), bbm(nlev+1)
+      real*8      :: pap(nlev), paph(nlev+1)
+
+      data aam / &
+         0.000000, &
+         2.000365, &
+         3.102241, &
+         4.666084, &
+         6.827977, &
+         9.746966, &
+        13.605424, &
+        18.608931, &
+        24.985718, &
+        32.985710, &
+        42.879242, &
+        54.955463, &
+        69.520576, &
+        86.895882, &
+       107.415741, &
+       131.425507, &
+       159.279404, &
+       191.338562, &
+       227.968948, &
+       269.539581, &
+       316.420746, &
+       368.982361, &
+       427.592499, &
+       492.616028, &
+       564.413452, &
+       643.339905, &
+       729.744141, &
+       823.967834, &
+       926.344910, &
+      1037.201172, &
+      1156.853638, &
+      1285.610352, &
+      1423.770142, &
+      1571.622925, &
+      1729.448975, &
+      1897.519287, &
+      2076.095947, &
+      2265.431641, &
+      2465.770508, &
+      2677.348145, &
+      2900.391357, &
+      3135.119385, &
+      3381.743652, &
+      3640.468262, &
+      3911.490479, &
+      4194.930664, &
+      4490.817383, &
+      4799.149414, &
+      5119.895020, &
+      5452.990723, &
+      5798.344727, &
+      6156.074219, &
+      6526.946777, &
+      6911.870605, &
+      7311.869141, &
+      7727.412109, &
+      8159.354004, &
+      8608.525391, &
+      9076.400391, &
+      9562.682617, &
+     10065.978516, &
+     10584.631836, &
+     11116.662109, &
+     11660.067383, &
+     12211.547852, &
+     12766.873047, &
+     13324.668945, &
+     13881.331055, &
+     14432.139648, &
+     14975.615234, &
+     15508.256836, &
+     16026.115234, &
+     16527.322266, &
+     17008.789062, &
+     17467.613281, &
+     17901.621094, &
+     18308.433594, &
+     18685.718750, &
+     19031.289062, &
+     19343.511719, &
+     19620.042969, &
+     19859.390625, &
+     20059.931641, &
+     20219.664062, &
+     20337.863281, &
+     20412.308594, &
+     20442.078125, &
+     20425.718750, &
+     20361.816406, &
+     20249.511719, &
+     20087.085938, &
+     19874.025391, &
+     19608.572266, &
+     19290.226562, &
+     18917.460938, &
+     18489.707031, &
+     18006.925781, &
+     17471.839844, &
+     16888.687500, &
+     16262.046875, &
+     15596.695312, &
+     14898.453125, &
+     14173.324219, &
+     13427.769531, &
+     12668.257812, &
+     11901.339844, &
+     11133.304688, &
+     10370.175781, &
+      9617.515625, &
+      8880.453125, &
+      8163.375000, &
+      7470.343750, &
+      6804.421875, &
+      6168.531250, &
+      5564.382812, &
+      4993.796875, &
+      4457.375000, &
+      3955.960938, &
+      3489.234375, &
+      3057.265625, &
+      2659.140625, &
+      2294.242188, &
+      1961.500000, &
+      1659.476562, &
+      1387.546875, &
+      1143.250000, &
+       926.507812, &
+       734.992188, &
+       568.062500, &
+       424.414062, &
+       302.476562, &
+       202.484375, &
+       122.101562, &
+        62.781250, &
+        22.835938, &
+         3.757813, &
+         0.000000, &
+         0.000000 /
+
+      data bbm / &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000000, &
+     0.000007, &
+     0.000024, &
+     0.000059, &
+     0.000112, &
+     0.000199, &
+     0.000340, &
+     0.000562, &
+     0.000890, &
+     0.001353, &
+     0.001992, &
+     0.002857, &
+     0.003971, &
+     0.005378, &
+     0.007133, &
+     0.009261, &
+     0.011806, &
+     0.014816, &
+     0.018318, &
+     0.022355, &
+     0.026964, &
+     0.032176, &
+     0.038026, &
+     0.044548, &
+     0.051773, &
+     0.059728, &
+     0.068448, &
+     0.077958, &
+     0.088286, &
+     0.099462, &
+     0.111505, &
+     0.124448, &
+     0.138313, &
+     0.153125, &
+     0.168910, &
+     0.185689, &
+     0.203491, &
+     0.222333, &
+     0.242244, &
+     0.263242, &
+     0.285354, &
+     0.308598, &
+     0.332939, &
+     0.358254, &
+     0.384363, &
+     0.411125, &
+     0.438391, &
+     0.466003, &
+     0.493800, &
+     0.521619, &
+     0.549301, &
+     0.576692, &
+     0.603648, &
+     0.630036, &
+     0.655736, &
+     0.680643, &
+     0.704669, &
+     0.727739, &
+     0.749797, &
+     0.770798, &
+     0.790717, &
+     0.809536, &
+     0.827256, &
+     0.843881, &
+     0.859432, &
+     0.873929, &
+     0.887408, &
+     0.899900, &
+     0.911448, &
+     0.922096, &
+     0.931881, &
+     0.940860, &
+     0.949064, &
+     0.956550, &
+     0.963352, &
+     0.969513, &
+     0.975078, &
+     0.980072, &
+     0.984542, &
+     0.988500, &
+     0.991984, &
+     0.995003, &
+     0.997630, &
+     1.000000 /
+
+      do jk=1,nlev+1
+        paph(jk)=aam(jk)+bbm(jk)*spres
+      end do
+      do jk=1,nlev
+        pap(jk)=0.5*(paph(jk)+paph(jk+1))
+      end do
+
+      return
+  end subroutine ec_p137l
diff --git a/planets/Earth/ECMWF/IFS/Eresmaa_137L/nwp_saf_profiles137.pdf b/planets/Earth/ECMWF/IFS/Eresmaa_137L/nwp_saf_profiles137.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..2c8a1a7064a16f28035335ce160abaedb6ccf287
Binary files /dev/null and b/planets/Earth/ECMWF/IFS/Eresmaa_137L/nwp_saf_profiles137.pdf differ