Index: /palm/trunk/SCRIPTS/.cvd.config.default
===================================================================
--- /palm/trunk/SCRIPTS/.cvd.config.default	(revision 4400)
+++ /palm/trunk/SCRIPTS/.cvd.config.default	(revision 4400)
@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM 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 3 of the License, or (at your option) any later version.
+#
+# PALM 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
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2020  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revision
+# 
+# 
+# 3668 2019-01-14 12:49:24Z maronga
+#
+# Description:
+# ------------
+# Processing tool for creating PIDS conform virtual measurement setup file
+# from UC2 data-standard conform observational data or from prescribed input
+# coordinates.
+#
+# @Author Matthias Suehring (suehring@muk.uni-hannover.de)
+#
+#--------------------------------------------------------------------------------#
+
+# Section where some global file attributes can be specified.
+[global]
+
+author = Suehring, Matthias, suehring@muk.uni-hannover.de
+contact_person = Suehring, Matthias, suehring@muk.uni-hannover.de
+acronym = LUHimuk
+comment = created with palm_cvd
+data_content = measurement coordinates
+dependencies = none
+keywords = virtual measurement setup file
+source = none
+campaign = winter iop 1
+location = Berlin
+site = Ernst Reuter Platz
+institution = Institute of Meteorology and Climatology, Leibniz University Hannover
+palm_version = 6.0
+references = none
+
+# Section for observational input data. If non observational input data is available
+# just omit this section.
+[input]
+
+data_path = < path to your observational input data >
+
+# Section to specify the output path and the output filename.
+[output]
+
+output_path     = < your output path >
+output_filename = vm_driver
+
+# Section to specify further sampling locations and sampling variables. Please note,
+# at the moment only timeseries can be specified. 
+# If non customized sampling location is desired, just set number_positions to zero.
+[custom_positions]
+
+# Give the number of custom sampling coordinates
+number_positions = 4
+
+# Give site coordinates: E_UTM (m), N_UTM (m), z (m).
+# Please give z in meters above ground.
+#
+#             [ E_UTM,    N_UTM,      z (agl) ]
+
+coordinates1: [ 387670.0, 5819090.0,  4.0     ]
+coordinates2: [ 387676.0, 5819139.0,  1.0     ]
+coordinates3: [ 387630.0, 5819065.0,  41.0    ]
+coordinates4: [ 387270.0, 5819191.0,  8.5     ]
+# ...
+
+# Give the variables that shall be sampled. Note, there is a 
+# number of default variables that will be always sampled, which is:
+# u, v, w, theta and hus (mixing ratio).
+# On top of these, you can individually prescribe further variables
+# for each site. In order to prescribe further variables, please
+# follow the naming convention defined in the [UC]2 data standard.
+
+vars_to_be_measured1: [ "u", "v", "w", "theta", "rlus", "rsus" ]
+vars_to_be_measured2: [ "rnds", "wdir" ]
+vars_to_be_measured3: [ "us" ]
+vars_to_be_measured4: [ "utheta", "vtheta", "wtheta", "wspeed" ]
+# ...
Index: /palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- /palm/trunk/SCRIPTS/.palm.iofiles	(revision 4399)
+++ /palm/trunk/SCRIPTS/.palm.iofiles	(revision 4400)
@@ -63,6 +63,4 @@
 SURFACE_DATA_AV_NETCDF*    out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   nc
 #
-VIRTUAL_MEAS_BIN*          out:lnpe *         $base_data/$run_identifier/OUTPUT         _vmeas
-#
 DATA_1D_FL_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _fl        nc
 DATA_1D_PR_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _pr        nc
@@ -81,4 +79,5 @@
 DATA_MASK_AV_NETCDF*       out:tr   *         $base_data/$run_identifier/OUTPUT         _av_masked nc
 DATA_AGT_NETCDF*           out:tr   *         $base_data/$run_identifier/OUTPUT         _agt       nc
+VM_OUTPUT*                 out:pe   *         $base_data/$run_identifier/OUTPUT         _vm
 #
 DATA_PRT_NETCDF*           out:pe   *         $base_data/$run_identifier/OUTPUT         _prt
Index: /palm/trunk/SCRIPTS/palm_cvd
===================================================================
--- /palm/trunk/SCRIPTS/palm_cvd	(revision 4400)
+++ /palm/trunk/SCRIPTS/palm_cvd	(revision 4400)
@@ -0,0 +1,551 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+#--------------------------------------------------------------------------------#
+# This file is part of the PALM model system.
+#
+# PALM 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 3 of the License, or (at your option) any later version.
+#
+# PALM 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
+# PALM. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright 1997-2020  Leibniz Universitaet Hannover
+#--------------------------------------------------------------------------------#
+#
+# Current revisions:
+# -----------------
+# 
+# 
+# Former revisions:
+# -----------------
+# $Id$
+# Initial revision
+# 
+# 
+# $
+#
+# Description:
+# ------------
+# Processing tool for creating PIDS conform virtual measurement setup file
+# from UC2 data-standard conform observational data or from prescribed input
+# coordinates.
+#
+# @Authors Matthias Suehring (suehring@muk.uni-hannover.de)
+#          Tobias Gronemeier (gronemeier@muk.uni-hannover.de)
+#
+# @todo Add further feature tpyes for customized observations. At the moment only
+#       timeSeries is possible.
+#--------------------------------------------------------------------------------#
+
+
+import netCDF4
+from netCDF4 import Dataset, stringtochar
+import os
+import numpy as np
+
+
+# Function to read the config file
+def read_config_file():
+
+   import configparser
+   import os
+   import sys
+   import json
+
+   # Definition of global configuration parameters
+   global global_acronym     
+   global global_author      
+   global global_campaign    
+   global global_comment     
+   global global_contact     
+   global global_data_content
+   global global_dependencies
+   global global_institution 
+   global global_keywords    
+   global global_location    
+   global global_references  
+   global global_site        
+   global global_source      
+   global global_palm_version
+   global data_path          
+   global output_path        
+   global output_filename    
+   global number_positions
+   global input_from_observations
+   global coordinates
+   global vars_to_be_measured
+   global custom_coordinates
+
+   global_acronym          = " "
+   global_author           = " "
+   global_campaign         = " "
+   global_comment          = " "
+   global_contact          = " "
+   global_data_content     = " "
+   global_dependencies     = " "
+   global_institution      = " "
+   global_keywords         = " "
+   global_location         = " "
+   global_references       = " "
+   global_site             = " "
+   global_source           = " "
+   global_palm_version     = 6.0
+   data_path               = " "
+   output_path             = " "
+   output_filename         = "none"
+   number_positions        = -999
+   input_from_observations = False
+   coordinates             = []
+   vars_to_be_measured     = []
+   custom_coordinates      = False
+
+   # Check if configuration files exists and quit otherwise
+   input_config = ".cvd.config.default"
+   for i in range(1,len(sys.argv)): 
+      input_config = str(sys.argv[i])
+
+   # Allow empty settings
+   config = configparser.RawConfigParser(allow_no_value=True)
+
+   # Check if a config file exists.
+   if ( os.path.isfile(input_config) == False ):
+      print ("Error. No configuration file " + input_config + " found.")
+      quit()
+
+   config.read(input_config)
+   
+   for section in range( 0, len( config.sections() ) ):
+
+      current_section = config.sections()[section]
+
+      # read global attributes which are written into the output file header
+      if ( current_section == 'global' ):
+
+         global_acronym      = config.get( current_section, 'acronym'        )
+         global_author       = config.get( current_section, 'author'         )
+         global_campaign     = config.get( current_section, 'campaign'       )
+         global_comment      = config.get( current_section, 'comment'        )
+         global_contact      = config.get( current_section, 'contact_person' )
+         global_data_content = config.get( current_section, 'data_content'   )
+         global_dependencies = config.get( current_section, 'dependencies'   )
+         global_institution  = config.get( current_section, 'institution'    )
+         global_keywords     = config.get( current_section, 'keywords'       )
+         global_location     = config.get( current_section, 'location'       )
+         global_references   = config.get( current_section, 'references'     )
+         global_site         = config.get( current_section, 'site'           )
+         global_source       = config.get( current_section, 'source'         )
+         global_palm_version = float( config.get( current_section, 'palm_version' ) )
+
+      # Read data input path for observational data
+      elif ( current_section == 'input' ):
+
+         data_path = config.get( current_section, 'data_path' )
+         input_from_observations = True
+
+      # Read output path and filename for the VM driver
+      elif ( current_section == 'output' ):
+
+         output_path     = config.get( current_section, 'output_path' )
+         output_filename = config.get( current_section, 'output_filename' )
+
+      # Read customized coordinates where virtual measurements shall be taken, 
+      # as well as the variables that should be sampled.
+      elif ( current_section == 'custom_positions' ):
+
+         number_positions = config.get( current_section, 'number_positions' )
+
+         for count in range( 0, int( number_positions ) ):
+            coordinates.append( json.loads( config.get( current_section, \
+                                                        "coordinates" + str( count + 1 ) ) ) )
+            # If coordinates are given, set a global flag.
+            custom_coordinates = True
+
+         for count in range( 0, int( number_positions ) ):
+            vars_to_be_measured.append( json.loads( config.get( current_section, \
+                                                    "vars_to_be_measured" + str( count + 1 ) ) ) )
+
+
+   return 0
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Main program:
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+# Define strings
+name_featuretype   = "featureType"
+name_ts            = "timeSeries"
+name_traj          = "trajectory"
+name_ntime         = "ntime"
+name_time          = "time"
+name_station       = "station"
+name_traj_dim      = "traj"
+name_nz            = "nz"
+name_datacontent   = "data_content"
+name_eutm          = "E_UTM"
+name_nutm          = "N_UTM"
+name_hao           = "height_above_origin"
+name_station_h     = "station_h"
+name_z             = "z"
+name_soil_sampling = "soil_sample"
+name_num_stat      = "number_of_stations"
+name_fill          = "_FillValue"
+name_site          = "site"
+name_orig_x        = "origin_x"
+name_orig_y        = "origin_y"
+name_orig_z        = "origin_z"
+
+max_string_len     = 50
+
+name_measvars      = "measured_variables"
+
+non_measurable_vars = ['station_name', 'time', 'time_bounds', 'crs', \
+                       'vrs', 'x', 'y', 'z', 'lon', 'lat', 'ntime', 'station', 'traj', \
+                       'E_UTM', 'N_UTM', 'height_above_origin', 'station_h', \
+                       'traj_name', 'height', 'band_pm_size', 'bands_pm', 'bands_pm_size_bounds' \
+                       'bands_pm_size', 'ancillary_detected_layer' ]
+
+soil_vars            = [ 't_soil', 'm_soil', 'lwc', 'lwcs', 'smp' ]
+
+dims_out             = [ name_eutm, name_nutm, name_hao, name_z, name_station_h ]
+
+# Define list of attributes which need to be of type float. In the data set this is not 
+# necessarily guranteed.
+atts_float           = [ 'origin_x', 'origin_y', 'origin_z', 'origin_lon', 'origin_lat', 'rotation_angle' ]
+
+# Define list of default variables that shall be measured at each site
+vars_default         = [ 'u', 'v', 'w', 'theta', 'hus' ]
+
+
+#Read config file
+read_config_file()
+
+# Initialize counter variable for the number of sites
+num_sites = 0
+
+# Set the output path for the data
+output_filename = output_path + output_filename
+
+# Open output file
+ncfile_out = Dataset( output_filename, "w", format="NETCDF4" )
+
+# First, add global attributes
+ncfile_out.setncattr( 'acronym',        global_acronym      )
+ncfile_out.setncattr( 'author',         global_author       )
+ncfile_out.setncattr( 'campaign',       global_campaign     )
+ncfile_out.setncattr( 'comment',        global_comment      )
+ncfile_out.setncattr( 'contact_person', global_contact      )
+ncfile_out.setncattr( 'data_content',   global_data_content )
+ncfile_out.setncattr( 'dependencies',   global_dependencies )
+ncfile_out.setncattr( 'institution',    global_institution  )
+ncfile_out.setncattr( 'keywords',       global_keywords     )
+ncfile_out.setncattr( 'location',       global_location     )
+ncfile_out.setncattr( 'references',     global_references   )
+ncfile_out.setncattr( 'site',           global_site         )
+ncfile_out.setncattr( 'source',         global_source       )
+ncfile_out.setncattr( 'palm_version',   global_palm_version )
+
+# Create universal dimension for the string length.
+ncfile_out.createDimension("string_len", max_string_len)
+
+
+# Check if observational data is available. This case, 
+# obtain an alphabetically sorted list of input data. List is sorted 
+# just for the sake of clarity in the resulting setup file.
+if ( input_from_observations == True ):
+   list_input_data = sorted( os.listdir( data_path ) )
+
+if ( input_from_observations ):
+
+   # Run loop over all subdirectories, detect the files and extract a list of sites.
+   # This is done to reduce the number of virtual measurements in the model. Each
+   # virtual measurement has an overhead and consumes memory.
+   sites = []
+   for dirname in list_input_data:
+       data_file = data_path + dirname
+
+       # Directory may contain various file versions.
+       # Take the one with highest cycle number.
+       highest_cycle_nr = 0
+       for filename in os.listdir(data_file):
+          start_seq = len( filename ) - 6
+          end_seq   = len( filename ) - 3
+          if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
+             highest_cycle_nr = int(filename[start_seq:end_seq])
+             latest_file      = filename
+
+       # Open the NetCDF file
+       input_file = data_file + "/" + latest_file
+       ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii')
+
+       # Read global attributes and write them immediately into the output file
+       for att in ncfile_in.ncattrs():
+          if ( att == name_site ):
+             site = ncfile_in.getncattr(att)
+
+       if ( site not in sites ):
+          sites.append(site)
+
+   # Define a flag array that is used to identify whether site dimensions are already
+   # defined or not.
+   create_metadata_for_site = [None] * len(sites)
+
+   # Define a nested list of default variables that shall be measured. Based on this list,
+   # the final number of measured variables is determined.
+   measured_variables_all_sites = [ ['u', 'v', 'w', 'theta', 'hus'] for var in range(0, len(sites))]
+
+   # Run loop over all subdirectories that contain observational data
+   for dirname in list_input_data:
+      data_file = data_path + dirname
+
+      # Directory may contain various file versions.
+      # Take the one with highest cycle number.
+      highest_cycle_nr = 0
+      for filename in os.listdir(data_file):
+         start_seq = len( filename ) - 6
+         end_seq   = len( filename ) - 3
+         if int( filename[start_seq:end_seq] ) > highest_cycle_nr:
+            highest_cycle_nr = int(filename[start_seq:end_seq])
+            latest_file      = filename
+
+      # Open the NetCDF input file
+      input_file = data_file + "/" + latest_file
+      ncfile_in = Dataset( input_file, "r", format="NETCDF4", encoding='ascii' )
+
+      # Read site attribue first
+      for att in ncfile_in.ncattrs():
+         if ( att == name_site ):
+            site = ncfile_in.getncattr(att)
+
+      # Determine index for the treated site
+      num_vmeas = sites.index( site ) + 1
+
+      # Check whether metadata for this site has been already created
+      if ( create_metadata_for_site[sites.index( site )] != "Done" ):
+
+         # Read global attributes and write them immediately into the output file
+         for att in ncfile_in.ncattrs():
+            if ( att == name_featuretype ):
+               feature = ncfile_in.getncattr(att)
+            if ( att == name_datacontent ):
+               content = ncfile_in.getncattr(att)
+            if ( att == name_site ):
+               site = ncfile_in.getncattr(att)
+
+            if ( att in atts_float ):
+               ncfile_out.setncattr( att + str(num_vmeas), np.double(ncfile_in.getncattr(att)) )
+            else:
+               ncfile_out.setncattr( att + str(num_vmeas), ncfile_in.getncattr(att) )
+
+         #timeSeries
+         if ( feature == name_ts ):
+            ntime = len( ncfile_in.dimensions[name_ntime]   )
+            nstat = len( ncfile_in.dimensions[name_station] )
+            ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
+            ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
+
+         #trajectory
+         elif ( feature == name_traj ):
+            ntime = len( ncfile_in.dimensions[name_ntime]   )
+            ntraj = len( ncfile_in.dimensions[name_traj_dim] )
+            ncfile_out.createDimension( name_ntime    + str(num_vmeas), ntime )
+            ncfile_out.createDimension( name_traj_dim + str(num_vmeas), ntraj )
+
+         #timeseriesProfile
+         else:
+            ntime = len( ncfile_in.dimensions[name_ntime]   )
+            nstat = len( ncfile_in.dimensions[name_station] )
+            nz    = len( ncfile_in.dimensions[name_nz]      )
+            ncfile_out.createDimension( name_ntime   + str(num_vmeas), ntime )
+            ncfile_out.createDimension( name_station + str(num_vmeas), nstat )
+            ncfile_out.createDimension( name_nz      + str(num_vmeas), nz    )
+
+         for var in ncfile_in.variables.keys():
+            if ( var in dims_out ):
+               # Create a variable and write it to file after it is read. In order to 
+               # avoid fill values in the dimensions, these are converted to zero 
+               # before written to file. Depending on the featureType of the measurement,
+               # the array shape is different. For more informations, please see 
+               # [UC]2 data standard.
+               # Timeseries
+               if ( feature == name_ts  ): 
+                  temp_ts = ncfile_out.createVariable( var + str(num_vmeas), float, \
+                                                       name_station + str(num_vmeas))
+                  temp_ts[:] = np.nan_to_num( ncfile_in.variables[var][:] )
+
+               # Trajectories
+               elif ( feature == name_traj ):
+                  temp_traj = ncfile_out.createVariable( var + str(num_vmeas), float, \
+                                                         ( name_traj_dim + str(num_vmeas), \
+                                                           name_ntime + str(num_vmeas) ) )
+                  temp_traj[:,:] = np.nan_to_num( ncfile_in.variables[var][:,:] )
+
+               # TimeseriesProfiles
+               else:
+                  if ( var == 'z' ):
+                     temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
+                                                         ( name_station + str(num_vmeas), \
+                                                           name_nz + str(num_vmeas) ) )
+                     temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:,0,:] )
+                  else:
+                     temp_pr = ncfile_out.createVariable( var + str(num_vmeas), float, \
+                                                          name_station + str(num_vmeas))
+                     temp_pr[:] = np.nan_to_num( ncfile_in.variables[var][:] )
+
+      # Search for variables to be measured. In case the variable isn't already defined,
+      # append the variable to the list.
+      for var in ncfile_in.variables.keys():
+         if ( var not in non_measurable_vars  and  \
+              var not in vars_default         and  \
+              var not in measured_variables_all_sites[sites.index( site )] ):
+
+            measured_variables_all_sites[sites.index( site )].append(var)
+
+      # Close the NetCDF input file
+      ncfile_in.close()
+
+      # Set flag to indicate that for this specific site dimensions have been 
+      # already created and attributes are already set.
+      if ( create_metadata_for_site[sites.index( site )] != "Done" ):
+         create_metadata_for_site[sites.index( site )] = "Done"
+
+   # After variables are gathered and dimensions / attributes are already written to file,
+   # the list of measured variables is written to file.
+   for site in sites:
+
+      num_vmeas = sites.index( site ) + 1
+
+      ncfile_out.createDimension( "nvar"+ str(num_vmeas), \
+                                  len( measured_variables_all_sites[sites.index( site )] ) )
+
+      measured = ncfile_out.createVariable( 'measured_variables' + str(num_vmeas), 'S1', \
+                                            ("nvar" + str(num_vmeas), "string_len")) # must be NC_CHAR
+
+      for counter, meas in enumerate( measured_variables_all_sites[sites.index( site )] ):
+         measured[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
+
+      # Check if any of the measured variables is a soil variable. Set flag accordingly.
+      soil = False
+      for var in measured_variables_all_sites[sites.index( site )]:
+         if ( var in soil_vars ):
+            soil = True
+
+      # Write soil flag
+      ncfile_out.setncattr( name_soil_sampling + str( num_vmeas), np.int8(soil) )
+
+
+
+    # Store the number of observational sites
+   num_sites += len( sites )
+
+
+#  Now process the customized input data. Please note, at the moment only timeseries are
+# are possible.
+if ( custom_coordinates ):
+
+   count_site = num_sites + 1
+   for coord in coordinates:
+      # Define mandatory attributes
+      ncfile_out.setncattr( name_featuretype + str(count_site),  \
+                            name_ts )
+      ncfile_out.setncattr( name_site        + str(count_site),  \
+                            "custom"         + str(count_site - num_sites) )
+      ncfile_out.setncattr( name_orig_x      + str(count_site),  \
+                            coord[0] )
+      ncfile_out.setncattr( name_orig_y      + str(count_site),  \
+                            coord[1] )
+      ncfile_out.setncattr( name_orig_z      + str(count_site),  \
+                            0.0 )
+
+      # Define dimensions
+      ntime = 1
+      nstat = 1
+      ncfile_out.createDimension( name_ntime   + str(count_site), ntime )
+      ncfile_out.createDimension( name_station + str(count_site), nstat )
+
+      # Define coordinate variables
+      temp_ts = ncfile_out.createVariable( name_eutm      + str(count_site), \
+                                           float,                            \
+                                           name_station   + str(count_site) )
+      temp_ts[:] = np.array( coord[0] )
+
+      temp_ts = ncfile_out.createVariable( name_nutm      + str(count_site), \
+                                           float,                            \
+                                           name_station   + str(count_site) )
+      temp_ts[:] = np.array( coord[1] )
+
+      temp_ts = ncfile_out.createVariable( name_z         + str(count_site), \
+                                           float,                            \
+                                           name_station   + str(count_site) )
+      temp_ts[:] = np.array( coord[2] )
+
+      temp_ts = ncfile_out.createVariable( name_station_h + str(count_site), \
+                                           float,                            \
+                                           name_station   + str(count_site) )
+      temp_ts[:] = np.array( 0.0 )
+
+
+      count_site += 1
+
+   # Reset counter variable
+   count_site = num_sites + 1
+
+   # check if variables are prescribed. If so, prepare final output string
+   # stored in measured_variables.
+   if ( vars_to_be_measured ):
+
+      for custom_vars in vars_to_be_measured:
+
+         measured_variables = []
+         for var in vars_default:
+            measured_variables.append(var)
+
+         # Check if given variables are already in the default variables. 
+         # If not, extend.
+         for var in custom_vars:
+             if ( var  not in  measured_variables ):
+
+                measured_variables.append(var)
+
+         ncfile_out.createDimension( "nvar"+ str(count_site), \
+                                     len( measured_variables ) )
+
+         measured_var = ncfile_out.createVariable( 'measured_variables' + str(count_site), 'S1', \
+                                                  ("nvar" + str(count_site), "string_len") ) # must be NC_CHAR
+
+         # Write the variables to the file
+         for counter, meas in enumerate( measured_variables ):
+            measured_var[counter] = stringtochar( np.array( meas,"S%s"%(max_string_len) ) )
+
+         # Add soil attribute for the current measurement.
+         soil = False
+         if ( any( var == soil_vars for var in measured_variables) ):
+            soil = True
+
+         # Write soil flag
+         ncfile_out.setncattr( name_soil_sampling + str( count_site), np.int8(soil) )
+
+         # Increment counter variable
+         count_site += 1
+
+         del ( measured_variables[:] )
+
+   # Add the number of customized sites.
+   num_sites += int( number_positions )
+
+
+# Finally, write the total number of sites to the output file
+ncfile_out.setncattr( name_num_stat, num_sites )
+
+
+print( "*** palm_cvd has been finished. You can find the output file under: " )
+print( "    " + output_filename )
+
+quit()
+
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 4399)
+++ /palm/trunk/SOURCE/Makefile	(revision 4400)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Add dependency for data-output module
+# 
+# 4392 2020-01-31 16:14:57Z pavelkrc
 # add dependency on fft for transpose
 # 
@@ -750,4 +753,5 @@
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o \
+	data_output_module.o \
 	diagnostic_output_quantities_mod.o \
 	dynamics_mod.o \
@@ -1273,12 +1277,16 @@
 	user_init_flight.o
 virtual_measurement_mod.o: \
+	basic_constants_and_equations_mod.o \
 	cpulog_mod.o \
 	chem_modules.o \
 	chem_gasphase_mod.o \
-	mod_kinds.o \
-	modules.o \
-	netcdf_data_input_mod.o \
-	netcdf_interface_mod.o \
-	radiation_model_mod.o
+	data_output_module.o \
+	land_surface_model_mod.o \
+	mod_kinds.o \
+	modules.o \
+	netcdf_data_input_mod.o \
+	radiation_model_mod.o \
+	surface_mod.o \
+	urban_surface_mod.o
 wind_turbine_model_mod.o: \
 	basic_constants_and_equations_mod.o \
Index: /palm/trunk/SOURCE/basic_constants_and_equations_mod.f90
===================================================================
--- /palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 4399)
+++ /palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 4400)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Move routine to transform coordinates from netcdf_interface_mod to 
+! basic_constants_and_equations_mod
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -113,4 +117,9 @@
             barometric_formula_1d
 
+
+    INTERFACE convert_utm_to_geographic
+       MODULE PROCEDURE convert_utm_to_geographic
+    END INTERFACE convert_utm_to_geographic
+
     INTERFACE magnus
        MODULE PROCEDURE magnus_0d
@@ -142,6 +151,146 @@
        MODULE PROCEDURE barometric_formula_1d
     END INTERFACE barometric_formula
+!
+!-- Public routines
+    PUBLIC convert_utm_to_geographic
 
  CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Convert UTM coordinates into geographic latitude and longitude. Conversion
+!> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
+!> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
+!> Based on a JavaScript of the geodesy function library written by chrisveness
+!> https://github.com/chrisveness/geodesy
+!------------------------------------------------------------------------------!
+    SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
+
+       INTEGER(iwp) ::  j   !< loop index
+
+       REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
+       REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
+       REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
+       REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
+
+       REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
+       REAL(wp) ::  cos_eta_s   !< cos(eta_s)
+       REAL(wp) ::  delta_i     !<
+       REAL(wp) ::  delta_tau_i !<
+       REAL(wp) ::  e           !< eccentricity
+       REAL(wp) ::  eta         !<
+       REAL(wp) ::  eta_s       !<
+       REAL(wp) ::  n           !< 3rd flattening
+       REAL(wp) ::  n2          !< n^2
+       REAL(wp) ::  n3          !< n^3
+       REAL(wp) ::  n4          !< n^4
+       REAL(wp) ::  n5          !< n^5
+       REAL(wp) ::  n6          !< n^6
+       REAL(wp) ::  nu          !<
+       REAL(wp) ::  nu_s        !<
+       REAL(wp) ::  sin_eta_s   !< sin(eta_s)
+       REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
+       REAL(wp) ::  tau_i       !<
+       REAL(wp) ::  tau_i_s     !<
+       REAL(wp) ::  tau_s       !<
+       REAL(wp) ::  x           !< adjusted easting
+       REAL(wp) ::  y           !< adjusted northing
+
+       REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
+
+       REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
+                                                  !< (/semi_major_axis,
+                                                  !<   inverse_flattening,
+                                                  !<   longitude_of_prime_meridian,
+                                                  !<   longitude_of_central_meridian,
+                                                  !<   scale_factor_at_central_meridian,
+                                                  !<   latitude_of_projection_origin,
+                                                  !<   false_easting,
+                                                  !<   false_northing /)
+
+       x = eutm - crs(7)  ! remove false easting
+       y = nutm - crs(8)  ! remove false northing
+!
+!--    from Karney 2011 Eq 15-22, 36:
+       e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
+       n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
+       n2 = n * n
+       n3 = n * n2
+       n4 = n * n3
+       n5 = n * n4
+       n6 = n * n5
+
+       a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
+                                              + 0.015625_wp * n4   &
+                                              + 3.90625E-3_wp * n6 )
+
+       nu  = x / ( crs(5) * a )
+       eta = y / ( crs(5) * a )
+
+!--    According to KrÃŒger (1912), eq. 26*
+       beta(1) =        0.5_wp                  * n  &
+               -        2.0_wp /         3.0_wp * n2 &
+               +       37.0_wp /        96.0_wp * n3 &
+               -        1.0_wp /       360.0_wp * n4 &
+               -       81.0_wp /       512.0_wp * n5 &
+               +    96199.0_wp /    604800.0_wp * n6
+
+       beta(2) =        1.0_wp /        48.0_wp * n2 &
+               +        1.0_wp /        15.0_wp * n3 &
+               -      437.0_wp /      1440.0_wp * n4 &
+               +       46.0_wp /       105.0_wp * n5 &
+               -  1118711.0_wp /   3870720.0_wp * n6
+
+       beta(3) =       17.0_wp /       480.0_wp * n3 &
+               -       37.0_wp /       840.0_wp * n4 &
+               -      209.0_wp /      4480.0_wp * n5 &
+               +     5569.0_wp /     90720.0_wp * n6
+
+       beta(4) =     4397.0_wp /    161280.0_wp * n4 &
+               -       11.0_wp /       504.0_wp * n5 &
+               -   830251.0_wp /   7257600.0_wp * n6
+
+       beta(5) =     4583.0_wp /    161280.0_wp * n5 &
+               -   108847.0_wp /   3991680.0_wp * n6
+
+       beta(6) = 20648693.0_wp / 638668800.0_wp * n6
+
+       eta_s = eta
+       nu_s  = nu
+       DO  j = 1, 6
+         eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
+         nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
+       ENDDO
+
+       sinh_nu_s = SINH( nu_s )
+       sin_eta_s = SIN( eta_s )
+       cos_eta_s = COS( eta_s )
+
+       tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
+
+       tau_i = tau_s
+       delta_tau_i = 1.0_wp
+
+       DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
+
+         delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
+
+         tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
+                  - delta_i * SQRT( 1.0_wp + tau_i**2 )
+
+         delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
+                      * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
+                      / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
+
+         tau_i = tau_i + delta_tau_i
+
+       ENDDO
+
+       lat = ATAN( tau_i ) / pi * 180.0_wp
+       lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
+
+    END SUBROUTINE convert_utm_to_geographic
 
 !------------------------------------------------------------------------------!
@@ -161,5 +310,5 @@
 !--    Saturation vapor pressure for a specific temperature:
        magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) /             &
-                                                   ( t - 29.65_wp  ) )
+                                               ( t - 29.65_wp  ) )
 
     END FUNCTION magnus_0d
Index: /palm/trunk/SOURCE/check_open.f90
===================================================================
--- /palm/trunk/SOURCE/check_open.f90	(revision 4399)
+++ /palm/trunk/SOURCE/check_open.f90	(revision 4400)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Remove binary output for virtual measurements
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -444,37 +447,4 @@
           ENDIF
 
-       CASE ( 27 )
-!
-!--       Binary files for virtual measurement data
-          IF ( myid_char == '' )  THEN
-             OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM( coupling_char )//       &
-                             myid_char, FORM='UNFORMATTED', POSITION='APPEND' )
-          ELSE
-
-             IF ( myid == 0  .AND. .NOT. openfile(file_id)%opened_before )  THEN
-                CALL local_system( 'mkdir  VIRTUAL_MEAS_BIN' //                &
-                                   TRIM( coupling_char ) )
-             ENDIF
-#if defined( __parallel )
-!
-!--          Set a barrier in order to allow that all other processors in the 
-!--          directory created by PE0 can open their file
-             CALL MPI_BARRIER( comm2d, ierr )
-#endif
-             ioerr = 1
-             DO WHILE ( ioerr /= 0 )
-                OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM(coupling_char)//      &
-                                '/'//myid_char,                                &
-                           FORM='UNFORMATTED', IOSTAT=ioerr )
-                IF ( ioerr /= 0 )  THEN
-                   WRITE( 9, * )  '*** could not open "VIRTUAL_MEAS_BIN'//     &
-                                  TRIM(coupling_char)//'/'//myid_char//        &
-                                  '"! Trying again in 1 sec.'
-                   CALL fortran_sleep( 1 )
-                ENDIF
-             ENDDO
-
-          ENDIF
-
        CASE ( 30 )
 
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 4399)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 4400)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - Use data-output module for virtual measurement output
+! - Remove deprecated routines for virtual measurement module
+! 
+! 4361 2020-01-07 12:22:38Z suehring
 ! Remove unused arrays in pmc_rrd_local
 ! 
@@ -172,7 +176,14 @@
     USE kinds
 
+    USE pegrid,                                                                &
+        ONLY:  comm2d
+
 !
 !-- load module-specific control parameters.
 !-- ToDo: move all of them to respective module or a dedicated central module
+    USE data_output_module,                                                    &
+        ONLY:  dom_def_end,                                                    &
+               dom_finalize_output,                                            &
+               dom_init
 
     USE dynamics_mod, &
@@ -220,4 +231,5 @@
         ONLY:  air_chemistry,                                                  &
                biometeorology,                                                 &
+               coupling_char,                                                  &
                debug_output,                                                   &
                debug_output_timestep,                                          &
@@ -493,5 +505,5 @@
         ONLY:  vm_check_parameters,                                            &
                vm_init,                                                        &
-               vm_last_actions,                                                &
+               vm_init_output,                                                 &
                vm_parin
 
@@ -543,4 +555,5 @@
        module_interface_init,                                                  &
        module_interface_init_checks,                                           &
+       module_interface_init_output,                                           &
        module_interface_header,                                                &
        module_interface_actions,                                               &
@@ -600,4 +613,8 @@
        MODULE PROCEDURE module_interface_init_checks
     END INTERFACE module_interface_init_checks
+
+    INTERFACE module_interface_init_output
+       MODULE PROCEDURE module_interface_init_output
+    END INTERFACE module_interface_init_output
 
     INTERFACE module_interface_header
@@ -1081,7 +1098,29 @@
     IF ( debug_output )  CALL debug_message( 'module-specific initialization', 'end' )
 
-
  END SUBROUTINE module_interface_init
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize data output
+!------------------------------------------------------------------------------!
+ SUBROUTINE module_interface_init_output
+
+    INTEGER(iwp) ::  return_value  !< returned status value of called function
+
+!
+!-- Initialize data-output module
+    CALL dom_init( file_suffix_of_output_group=coupling_char,                  &
+                   mpi_comm_of_output_group=comm2d,                            &
+                   program_debug_output_unit=6,                                &
+                   debug_output=debug_output )
+!
+!-- Define module-specific output quantities
+    IF ( virtual_measurement )  CALL vm_init_output
+!
+!-- Leave output-definition state
+    return_value = dom_def_end()
+
+ END SUBROUTINE module_interface_init_output
 
 !------------------------------------------------------------------------------!
@@ -1890,11 +1929,13 @@
  SUBROUTINE module_interface_last_actions
 
+    INTEGER ::  return_value  !< returned status value of a called function
+
 
     IF ( debug_output )  CALL debug_message( 'module-specific last actions', 'start' )
 
+    return_value = dom_finalize_output()
+
     CALL dynamics_last_actions
 
-    IF ( virtual_measurement )  CALL vm_last_actions
-
     IF ( user_module_enabled )  CALL user_last_actions
 
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4399)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4400)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - Routine to inquire default fill values added
+! - netcdf_data_input_att and netcdf_data_input_var routines removed
+! 
+! 4392 2020-01-31 16:14:57Z pavelkrc
 ! (resler) Decrease length of reading buffer (fix problem of ifort/icc compilers)
 ! 
@@ -542,4 +546,5 @@
        INTEGER(iwp) ::  version              !< version of data set
 
+       REAL(wp) ::  fillvalue = -9999.0      !< default fill value
        REAL(wp) ::  origin_lat               !< latitude of lower left corner
        REAL(wp) ::  origin_lon               !< longitude of lower left corner
@@ -652,4 +657,6 @@
     LOGICAL ::  collective_read = .FALSE.      !< Enable NetCDF collective read
 
+    REAL(wp), DIMENSION(8) ::  crs_list        !< list of coord_ref_sys values
+
     TYPE(global_atts_type) ::  input_file_atts !< global attributes of input file
 
@@ -666,11 +673,16 @@
     END INTERFACE netcdf_data_input_check_static
 
-    INTERFACE netcdf_data_input_chemistry_data                        
+    INTERFACE netcdf_data_input_chemistry_data
        MODULE PROCEDURE netcdf_data_input_chemistry_data
     END INTERFACE netcdf_data_input_chemistry_data
-    
-    INTERFACE get_dimension_length                        
+
+    INTERFACE get_dimension_length
        MODULE PROCEDURE get_dimension_length
     END INTERFACE get_dimension_length
+
+    INTERFACE inquire_fill_value
+       MODULE PROCEDURE inquire_fill_value_int
+       MODULE PROCEDURE inquire_fill_value_real
+    END INTERFACE inquire_fill_value
 
     INTERFACE netcdf_data_input_inquire_file
@@ -681,11 +693,4 @@
        MODULE PROCEDURE netcdf_data_input_init
     END INTERFACE netcdf_data_input_init
-    
-    INTERFACE netcdf_data_input_att
-       MODULE PROCEDURE netcdf_data_input_att_int8
-       MODULE PROCEDURE netcdf_data_input_att_int32
-       MODULE PROCEDURE netcdf_data_input_att_real
-       MODULE PROCEDURE netcdf_data_input_att_string
-    END INTERFACE netcdf_data_input_att
 
     INTERFACE netcdf_data_input_init_3d
@@ -696,10 +701,4 @@
        MODULE PROCEDURE netcdf_data_input_surface_data
     END INTERFACE netcdf_data_input_surface_data
-
-    INTERFACE netcdf_data_input_var
-       MODULE PROCEDURE netcdf_data_input_var_char
-       MODULE PROCEDURE netcdf_data_input_var_real_1d
-       MODULE PROCEDURE netcdf_data_input_var_real_2d
-    END INTERFACE netcdf_data_input_var
 
     INTERFACE netcdf_data_input_uvem
@@ -752,4 +751,5 @@
            chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
            coord_ref_sys,                                                      &
+           crs_list,                                                           &
            init_3d, init_model, input_file_atts,                               &
            input_file_dynamic,                                                 &
@@ -781,8 +781,6 @@
            netcdf_data_input_init,                                             &
            netcdf_data_input_init_3d,                                          &
-           netcdf_data_input_att,                                              &
            netcdf_data_input_surface_data,                                     &
            netcdf_data_input_topo,                                             &
-           netcdf_data_input_var,                                              &
            get_attribute,                                                      &
            get_variable,                                                       &
@@ -790,4 +788,5 @@
            open_read_file,                                                     &
            check_existence,                                                    &
+           inquire_fill_value,                                                 &
            inquire_num_variables,                                              &
            inquire_variable_names,                                             &
@@ -817,5 +816,5 @@
        INQUIRE( FILE = TRIM( input_file_chem )    // TRIM( coupling_char ),    &
                 EXIST = input_pids_chem )
-       INQUIRE( FILE = TRIM( input_file_uvem ) // TRIM( coupling_char ),       &
+       INQUIRE( FILE = TRIM( input_file_uvem )    // TRIM( coupling_char ),    &
                 EXIST = input_pids_uvem  )
        INQUIRE( FILE = TRIM( input_file_vm )      // TRIM( coupling_char ),    &
@@ -962,5 +961,16 @@
        init_model%origin_z        = input_file_atts%origin_z  
        init_model%rotation_angle  = input_file_atts%rotation_angle  
-            
+
+!
+!--    Define list of crs attributes. This is required for coordinate 
+!--    transformation. 
+       crs_list = (/ coord_ref_sys%semi_major_axis,                            &
+                     coord_ref_sys%inverse_flattening,                         &
+                     coord_ref_sys%longitude_of_prime_meridian,                &
+                     coord_ref_sys%longitude_of_central_meridian,              &
+                     coord_ref_sys%scale_factor_at_central_meridian,           &
+                     coord_ref_sys%latitude_of_projection_origin,              &
+                     coord_ref_sys%false_easting,                              &
+                     coord_ref_sys%false_northing /)
 !
 !--    In case of nested runs, each model domain might have different longitude
@@ -977,258 +987,5 @@
 
     END SUBROUTINE netcdf_data_input_init
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read an array of characters.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_var_char( val, search_string, id_mod )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string     !< name of the variable
-       CHARACTER(LEN=*), DIMENSION(:) ::  val !< variable which should be read
-       
-       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
-
-#if defined ( __netcdf )
-!
-!--    Read variable
-       CALL get_variable( id_mod, search_string, val )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_var_char
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read an 1D array of REAL values.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_var_real_1d( val, search_string, id_mod )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string     !< name of the variable      
-       
-       INTEGER(iwp) ::  id_mod        !< NetCDF id of input file
-       
-       REAL(wp), DIMENSION(:) ::  val !< variable which should be read
-
-#if defined ( __netcdf )
-!
-!--    Read variable
-       CALL get_variable( id_mod, search_string, val )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_var_real_1d
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read an 1D array of REAL values.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_var_real_2d( val, search_string,              &
-                                              id_mod, d1s, d1e, d2s, d2e )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string     !< name of the variable     
-       
-       INTEGER(iwp) ::  id_mod  !< NetCDF id of input file
-       INTEGER(iwp) ::  d1e     !< end index of first dimension to be read
-       INTEGER(iwp) ::  d2e     !< end index of second dimension to be read
-       INTEGER(iwp) ::  d1s     !< start index of first dimension to be read
-       INTEGER(iwp) ::  d2s     !< start index of second dimension to be read
-       
-       REAL(wp), DIMENSION(:,:) ::  val !< variable which should be read
-
-#if defined ( __netcdf )
-!
-!--    Read character variable
-       CALL get_variable( id_mod, search_string, val, d1s, d1e, d2s, d2e )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_var_real_2d
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read a global string attribute 
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_att_string( val, search_string, id_mod,       &
-                                             input_file, global, openclose,    &
-                                             variable_name )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string !< name of the attribue
-       CHARACTER(LEN=*) ::  val           !< attribute
-       
-       CHARACTER(LEN=*) ::  input_file    !< name of input file
-       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
-       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed  
-       
-       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
-       
-       LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
-
-#if defined ( __netcdf )
-!
-!--    Open file in read-only mode if necessary
-       IF ( openclose == 'open' )  THEN
-          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
-                                  id_mod )
-       ENDIF
-!
-!--    Read global attribute
-       IF ( global )  THEN
-          CALL get_attribute( id_mod, search_string, val, global )
-!
-!--    Read variable attribute
-       ELSE
-          CALL get_attribute( id_mod, search_string, val, global, variable_name )
-       ENDIF
-!
-!--    Close input file 
-       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_att_string
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read a global 8-bit integer attribute 
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_att_int8( val, search_string, id_mod,         &
-                                           input_file, global, openclose,      &
-                                           variable_name )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string !< name of the attribue
-       
-       CHARACTER(LEN=*) ::  input_file    !< name of input file
-       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
-       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
-       
-       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
-       INTEGER(KIND=1) ::  val      !< value of the attribute
-       
-       LOGICAL ::  global        !< flag indicating a global or a variable's attribute
-
-#if defined ( __netcdf )
-!
-!--    Open file in read-only mode
-       IF ( openclose == 'open' )  THEN
-          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
-                                  id_mod )
-       ENDIF
-!
-!--    Read global attribute
-       IF ( global )  THEN
-          CALL get_attribute( id_mod, search_string, val, global )
-!
-!--    Read variable attribute
-       ELSE
-          CALL get_attribute( id_mod, search_string, val, global, variable_name )
-       ENDIF
-!
-!--    Finally, close input file
-       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_att_int8
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read a global 32-bit integer attribute 
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_att_int32( val, search_string, id_mod,        &
-                                            input_file, global, openclose,     &
-                                            variable_name )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string !< name of the attribue
-       
-       CHARACTER(LEN=*) ::  input_file    !< name of input file
-       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
-       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
-       
-       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
-       INTEGER(iwp) ::  val      !< value of the attribute
-       
-       LOGICAL ::  global        !< flag indicating a global or a variable's attribute
-
-#if defined ( __netcdf )
-!
-!--    Open file in read-only mode
-       IF ( openclose == 'open' )  THEN
-          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
-                                  id_mod )
-       ENDIF
-!
-!--    Read global attribute
-       IF ( global )  THEN
-          CALL get_attribute( id_mod, search_string, val, global )
-!
-!--    Read variable attribute
-       ELSE
-          CALL get_attribute( id_mod, search_string, val, global, variable_name )
-       ENDIF
-!
-!--    Finally, close input file
-       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_att_int32
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Read a global real attribute 
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_att_real( val, search_string, id_mod,         &
-                                           input_file, global, openclose,      &
-                                           variable_name )
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*) ::  search_string !< name of the attribue
-       
-       CHARACTER(LEN=*) ::  input_file    !< name of input file
-       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
-       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
-       
-       INTEGER(iwp) ::  id_mod            !< NetCDF id of input file
-       
-       LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
-       
-       REAL(wp) ::  val                   !< value of the attribute
-
-#if defined ( __netcdf )
-!
-!--    Open file in read-only mode
-       IF ( openclose == 'open' )  THEN
-          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
-                                  id_mod )
-       ENDIF
-!
-!--    Read global attribute
-       IF ( global )  THEN
-          CALL get_attribute( id_mod, search_string, val, global )
-!
-!--    Read variable attribute
-       ELSE
-          CALL get_attribute( id_mod, search_string, val, global, variable_name )
-       ENDIF
-!
-!--    Finally, close input file
-       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
-#endif            
-
-    END SUBROUTINE netcdf_data_input_att_real
+
 
 !------------------------------------------------------------------------------!
@@ -5823,4 +5580,53 @@
 ! Description:
 ! ------------
+!> Inquires the _FillValue settings of an integer variable.
+!------------------------------------------------------------------------------!
+    SUBROUTINE inquire_fill_value_int( id, var_name, nofill, fill_value )
+#if defined( __netcdf )
+       CHARACTER(LEN=*), INTENT(IN) ::  var_name    !< variable name
+
+       INTEGER(iwp), INTENT(IN)  ::  id          !< file id
+       INTEGER(iwp)              ::  nofill      !< flag indicating whether fill values are set or not
+       INTEGER(iwp)              ::  fill_value  !< fill value
+       INTEGER(iwp)              ::  id_var      !< netCDF variable id (varid)
+
+       nc_stat = NF90_INQ_VARID( id, TRIM( var_name ), id_var )
+       nc_stat = NF90_INQ_VAR_FILL(id, id_var, no_fill, fill_value )
+#endif
+!
+!--    Further line is just to avoid compiler warnings. nofill might be used
+!--    in future.
+       IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
+
+    END SUBROUTINE inquire_fill_value_int
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Inquires the _FillValue settings of a real variable.
+!------------------------------------------------------------------------------!
+    SUBROUTINE inquire_fill_value_real( id, var_name, nofill, fill_value )
+#if defined( __netcdf )
+       CHARACTER(LEN=*), INTENT(IN) ::  var_name    !< variable name
+
+       INTEGER(iwp), INTENT(IN)  ::  id          !< file id
+       INTEGER(iwp)              ::  nofill      !< flag indicating whether fill values are set or not
+       INTEGER(iwp)              ::  id_var      !< netCDF variable id (varid)
+
+       REAL(wp), INTENT(OUT)     ::  fill_value  !< fill value
+
+       nc_stat = NF90_INQ_VARID( id, TRIM( var_name ), id_var )
+       nc_stat = NF90_INQ_VAR_FILL(id, id_var, no_fill, fill_value )
+#endif
+!
+!--    Further line is just to avoid compiler warnings. nofill might be used
+!--    in future.
+       IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
+
+    END SUBROUTINE inquire_fill_value_real
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Prints out a text message corresponding to the current status.
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4399)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4400)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Move routine to transform coordinates from netcdf_interface_mod to 
+! basic_constants_and_equations_mod
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Adjusted output of multi-agent system for biometeorology
 ! 
@@ -135,5 +139,7 @@
 
     USE netcdf_data_input_mod,                                                 &
-        ONLY: coord_ref_sys, init_model
+        ONLY: coord_ref_sys,                                                   &
+              crs_list,                                                        &
+              init_model
 
     PRIVATE
@@ -417,5 +423,6 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  pi
+        ONLY:  convert_utm_to_geographic,                                      &
+               pi
 
     USE control_parameters,                                                    &
@@ -564,5 +571,4 @@
 
     REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
-    REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
 
     REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
@@ -662,14 +668,4 @@
 
     ENDIF
-!
-!-- Convert coord_ref_sys into vector (used for lat/lon calculation)
-    crs_list = (/ coord_ref_sys%semi_major_axis,                  &
-                  coord_ref_sys%inverse_flattening,               &
-                  coord_ref_sys%longitude_of_prime_meridian,      &
-                  coord_ref_sys%longitude_of_central_meridian,    &
-                  coord_ref_sys%scale_factor_at_central_meridian, &
-                  coord_ref_sys%latitude_of_projection_origin,    &
-                  coord_ref_sys%false_easting,                    &
-                  coord_ref_sys%false_northing /)
 
 !
@@ -7108,144 +7104,3 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Convert UTM coordinates into geographic latitude and longitude. Conversion
-!> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
-!> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
-!> Based on a JavaScript of the geodesy function library written by chrisveness
-!> https://github.com/chrisveness/geodesy
-!------------------------------------------------------------------------------!
- SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
-
-    USE basic_constants_and_equations_mod,                                     &
-        ONLY:  pi
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  j   !< loop index
-
-    REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
-    REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
-    REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
-    REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
-
-    REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
-    REAL(wp) ::  cos_eta_s   !< cos(eta_s)
-    REAL(wp) ::  delta_i     !<
-    REAL(wp) ::  delta_tau_i !<
-    REAL(wp) ::  e           !< eccentricity
-    REAL(wp) ::  eta         !<
-    REAL(wp) ::  eta_s       !<
-    REAL(wp) ::  n           !< 3rd flattening
-    REAL(wp) ::  n2          !< n^2
-    REAL(wp) ::  n3          !< n^3
-    REAL(wp) ::  n4          !< n^4
-    REAL(wp) ::  n5          !< n^5
-    REAL(wp) ::  n6          !< n^6
-    REAL(wp) ::  nu          !<
-    REAL(wp) ::  nu_s        !<
-    REAL(wp) ::  sin_eta_s   !< sin(eta_s)
-    REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
-    REAL(wp) ::  tau_i       !<
-    REAL(wp) ::  tau_i_s     !<
-    REAL(wp) ::  tau_s       !<
-    REAL(wp) ::  x           !< adjusted easting
-    REAL(wp) ::  y           !< adjusted northing
-
-    REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
-
-    REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
-                                               !< (/semi_major_axis,
-                                               !<   inverse_flattening,
-                                               !<   longitude_of_prime_meridian,
-                                               !<   longitude_of_central_meridian,
-                                               !<   scale_factor_at_central_meridian,
-                                               !<   latitude_of_projection_origin,
-                                               !<   false_easting,
-                                               !<   false_northing /)
-
-    x = eutm - crs(7)  ! remove false easting
-    y = nutm - crs(8)  ! remove false northing
-
-!-- from Karney 2011 Eq 15-22, 36:
-    e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
-    n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
-    n2 = n * n
-    n3 = n * n2
-    n4 = n * n3
-    n5 = n * n4
-    n6 = n * n5
-
-    a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
-                                           + 0.015625_wp * n4   &
-                                           + 3.90625E-3_wp * n6 )
-
-    nu  = x / ( crs(5) * a )
-    eta = y / ( crs(5) * a )
-
-!-- According to KrÃŒger (1912), eq. 26*
-    beta(1) =        0.5_wp                  * n  &
-            -        2.0_wp /         3.0_wp * n2 &
-            +       37.0_wp /        96.0_wp * n3 &
-            -        1.0_wp /       360.0_wp * n4 &
-            -       81.0_wp /       512.0_wp * n5 &
-            +    96199.0_wp /    604800.0_wp * n6
-
-    beta(2) =        1.0_wp /        48.0_wp * n2 &
-            +        1.0_wp /        15.0_wp * n3 &
-            -      437.0_wp /      1440.0_wp * n4 &
-            +       46.0_wp /       105.0_wp * n5 &
-            -  1118711.0_wp /   3870720.0_wp * n6
-
-    beta(3) =       17.0_wp /       480.0_wp * n3 &
-            -       37.0_wp /       840.0_wp * n4 &
-            -      209.0_wp /      4480.0_wp * n5 &
-            +     5569.0_wp /     90720.0_wp * n6
-
-    beta(4) =     4397.0_wp /    161280.0_wp * n4 &
-            -       11.0_wp /       504.0_wp * n5 &
-            -   830251.0_wp /   7257600.0_wp * n6
-
-    beta(5) =     4583.0_wp /    161280.0_wp * n5 &
-            -   108847.0_wp /   3991680.0_wp * n6
-
-    beta(6) = 20648693.0_wp / 638668800.0_wp * n6
-
-    eta_s = eta
-    nu_s  = nu
-    DO  j = 1, 6
-      eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
-      nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
-    ENDDO
-
-    sinh_nu_s = SINH( nu_s )
-    sin_eta_s = SIN( eta_s )
-    cos_eta_s = COS( eta_s )
-
-    tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
-
-    tau_i = tau_s
-    delta_tau_i = 1.0_wp
-
-    DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
-
-      delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
-
-      tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
-               - delta_i * SQRT( 1.0_wp + tau_i**2 )
-
-      delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
-                   * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
-                   / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
-
-      tau_i = tau_i + delta_tau_i
-
-    ENDDO
-
-    lat = ATAN( tau_i ) / pi * 180.0_wp
-    lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
-
- END SUBROUTINE convert_utm_to_geographic
-
  END MODULE netcdf_interface
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 4399)
+++ /palm/trunk/SOURCE/palm.f90	(revision 4400)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add interface to initialize data output with dom
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! implement new palm_date_time_mod
 ! 
@@ -90,5 +93,7 @@
 
     USE module_interface,                                                      &
-        ONLY:  module_interface_last_actions
+        ONLY:  module_interface_init_output,                                   &
+               module_interface_last_actions
+
 
     USE multi_agent_system_mod,                                                &
@@ -269,4 +274,6 @@
 
     CALL init_3d_model
+
+    CALL module_interface_init_output
 
 !
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4399)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4400)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Initialize radiation arrays with zero
+! 
+! 4392 2020-01-31 16:14:57Z pavelkrc
 ! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
 ! - Print exact counts of SVF and CSF if debut_output is enabled
@@ -1782,14 +1785,18 @@
           IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
              ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_sw_in = 0.0_wp
           ENDIF
           IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
              ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_sw_out = 0.0_wp
           ENDIF
 
           IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
              ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_lw_in = 0.0_wp
           ENDIF
           IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
              ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_lw_out = 0.0_wp
           ENDIF
 
Index: /palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- /palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4399)
+++ /palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4400)
@@ -25,58 +25,65 @@
 ! -----------------
 ! $Id$
+! Revision of the module:
+! - revised input from NetCDF setup file
+! - parallel NetCDF output via data-output module ( Tobias Gronemeier )
+! - variable attributes added
+! - further variables defined
+!
+! 4346 2019-12-18 11:55:56Z motisi
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
-! 
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4226 2019-09-10 17:03:24Z suehring
 ! Netcdf input routine for dimension length renamed
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4168 2019-08-16 13:50:17Z suehring
 ! Replace function get_topography_top_index by topo_top_ind
-! 
+!
 ! 3988 2019-05-22 11:32:37Z kanani
 ! Add variables to enable steering of output interval for virtual measurements
-! 
+!
 ! 3913 2019-04-17 15:12:28Z gronemeier
 ! Bugfix: rotate positions of measurements before writing them into file
-! 
+!
 ! 3910 2019-04-17 11:46:56Z suehring
 ! Bugfix in rotation of UTM coordinates
-! 
+!
 ! 3904 2019-04-16 18:22:51Z gronemeier
 ! Rotate coordinates of stations by given rotation_angle
-! 
+!
 ! 3876 2019-04-08 18:41:49Z knoop
 ! Remove print statement
-! 
+!
 ! 3854 2019-04-02 16:59:33Z suehring
 ! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
-! nspec by nvar 
-! 
+! nspec by nvar
+!
 ! 3766 2019-02-26 16:23:41Z raasch
 ! unused variables removed
-! 
+!
 ! 3718 2019-02-06 11:08:28Z suehring
 ! Adjust variable name connections between UC2 and chemistry variables
-! 
+!
 ! 3717 2019-02-05 17:21:16Z suehring
 ! Additional check + error numbers adjusted
-! 
+!
 ! 3706 2019-01-29 20:02:26Z suehring
 ! unused variables removed
-! 
+!
 ! 3705 2019-01-29 19:56:39Z suehring
 ! - initialization revised
 ! - binary data output
 ! - list of allowed variables extended
-! 
+!
 ! 3704 2019-01-29 19:51:41Z suehring
 ! Sampling of variables
-! 
+!
 ! 3473 2018-10-30 20:50:15Z suehring
 ! Initial revision
@@ -85,29 +92,45 @@
 ! --------
 ! @author Matthias Suehring
+! @author Tobias Gronemeier
 !
 ! Description:
 ! ------------
-!> The module acts as an interface between 'real-world' observations and 
+!> The module acts as an interface between 'real-world' observations and
 !> model simulations. Virtual measurements will be taken in the model at the
-!> coordinates representative for the 'real-world' observation coordinates. 
+!> coordinates representative for the 'real-world' observation coordinates.
 !> More precisely, coordinates and measured quanties will be read from a
-!> NetCDF file which contains all required information. In the model, 
+!> NetCDF file which contains all required information. In the model,
 !> the same quantities (as long as all the required components are switched-on)
 !> will be sampled at the respective positions and output into an extra file,
-!> which allows for straight-forward comparison of model results with 
-!> observations. 
+!> which allows for straight-forward comparison of model results with
+!> observations.
 !>
-!> @todo list_of_allowed variables needs careful checking
-!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
-!>       the (UC)2 standard 
-!> @note Fluxes are not processed
+!> @todo Check why there is an error when _FillValue attributes are added via
+!>       dom.
+!> @todo Output of character variable station_name (dom hasn't this feature
+!>       yet implemented).
 !------------------------------------------------------------------------------!
  MODULE virtual_measurement_mod
 
     USE arrays_3d,                                                             &
-        ONLY:  q, pt, u, v, w, zu, zw
+        ONLY:  dzw,                                                            &
+               exner,                                                          &
+               hyp,                                                            &
+               q,                                                              &
+               ql,                                                             &
+               pt,                                                             &
+               rho_air,                                                        &
+               u,                                                              &
+               v,                                                              &
+               w,                                                              &
+               zu,                                                             &
+               zw
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  pi
+        ONLY:  convert_utm_to_geographic,                                      &
+               degc_to_k,                                                      &
+               magnus,                                                         &
+               pi,                                                             &
+               rd_d_rv
 
     USE chem_gasphase_mod,                                                     &
@@ -116,98 +139,175 @@
     USE chem_modules,                                                          &
         ONLY:  chem_species
-        
+
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, dz, humidity, io_blocks, io_group, neutral,      &
-               message_string, time_since_reference_point, virtual_measurement
+        ONLY:  air_chemistry,                                                  &
+               child_domain,                                                   &
+               coupling_char,                                                  &
+               dz,                                                             &
+               end_time,                                                       &
+               humidity,                                                       &
+               message_string,                                                 &
+               neutral,                                                        &
+               origin_date_time,                                               &
+               rho_surface,                                                    &
+               surface_pressure,                                               &
+               time_since_reference_point,                                     &
+               virtual_measurement,                                            &
+               write_binary
 
     USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point
-        
+        ONLY:  cpu_log,                                                        &
+               log_point
+
+    USE data_output_module
+
     USE grid_variables,                                                        &
-        ONLY:  dx, dy
+        ONLY:  ddx,                                                            &
+               ddy,                                                            &
+               dx,                                                             &
+               dy
 
     USE indices,                                                               &
-        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind,             &
+        ONLY:  nbgp,                                                           &
+               nzb,                                                            &
+               nzt,                                                            &
+               nxl,                                                            &
+               nxlg,                                                           &
+               nxr,                                                            &
+               nxrg,                                                           &
+               nys,                                                            &
+               nysg,                                                           &
+               nyn,                                                            &
+               nyng,                                                           &
+               topo_top_ind,                                                   &
                wall_flags_total_0
 
     USE kinds
-    
+
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  init_model
-        
+        ONLY:  close_input_file,                                               &
+               coord_ref_sys,                                                  &
+               crs_list,                                                       &
+               get_attribute,                                                  &
+               get_dimension_length,                                           &
+               get_variable,                                                   &
+               init_model,                                                     &
+               input_file_atts,                                                &
+               input_file_vm,                                                  &
+               input_pids_static,                                              &
+               input_pids_vm,                                                  &
+               inquire_fill_value,                                             &
+               open_read_file,                                                 &
+               pids_id
+
     USE pegrid
-    
+
     USE surface_mod,                                                           &
-        ONLY:  surf_lsm_h, surf_usm_h
-        
+        ONLY:  surf_lsm_h,                                                     &
+               surf_usm_h
+
     USE land_surface_model_mod,                                                &
-        ONLY:  nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h 
-        
-    USE radiation_model_mod
-        
+        ONLY:  m_soil_h,                                                       &
+               nzb_soil,                                                       &
+               nzt_soil,                                                       &
+               t_soil_h,                                                       &
+               zs
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  rad_lw_in,                                                      &
+               rad_lw_out,                                                     &
+               rad_sw_in,                                                      &
+               rad_sw_in_diff,                                                 &
+               rad_sw_out,                                                     &
+               radiation_scheme
+
     USE urban_surface_mod,                                                     &
-        ONLY:  nzb_wall, nzt_wall, t_wall_h 
+        ONLY:  nzb_wall,                                                       &
+               nzt_wall,                                                       &
+               t_wall_h
 
 
     IMPLICIT NONE
-    
+
     TYPE virt_general
-       INTEGER(iwp) ::  id_vm     !< NetCDF file id for virtual measurements
        INTEGER(iwp) ::  nvm = 0   !< number of virtual measurements
     END TYPE virt_general
 
+    TYPE virt_var_atts
+       CHARACTER(LEN=100) ::  coordinates          !< defined longname of the variable
+       CHARACTER(LEN=100) ::  grid_mapping         !< defined longname of the variable
+       CHARACTER(LEN=100) ::  long_name            !< defined longname of the variable
+       CHARACTER(LEN=100) ::  name                 !< variable name
+       CHARACTER(LEN=100) ::  standard_name        !< defined standard name of the variable
+       CHARACTER(LEN=100) ::  units                !< unit of the output variable
+
+       REAL(wp)           ::  fill_value = -9999.0 !< _FillValue attribute
+    END TYPE virt_var_atts
+
     TYPE virt_mea
-    
-       CHARACTER(LEN=100)  ::  feature_type      !< type of the measurement
-       CHARACTER(LEN=100)  ::  filename_original !< name of the original file
-       CHARACTER(LEN=100)  ::  site              !< name of the measurement site
-   
-       CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE ::  measured_vars_name !< name of the measured variables
-    
-       INTEGER(iwp) ::  ns = 0          !< number of observation coordinates on subdomain, for atmospheric measurements
-       INTEGER(iwp) ::  ns_tot = 0      !< total number of observation coordinates, for atmospheric measurements
-       INTEGER(iwp) ::  ntraj           !< number of trajectories of a measurement
-       INTEGER(iwp) ::  nmeas           !< number of measured variables (atmosphere + soil)
-       
-       INTEGER(iwp) ::  ns_soil = 0     !< number of observation coordinates on subdomain, for soil measurements
-       INTEGER(iwp) ::  ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
-       
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
-       
+
+       CHARACTER(LEN=100)  ::  feature_type                     !< type of the real-world measurement
+       CHARACTER(LEN=100)  ::  feature_type_out = 'timeSeries'  !< type of the virtual measurement
+                                                                !< (all will be timeSeries, even trajectories)
+       CHARACTER(LEN=100)  ::  nc_filename                      !< name of the NetCDF output file for the station
+       CHARACTER(LEN=100)  ::  site                             !< name of the measurement site
+
+       CHARACTER(LEN=1000) ::  data_content = REPEAT(' ', 1000) !< string of measured variables (data output only)
+
+       INTEGER(iwp) ::  end_coord_a = 0     !< end coordinate in NetCDF file for local atmosphere observations
+       INTEGER(iwp) ::  end_coord_s = 0     !< end coordinate in NetCDF file for local soil observations
+       INTEGER(iwp) ::  file_time_index = 0 !< time index in NetCDF output file
+       INTEGER(iwp) ::  ns = 0              !< number of observation coordinates on subdomain, for atmospheric measurements
+       INTEGER(iwp) ::  ns_tot = 0          !< total number of observation coordinates, for atmospheric measurements
+       INTEGER(iwp) ::  n_tr_st             !< number of trajectories / station of a measurement
+       INTEGER(iwp) ::  nmeas               !< number of measured variables (atmosphere + soil)
+       INTEGER(iwp) ::  ns_soil = 0         !< number of observation coordinates on subdomain, for soil measurements
+       INTEGER(iwp) ::  ns_soil_tot = 0     !< total number of observation coordinates, for soil measurements
+       INTEGER(iwp) ::  start_coord_a = 0   !< start coordinate in NetCDF file for local atmosphere observations
+       INTEGER(iwp) ::  start_coord_s = 0   !< start coordinate in NetCDF file for local soil observations
+
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory 
+                                                         !< or station that are no _FillValues
+
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< grid index for measurement position in x-direction
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< grid index for measurement position in y-direction
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< grid index for measurement position in k-direction
-       
+
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i_soil  !< grid index for measurement position in x-direction
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j_soil  !< grid index for measurement position in y-direction
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_soil  !< grid index for measurement position in k-direction
-           
+
        LOGICAL ::  trajectory         = .FALSE. !< flag indicating that the observation is a mobile observation
        LOGICAL ::  timseries          = .FALSE. !< flag indicating that the observation is a stationary point measurement
        LOGICAL ::  timseries_profile  = .FALSE. !< flag indicating that the observation is a stationary profile measurement
        LOGICAL ::  soil_sampling      = .FALSE. !< flag indicating that soil state variables were sampled
-       
-       REAL(wp) ::  fill_eutm          !< fill value for UTM coordinates in case of missing values
-       REAL(wp) ::  fill_nutm          !< fill value for UTM coordinates in case of missing values
-       REAL(wp) ::  fill_zag           !< fill value for heigth coordinates in case of missing values 
-       REAL(wp) ::  fillout = -999.9   !< fill value for output in case a observation is taken from inside a building 
-       REAL(wp) ::  origin_x_obs       !< origin of the observation in UTM coordiates in x-direction
-       REAL(wp) ::  origin_y_obs       !< origin of the observation in UTM coordiates in y-direction
-       
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  z_ag           !< measurement height above ground level
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth          !< measurement depth in soil
-             
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables 
-       
+
+       REAL(wp) ::  fill_eutm                            !< fill value for UTM coordinates in case of missing values
+       REAL(wp) ::  fill_nutm                            !< fill value for UTM coordinates in case of missing values
+       REAL(wp) ::  fill_zar                             !< fill value for heigth coordinates in case of missing values
+       REAL(wp) ::  fillout = -9999.0                    !< fill value for output in case a observation is taken 
+                                                         !< e.g. from inside a building
+       REAL(wp) ::  origin_x_obs                         !< origin of the observation in UTM coordiates in x-direction
+       REAL(wp) ::  origin_y_obs                         !< origin of the observation in UTM coordiates in y-direction
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zar           !< measurement height above ground level
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth         !< measurement depth in soil
+
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars       !< measured variables
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables
+
+       TYPE( virt_var_atts ), DIMENSION(:), ALLOCATABLE ::  var_atts !< variable attributes
+
     END TYPE virt_mea
 
     CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
     CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
-    
-    ! This need to generalized
-    CHARACTER(LEN=8)  ::  char_filename = "filename"               !< attribute name for filename
+
+    ! This need to be generalized
+!     CHARACTER(LEN=10) ::  char_fill = '_FillValue'
+    CHARACTER(LEN=9)  ::  char_long = 'long_name'                  !< attribute name for long_name
+    CHARACTER(LEN=13) ::  char_standard = 'standard_name'          !< attribute name for standard_name
+    CHARACTER(LEN=5)  ::  char_unit = 'units'                      !< attribute name for standard_name
     CHARACTER(LEN=11) ::  char_soil = "soil_sample"                !< attribute name for soil sampling indication
-    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue 
     CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
     CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
@@ -216,9 +316,10 @@
     CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
     CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
-    CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
+    CHARACTER(LEN=9)  ::  char_station_h = "station_h"             !< variable name indicating height of the site
+    CHARACTER(LEN=1)  ::  char_zar = "z"                           !< attribute name indicating height above reference level
     CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
     CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
     CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
-    
+
     CHARACTER(LEN=6), DIMENSION(1:5) ::  soil_vars       = (/                  & !< list of soil variables
                             't_soil',                                          &
@@ -227,106 +328,57 @@
                             'lwcs  ',                                          &
                             'smp   '                       /)
-                            
+
     CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
                                               'mcpm1     ', 'PM1       ',      &
                                               'mcpm2p5   ', 'PM2.5     ',      &
-                                              'mcpm25    ', 'PM25      ',      &
                                               'mcpm10    ', 'PM10      ',      &
                                               'mfno2     ', 'NO2       ',      &
                                               'mfno      ', 'NO        ',      &
-                                              'tro3      ', 'O3        ',      &
-                                              'mfco      ', 'CO        '       &
+                                              'mcno2     ', 'NO2       ',      &
+                                              'mcno      ', 'NO        ',      &
+                                              'tro3      ', 'O3        '       &
                                                                    /), (/ 2, 8 /) )
-!
-!-- MS: List requires careful revision!
-    CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
-       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
-          'hfss      ',  & ! surface sensible heat flux (W/m2)
-          'hur       ',  & ! relative humidity (-)
-          'hus       ',  & ! specific humidity (g/kg)
-          'haa       ',  & ! absolute atmospheric humidity (kg/m3)
-          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
-          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
-          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
-          'mcco      ',  & ! mass concentration of CO (kg/m3)
-          'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
-          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
-          'ncaa      ',  & ! number concentation of particles (1/m3)
-          'mfco      ',  & ! mole fraction of CO (mol/mol)
-          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
-          'mfch4     ',  & ! mole fraction of methane (mol/mol)
-          'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
-          'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
-          'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
-          'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
-          'mfh20     ',  & ! mole fraction of water (mol/mol)
-          'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
-          'rlds      ',  & ! surface downward longwave flux  (W/m2)
-          'rlus      ',  & ! surface upward longwave flux (W/m2)
-          'rsds      ',  & ! surface downward shortwave flux (W/m2)
-          'rsus      ',  & ! surface upward shortwave flux (W/m2)
-          'ta        ',  & ! air temperature (degree C)
-          't_va      ',  & ! virtual accoustic temperature (K)
-          'theta     ',  & ! potential temperature (K)
-          'tro3      ',  & ! mole fraction of ozone air (mol/mol)
-          'ts        ',  & ! scaling parameter of temperature (K)
-          'wspeed    ',  & ! ? wind speed - horizontal?
-          'wdir      ',  & ! wind direction
-          'us        ',  & ! friction velocity
-          'msoil     ',  & ! ? soil moisture - which depth?  
-          'tsoil     ',  & ! ? soil temperature - which depth?                                                                
-          'u         ',  & ! u-component
-          'utheta    ',  & ! total eastward kinematic heat flux
-          'ua        ',  & ! eastward wind (is there any difference to u?)
-          'v         ',  & ! v-component
-          'vtheta    ',  & ! total northward kinematic heat flux
-          'va        ',  & ! northward wind (is there any difference to v?)
-          'w         ',  & ! w-component
-          'wtheta    ',  & ! total vertical kinematic heat flux
-          'rld       ',  & ! downward longwave radiative flux (W/m2)
-          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
-          'rsd       ',  & ! downward shortwave radiative flux (W/m2)
-          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
-          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
-          'rnds      ',  & ! surface net downward radiative flux (W/m2)
-          't_soil    ',  &
-          'm_soil    ',  &
-          'lwc       ',  &
-          'lwcs      ',  &
-          'smp       '   &
-       /)
-                                                            
-    
-    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
-    LOGICAL ::  init = .TRUE.                     !< flag indicating initialization of data output
-    LOGICAL ::  use_virtual_measurement = .FALSE. !< Namelist parameter
-
-    REAL(wp) ::  dt_virtual_measurement = 0.0_wp    !< namelist parameter
-    REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last 3d output
-    REAL(wp) ::  vm_time_start = 0.0                !< time after virtual measurements should start
-
-    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass general variables
-    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
-    
+
+    LOGICAL ::  global_attribute = .TRUE.           !< flag indicating a global attribute
+    LOGICAL ::  initial_write_coordinates = .FALSE. !< flag indicating a global attribute
+    LOGICAL ::  use_virtual_measurement = .FALSE.   !< Namelist parameter
+
+    INTEGER(iwp) ::  maximum_name_length = 32 !< maximum name length of station names
+    INTEGER(iwp) ::  ntimesteps               !< number of timesteps defined in NetCDF output file
+    INTEGER(iwp) ::  ntimesteps_max = 80000   !< number of maximum timesteps defined in NetCDF output file
+    INTEGER(iwp) ::  off_pr = 1               !< number neighboring grid points (in each direction) where virtual profile
+                                              !< measurements shall be taken, in addition to the given coordinates in the driver
+    INTEGER(iwp) ::  off_ts = 1               !< number neighboring grid points (in each direction) where virtual timeseries
+                                              !< measurements shall be taken, in addition to the given coordinates in the driver
+    INTEGER(iwp) ::  off_tr = 1               !< number neighboring grid points (in each direction) where virtual trajectory
+                                              !< measurements shall be taken, in addition to the given coordinates in the driver
+
+    REAL(wp) ::  dt_virtual_measurement = 0.0_wp    !< sampling interval
+    REAL(wp) ::  time_virtual_measurement = 0.0_wp  !< time since last sampling
+    REAL(wp) ::  vm_time_start = 0.0                !< time after which sampling shall start
+
+    TYPE( virt_general )                        ::  vmea_general !< data structure which encompass global variables
+    TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE ::  vmea         !< data structure contain station-specific variables
+
     INTERFACE vm_check_parameters
        MODULE PROCEDURE vm_check_parameters
     END INTERFACE vm_check_parameters
-    
+
     INTERFACE vm_data_output
        MODULE PROCEDURE vm_data_output
     END INTERFACE vm_data_output
-    
+
     INTERFACE vm_init
        MODULE PROCEDURE vm_init
     END INTERFACE vm_init
-    
-    INTERFACE vm_last_actions
-       MODULE PROCEDURE vm_last_actions
-    END INTERFACE vm_last_actions
-    
+
+    INTERFACE vm_init_output
+       MODULE PROCEDURE vm_init_output
+    END INTERFACE vm_init_output
+
     INTERFACE vm_parin
        MODULE PROCEDURE vm_parin
     END INTERFACE vm_parin
-    
+
     INTERFACE vm_sampling
        MODULE PROCEDURE vm_sampling
@@ -339,10 +391,18 @@
 !
 !-- Public interfaces
-    PUBLIC  vm_check_parameters, vm_data_output, vm_init, vm_last_actions,     &
-            vm_parin, vm_sampling
+    PUBLIC  vm_check_parameters,                                               &
+            vm_data_output,                                                    &
+            vm_init,                                                           &
+            vm_init_output,                                                    &
+            vm_parin,                                                          &
+            vm_sampling
 
 !
 !-- Public variables
-    PUBLIC  dt_virtual_measurement, time_virtual_measurement, vmea, vmea_general, vm_time_start
+    PUBLIC  dt_virtual_measurement,                                            &
+            time_virtual_measurement,                                          &
+            vmea,                                                              &
+            vmea_general,                                                      &
+            vm_time_start
 
  CONTAINS
@@ -356,32 +416,370 @@
  SUBROUTINE vm_check_parameters
 
-    USE control_parameters,                                                    &
-        ONLY:  message_string, virtual_measurement
- 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  input_pids_static, input_pids_vm
-        
-    IMPLICIT NONE
-
-!
-!-- Virtual measurements require a setup file. 
-    IF ( virtual_measurement  .AND.  .NOT. input_pids_vm )  THEN
+    IF ( .NOT. virtual_measurement )  RETURN
+!
+!-- Virtual measurements require a setup file.
+    IF ( .NOT. input_pids_vm )  THEN
        message_string = 'If virtual measurements are taken, a setup input ' // &
                         'file for the site locations is mandatory.'
        CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
-    ENDIF   
+    ENDIF
 !
 !-- In case virtual measurements are taken, a static input file is required.
 !-- This is because UTM coordinates for the PALM domain origin are required
-!-- for correct mapping of the measurements. 
+!-- for correct mapping of the measurements.
 !-- ToDo: Revise this later and remove this requirement.
-    IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
+    IF ( .NOT. input_pids_static )  THEN
        message_string = 'If virtual measurements are taken, a static input ' //&
                         'file is mandatory.'
        CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
     ENDIF
- 
+
+#if !defined( __netcdf4_parallel )
+!
+!-- In case of non-parallel NetCDF the virtual measurement output is not
+!-- working. This is only designed for parallel NetCDF.
+    message_string = 'If virtual measurements are taken, parallel ' //         &
+                     'NetCDF is required.'
+    CALL message( 'vm_check_parameters', 'PA0708', 1, 2, 0, 6, 0 )
+#endif
+!
+!-- Check if the given number of neighboring grid points do not exceeds the number
+!-- of ghost points.
+    IF ( off_pr > nbgp - 1  .OR.  off_ts > nbgp - 1  .OR.  off_tr > nbgp - 1 ) &
+    THEN
+       WRITE(message_string,*)                                                 &
+                        'If virtual measurements are taken, the number ' //    &
+                        'of surrounding grid points must not be larger ' //    &
+                        'than the number of ghost points - 1, which is: ',     &
+                        nbgp - 1
+       CALL message( 'vm_check_parameters', 'PA0705', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Check if restart activation string is set. Virtual measurements may exceed
+!-- the maximum number of allowed timesteps in a NetCDF file. In order to handle
+!-- this, the module will trigger a programm abortion with writing the restart
+!-- data.
+    IF ( .NOT. write_binary )  THEN
+       message_string = 'virtual measurements require file activation ' //     &
+                        'string "restart"'
+       CALL message( 'check_parameters', 'PA0706', 1, 2, 0, 6, 0 )
+    ENDIF
+
  END SUBROUTINE vm_check_parameters
  
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine defines variable attributes according to UC2 standard. Note, later
+!> this list can be moved to the data-output module where it can be re-used also
+!> for other output.
+!------------------------------------------------------------------------------!
+  SUBROUTINE vm_set_attributes( output_variable )
+
+     TYPE( virt_var_atts ), INTENT(INOUT) ::  output_variable !< data structure with attributes that need to be set
+
+     output_variable%long_name     = 'none'
+     output_variable%standard_name = 'none'
+     output_variable%units         = 'none'
+     output_variable%coordinates   = 'lon lat E_UTM N_UTM x y z time station_name'
+     output_variable%grid_mapping  = 'crs'
+
+     SELECT CASE ( TRIM( output_variable%name ) )
+
+        CASE ( 'u' )
+           output_variable%long_name     = 'u wind component'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'ua' )
+           output_variable%long_name     = 'eastward wind'
+           output_variable%standard_name = 'eastward_wind'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'v' )
+           output_variable%long_name     = 'v wind component'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'va' )
+           output_variable%long_name     = 'northward wind'
+           output_variable%standard_name = 'northward_wind'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'w' )
+           output_variable%long_name     = 'w wind component'
+           output_variable%standard_name = 'upward_air_velocity'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'wspeed' )
+           output_variable%long_name     = 'wind speed'
+           output_variable%standard_name = 'wind_speed'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'wdir' )
+           output_variable%long_name     = 'wind from direction'
+           output_variable%standard_name = 'wind_from_direction'
+           output_variable%units         = 'degrees'
+
+        CASE ( 'theta' )
+           output_variable%long_name     = 'air potential temperature'
+           output_variable%standard_name = 'air_potential_temperature'
+           output_variable%units         = 'K'
+
+        CASE ( 'utheta' )
+           output_variable%long_name     = 'eastward kinematic sensible heat flux in air'
+           output_variable%units         = 'K m s-1'
+
+        CASE ( 'vtheta' )
+           output_variable%long_name     = 'northward kinematic sensible heat flux in air'
+           output_variable%units         = 'K m s-1'
+
+        CASE ( 'wtheta' )
+           output_variable%long_name     = 'upward kinematic sensible heat flux in air'
+           output_variable%units         = 'K m s-1'
+
+        CASE ( 'ta' )
+           output_variable%long_name     = 'air temperature'
+           output_variable%standard_name = 'air_temperature'
+           output_variable%units         = 'degree_C'
+
+        CASE ( 'tva' )
+           output_variable%long_name     = 'virtual acoustic temperature'
+           output_variable%units         = 'K'
+
+        CASE ( 'haa' )
+           output_variable%long_name     = 'absolute atmospheric humidity'
+           output_variable%units         = 'kg m-3'
+
+        CASE ( 'hus' )
+           output_variable%long_name     = 'specific humidity'
+           output_variable%standard_name = 'specific_humidity'
+           output_variable%units         = 'kg kg-1'
+
+        CASE ( 'hur' )
+           output_variable%long_name     = 'relative humidity'
+           output_variable%standard_name = 'relative_humidity'
+           output_variable%units         = '1'
+
+        CASE ( 'rlu' )
+           output_variable%long_name     = 'upwelling longwave flux in air'
+           output_variable%standard_name = 'upwelling_longwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rlus' )
+           output_variable%long_name     = 'surface upwelling longwave flux in air'
+           output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rld' )
+           output_variable%long_name     = 'downwelling longwave flux in air'
+           output_variable%standard_name = 'downwelling_longwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rsddif' )
+           output_variable%long_name     = 'diffuse downwelling shortwave flux in air'
+           output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rsd' )
+           output_variable%long_name     = 'downwelling shortwave flux in air'
+           output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rnds' )
+           output_variable%long_name     = 'surface net downward radiative flux'
+           output_variable%standard_name = 'surface_net_downward_radiative_flux'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rsu' )
+           output_variable%long_name     = 'upwelling shortwave flux in air'
+           output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rsus' )
+           output_variable%long_name     = 'surface upwelling shortwave flux in air'
+           output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'rsds' )
+           output_variable%long_name     = 'surface downwelling shortwave flux in air'
+           output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'hfss' )
+           output_variable%long_name     = 'surface upward sensible heat flux'
+           output_variable%standard_name = 'surface_upward_sensible_heat_flux'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'hfls' )
+           output_variable%long_name     = 'surface upward latent heat flux'
+           output_variable%standard_name = 'surface_upward_latent_heat_flux'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'ts' )
+           output_variable%long_name     = 'surface temperature'
+           output_variable%standard_name = 'surface_temperature'
+           output_variable%units         = 'K'
+
+        CASE ( 'thetas' )
+           output_variable%long_name     = 'surface layer temperature scale'
+           output_variable%units         = 'K'
+
+        CASE ( 'us' )
+           output_variable%long_name     = 'friction velocity'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'uw' )
+           output_variable%long_name     = 'upward eastward kinematic momentum flux in air'
+           output_variable%units         = 'm2 s-2'
+
+        CASE ( 'vw' )
+           output_variable%long_name     = 'upward northward kinematic momentum flux in air'
+           output_variable%units         = 'm2 s-2'
+
+        CASE ( 'uv' )
+           output_variable%long_name     = 'eastward northward kinematic momentum flux in air'
+           output_variable%units         = 'm2 s-2'
+
+        CASE ( 'plev' )
+           output_variable%long_name     = 'air pressure'
+           output_variable%standard_name = 'air_pressure'
+           output_variable%units         = 'Pa'
+
+        CASE ( 'm_soil' )
+           output_variable%long_name     = 'soil moisture volumetric'
+           output_variable%units         = 'm3 m-3'
+
+        CASE ( 't_soil' )
+           output_variable%long_name     = 'soil temperature'
+           output_variable%standard_name = 'soil_temperature'
+           output_variable%units         = 'degree_C'
+
+        CASE ( 'hfdg' )
+           output_variable%long_name     = 'downward heat flux at ground level in soil'
+           output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'hfds' )
+           output_variable%long_name     = 'downward heat flux in soil'
+           output_variable%standard_name = 'downward_heat_flux_in_soil'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'hfla' )
+           output_variable%long_name     = 'upward latent heat flux in air'
+           output_variable%standard_name = 'upward_latent_heat_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'hfsa' )
+           output_variable%long_name     = 'upward latent heat flux in air'
+           output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
+           output_variable%units         = 'W m-2'
+
+        CASE ( 'jno2' )
+           output_variable%long_name     = 'photolysis rate of nitrogen dioxide'
+           output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
+           output_variable%units         = 's-1'
+
+        CASE ( 'lwcs' )
+           output_variable%long_name     = 'liquid water content of soil layer'
+           output_variable%standard_name = 'liquid_water_content_of_soil_layer'
+           output_variable%units         = 'kg m-2'
+
+        CASE ( 'lwp' )
+           output_variable%long_name     = 'liquid water path'
+           output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
+           output_variable%units         = 'kg m-2'
+
+        CASE ( 'ps' )
+           output_variable%long_name     = 'surface air pressure'
+           output_variable%standard_name = 'surface_air_pressure'
+           output_variable%units         = 'hPa'
+
+        CASE ( 'pswrtg' )
+           output_variable%long_name     = 'platform speed wrt ground'
+           output_variable%standard_name = 'platform_speed_wrt_ground'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'pswrta' )
+           output_variable%long_name     = 'platform speed wrt air'
+           output_variable%standard_name = 'platform_speed_wrt_air'
+           output_variable%units         = 'm s-1'
+
+        CASE ( 'pwv' )
+           output_variable%long_name     = 'water vapor partial pressure in air'
+           output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
+           output_variable%units         = 'hPa'
+
+        CASE ( 'ssdu' )
+           output_variable%long_name     = 'duration of sunshine'
+           output_variable%standard_name = 'duration_of_sunshine'
+           output_variable%units         = 's'
+
+        CASE ( 't_lw' )
+           output_variable%long_name     = 'land water temperature'
+           output_variable%units         = 'degree_C'
+
+        CASE ( 'tb' )
+           output_variable%long_name     = 'brightness temperature'
+           output_variable%standard_name = 'brightness_temperature'
+           output_variable%units         = 'K'
+
+        CASE ( 'uqv' )
+           output_variable%long_name     = 'eastward kinematic latent heat flux in air'
+           output_variable%units         = 'g kg-1 m s-1'
+
+        CASE ( 'vqv' )
+           output_variable%long_name     = 'northward kinematic latent heat flux in air'
+           output_variable%units         = 'g kg-1 m s-1'
+
+        CASE ( 'wqv' )
+           output_variable%long_name     = 'upward kinematic latent heat flux in air'
+           output_variable%units         = 'g kg-1 m s-1'
+
+        CASE ( 'zcb' )
+           output_variable%long_name     = 'cloud base altitude'
+           output_variable%standard_name = 'cloud_base_altitude'
+           output_variable%units         = 'm'
+
+        CASE ( 'zmla' )
+           output_variable%long_name     = 'atmosphere boundary layer thickness'
+           output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
+           output_variable%units         = 'm'
+
+        CASE ( 'mcpm1' )
+           output_variable%long_name     = 'mass concentration of pm1 ambient aerosol particles in air'
+           output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
+           output_variable%units         = 'kg m-3'
+
+        CASE ( 'mcpm10' )
+           output_variable%long_name     = 'mass concentration of pm10 ambient aerosol particles in air'
+           output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
+           output_variable%units         = 'kg m-3'
+
+        CASE ( 'mcpm2p5' )
+           output_variable%long_name     = 'mass concentration of pm2p5 ambient aerosol particles in air'
+           output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
+           output_variable%units         = 'kg m-3'
+
+        CASE ( 'mfno', 'mcno'  )
+           output_variable%long_name     = 'mole fraction of nitrogen monoxide in air'
+           output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
+           output_variable%units         = 'ppm' !'mol mol-1'
+
+        CASE ( 'mfno2', 'mcno2'  )
+           output_variable%long_name     = 'mole fraction of nitrogen dioxide in air'
+           output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
+           output_variable%units         = 'ppm' !'mol mol-1'
+
+        CASE ( 'tro3'  )
+           output_variable%long_name     = 'mole fraction of ozone in air'
+           output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
+           output_variable%units         = 'ppm' !'mol mol-1'
+
+        CASE DEFAULT
+
+     END SELECT
+
+  END SUBROUTINE vm_set_attributes
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -390,13 +788,15 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE vm_parin
- 
+
     CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
- 
-    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,                             &
-                                               use_virtual_measurement,                            &
+
+    NAMELIST /virtual_measurement_parameters/  dt_virtual_measurement,         &
+                                               off_ts,                         &
+                                               off_pr,                         &
+                                               off_tr,                         &
+                                               use_virtual_measurement,        &
                                                vm_time_start
 
     line = ' '
-
 !
 !-- Try to find stg package
@@ -415,5 +815,5 @@
 !-- Set flag that indicates that the virtual measurement module is switched on
     IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
-    
+
     GOTO 20
 
@@ -423,5 +823,5 @@
 
  20 CONTINUE
- 
+
  END SUBROUTINE vm_parin
 
@@ -430,5 +830,5 @@
 ! Description:
 ! ------------
-!> Initialize virtual measurements: read coordiante arrays and measured 
+!> Initialize virtual measurements: read coordiante arrays and measured
 !> variables, set indicies indicating the measurement points, read further
 !> attributes, etc..
@@ -436,27 +836,9 @@
  SUBROUTINE vm_init
 
-    USE arrays_3d,                                                             &
-        ONLY:  zu, zw
-        
-    USE grid_variables,                                                        &
-        ONLY:  ddx, ddy, dx, dy
-        
-    USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys
- 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  get_dimension_length,                                           &
-               init_model,                                                     &
-               input_file_vm,                                                  &
-               netcdf_data_input_att,                                          &
-               netcdf_data_input_var
-        
-    IMPLICIT NONE
-    
-    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
-    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
-    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed   
-    
-    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate 
+    CHARACTER(LEN=5)                  ::  dum                          !< dummy string indicating station id
+    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
+    CHARACTER(LEN=100), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
+
+    INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate
     INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
     INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
@@ -471,52 +853,60 @@
     INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
     INTEGER(iwp) ::  ll        !< running index over all measured variables in file
-    INTEGER(iwp) ::  lll       !< running index over all allowed variables
+    INTEGER(iwp) ::  m         !< running index for surface elements
     INTEGER(iwp) ::  n         !< running index over trajectory coordinates
+    INTEGER(iwp) ::  nofill    !< dummy for nofill return value (not used)
     INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
+    INTEGER(iwp) ::  off       !< number of surrounding grid points to be sampled
     INTEGER(iwp) ::  t         !< running index over number of trajectories
-    INTEGER(iwp) ::  m
-    
-    INTEGER(KIND=1)::  soil_dum
-    
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_all !< dummy array used to sum-up the number of observation coordinates
-    
+
+    INTEGER(KIND=1)                             ::  soil_dum !< dummy variable to input a soil flag
+
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE     ::  ns_all !< dummy array used to sum-up the number of observation coordinates
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_atmos !< number of observation points for each station on each mpi rank
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  ns_soil  !< number of observation points for each station on each mpi rank
+
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
-    
-!    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
-    LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
-    
-    REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
-    REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
-    REAL(wp)     ::  fill_zag  !< _FillValue for height coordinate
-    
+
+    LOGICAL  ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
+
+    REAL(wp) ::  fill_eutm !< _FillValue for coordinate array E_UTM
+    REAL(wp) ::  fill_nutm !< _FillValue for coordinate array N_UTM
+    REAL(wp) ::  fill_zar  !< _FillValue for height coordinate
+
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm     !< easting UTM coordinate, temporary variable
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm     !< northing UTM coordinate, temporary variable
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm_tmp !< EUTM coordinate before rotation
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm_tmp !< NUTM coordinate before rotation
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag      !< height coordinate relative to origin_z, temporary variable
-!
-!-- Obtain number of sites. Also, pass the 'open' string, in order to initially
-!-- open the measurement driver. 
-    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
-                                vmea_general%id_vm, input_file_vm,             &
-                                global_attribute, 'open', '' )
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  station_h !< station height above reference
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zar       !< observation height above reference
+#if defined( __netcdf )
+!
+!-- Open the input file.
+    CALL open_read_file( input_file_vm, pids_id )
+!
+!-- Obtain number of sites.
+    CALL get_attribute( pids_id,                                               &
+                        char_numstations,                                      &
+                        vmea_general%nvm,                                      &
+                        global_attribute )
 !
 !-- Allocate data structure which encompass all required information, such as
-!-- grid points indicies, absolute UTM coordinates, the measured quantities, 
+!-- grid points indicies, absolute UTM coordinates, the measured quantities,
 !-- etc. .
     ALLOCATE( vmea(1:vmea_general%nvm) )
 !
 !-- Allocate flag array. This dummy array is used to identify grid points
-!-- where virtual measurements should be taken. Please note, at least one 
-!-- ghost point is required, in order to include also the surrounding 
-!-- grid points of the original coordinate.  
-    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+!-- where virtual measurements should be taken. Please note, in order to 
+!-- include also the surrounding grid points of the original coordinate
+!-- ghost points are required.
+    ALLOCATE( meas_flag(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     meas_flag = 0
 !
-!-- Loop over all sites.
+!-- Loop over all sites in the setup file.
     DO  l = 1, vmea_general%nvm
 !
 !--    Determine suffix which contains the ID, ordered according to the number
-!--    of measurements. 
+!--    of measurements.
        IF( l < 10 )  THEN
           WRITE( dum, '(I1)')  l
@@ -531,35 +921,35 @@
        ENDIF
 !
-!--    Read site coordinates (UTM).
-       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx //         &
-                                   TRIM( dum ), vmea_general%id_vm, '',        &
-                                   global_attribute, '', '' )
-       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy //         &
-                                   TRIM( dum ), vmea_general%id_vm, '',        &
-                                   global_attribute, '', '' )
-!
-!--    Read site name                  
-       CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ),     &
-                                   vmea_general%id_vm, '', global_attribute,   &
-                                   '', '' )
+!--    Read the origin site coordinates (UTM).
+       CALL get_attribute( pids_id,                                            &
+                           char_origx // TRIM( dum ),                          &
+                           vmea(l)%origin_x_obs,                               &
+                           global_attribute )
+       CALL get_attribute( pids_id,                                            &
+                           char_origy // TRIM( dum ),                          &
+                           vmea(l)%origin_y_obs,                               &
+                           global_attribute )
+!
+!--    Read site name.
+       CALL get_attribute( pids_id,                                            &
+                           char_site // TRIM( dum ),                           &
+                           vmea(l)%site,                                       &
+                           global_attribute )
+!
+!--    Read a flag which indicates that also soil quantities are take at the
+!--    respective site (is part of the virtual measurement driver).
+       CALL get_attribute( pids_id,                                            &
+                           char_soil // TRIM( dum ),                           &
+                           soil_dum,                                           &
+                           global_attribute )
+!
+!--    Set flag indicating soil-sampling.
+       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
 !
 !--    Read type of the measurement (trajectory, profile, timeseries).
-       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature //       &
-                                   TRIM( dum ), vmea_general%id_vm, '',        &
-                                   global_attribute, '', '' )
-!
-!--    Read the name of the original file where observational data is stored. 
-       CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
-                                   TRIM( dum ), vmea_general%id_vm, '',        &
-                                   global_attribute, '', '' )
-!
-!--    Read a flag which indicates that also soil quantities are take at the 
-!--    respective site (is part of the virtual measurement driver).  
-       CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ),         &
-                                   vmea_general%id_vm, '', global_attribute,   &
-                                   '', '' )
-!
-!--    Set flag for soil-sampling.
-       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
+       CALL get_attribute( pids_id,                                            &
+                           char_feature // TRIM( dum ),                        &
+                           vmea(l)%feature_type,                               &
+                           global_attribute )
 !
 !---   Set logicals depending on the type of the measurement
@@ -571,24 +961,24 @@
           vmea(l)%trajectory        = .TRUE.
 !
-!--   Give error message in case the type matches non of the pre-defined types.
+!--    Give error message in case the type matches non of the pre-defined types.
        ELSE
           message_string = 'Attribue featureType = ' //                        &
                            TRIM( vmea(l)%feature_type ) //                     &
-                           ' is not allowed.' 
+                           ' is not allowed.'
           CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
        ENDIF
 !
-!--    Read string with all measured variables at this site
+!--    Read string with all measured variables at this site.
        measured_variables_file = ''
-       CALL netcdf_data_input_var( measured_variables_file,                    &
-                                   char_mv // TRIM( dum ), vmea_general%id_vm )
-!
-!--    Count the number of measured variables. Only count variables that match
-!--    with the allowed variables. 
+       CALL get_variable( pids_id,                                             &
+                          char_mv // TRIM( dum ),                              &
+                          measured_variables_file )
+!
+!--    Count the number of measured variables.
 !--    Please note, for some NetCDF interal reasons characters end with a NULL,
 !--    i.e. also empty characters contain a NULL. Therefore, check the strings
-!--    for a NULL to get the correct character length in order to compare 
-!--    them with the list of allowed variables. 
-       vmea(l)%nmeas  = 0
+!--    for a NULL to get the correct character length in order to compare
+!--    them with the list of allowed variables.
+       vmea(l)%nmeas  = 1
        DO ll = 1, SIZE( measured_variables_file )
           IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
@@ -602,58 +992,39 @@
              ENDDO
              len_char = len_char - 1
-!
-!--          Now, compare the measured variable with the list of allowed 
-!--          variables. 
-             DO  lll= 1, SIZE( list_allowed_variables )
-                IF ( measured_variables_file(ll)(1:len_char) ==                &
-                     TRIM( list_allowed_variables(lll) ) )  THEN
-                   vmea(l)%nmeas = vmea(l)%nmeas + 1
-                   measured_variables(vmea(l)%nmeas) =                         &
+
+             measured_variables(vmea(l)%nmeas) =                               &
                                        measured_variables_file(ll)(1:len_char)
-                ENDIF
-             ENDDO
+             vmea(l)%nmeas = vmea(l)%nmeas + 1
+
           ENDIF
        ENDDO
-       
-       
-!
-!--    Allocate array for the measured variables names for the respective site.
-       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nmeas) )
-
+       vmea(l)%nmeas = vmea(l)%nmeas - 1
+!
+!--    Allocate data-type array for the measured variables names and attributes
+!--    at the respective site.
+       ALLOCATE( vmea(l)%var_atts(1:vmea(l)%nmeas) )
+!
+!--    Store the variable names in a data structures, which assigns further 
+!--    attributes to this name. Further, for data output reasons, create a 
+!--    string of output variables, which will be written into the attribute
+!--    data_content.
        DO  ll = 1, vmea(l)%nmeas
-          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
+          vmea(l)%var_atts(ll)%name = TRIM( measured_variables(ll) )
+
+          vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " //        &
+                                 TRIM( vmea(l)%var_atts(ll)%name )
        ENDDO
 !
-!--    In case of chemistry, check if species is considered in the modelled
-!--    chemistry mechanism.
-!        IF ( air_chemistry )  THEN
-!           DO  ll = 1, vmea(l)%nmeas
-!              chem_include = .FALSE.
-!              DO  n = 1, nvar
-!                 IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
-!                      TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
-!              ENDDO
-! !
-! !--  Revise this. It should only check for chemistry variables and not for all!
-!              IF ( .NOT. chem_include )  THEN
-!                 message_string = TRIM( vmea(l)%measured_vars_name(ll) ) //     &
-!                                  ' is not considered in the modelled '  //     &
-!                                  'chemistry mechanism'
-!                 CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
-!              ENDIF
-!           ENDDO
-!        ENDIF
-!
-!--    Read the UTM coordinates for the actual site. Based on the coordinates,
+!--    Read all the UTM coordinates for the site. Based on the coordinates,
 !--    define the grid-index space on each subdomain where virtual measurements
-!--    should be taken. Note, the entire coordinate arrays will not be stored 
-!--    as this would exceed memory requirements, particularly for trajectory 
-!--    measurements. 
+!--    should be taken. Note, the entire coordinate array (on the entire model
+!--    domain) won't be stored as this would exceed memory requirements, 
+!--    particularly for trajectories.
        IF ( vmea(l)%nmeas > 0 )  THEN
 !
-!--       For stationary measurements UTM coordinates are just one value and 
-!--       its dimension is "station", while for mobile measurements UTM 
+!--       For stationary measurements UTM coordinates are just one value and
+!--       its dimension is "station", while for mobile measurements UTM
 !--       coordinates are arrays depending on the number of trajectories and
-!--       time, according to (UC)2 standard. First, inquire dimension length 
+!--       time, according to (UC)2 standard. First, inquire dimension length
 !--       of the UTM coordinates.
           IF ( vmea(l)%trajectory )  THEN
@@ -661,78 +1032,110 @@
 !--          For non-stationary measurements read the number of trajectories
 !--          and the number of time coordinates.
-             CALL get_dimension_length( vmea_general%id_vm,                    &
-                                        vmea(l)%ntraj,                         &
-                                        "traj" //                              &
-                                        TRIM( dum ) )
-             CALL get_dimension_length( vmea_general%id_vm,                    &
+             CALL get_dimension_length( pids_id,                               &
+                                        vmea(l)%n_tr_st,                       &
+                                        "traj" // TRIM( dum ) )
+             CALL get_dimension_length( pids_id,                               &
                                         dim_ntime,                             &
-                                        "ntime" //                             &
-                                        TRIM( dum ) )
-!
-!--       For stationary measurements the dimension for UTM and time 
-!--       coordinates is 1. 
+                                        "ntime" // TRIM( dum ) )
+!
+!--       For stationary measurements the dimension for UTM is station
+!--       and for the time-coordinate it is one.
           ELSE
-             vmea(l)%ntraj  = 1
+             CALL get_dimension_length( pids_id,                               &
+                                        vmea(l)%n_tr_st,                       &
+                                        "station" // TRIM( dum ) )
              dim_ntime = 1
           ENDIF
 !
-!-        Allocate array which defines individual time frame for each 
+!-        Allocate array which defines individual time/space frame for each
 !--       trajectory or station.
-          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
-!
-!--       Allocate temporary arrays for UTM and height coordinates. Note, 
+          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
+!
+!--       Allocate temporary arrays for UTM and height coordinates. Note,
 !--       on file UTM coordinates might be 1D or 2D variables
-          ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
-          ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
-          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
-          
-          ALLOCATE( e_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
-          ALLOCATE( n_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
-!
-!--       Read _FillValue attributes of the coordinate dimensions.
-          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
-                                      vmea_general%id_vm, '',                  &
-                                      .NOT. global_attribute, '',              &
-                                      char_eutm // TRIM( dum ) )
-          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
-                                      vmea_general%id_vm, '',                  &
-                                      .NOT. global_attribute, '',              &
-                                      char_nutm // TRIM( dum ) )
-          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
-                                      vmea_general%id_vm, '',                  &
-                                      .NOT. global_attribute, '',              &
-                                      char_zag  // TRIM( dum ) )
-!
-!--       Read UTM and height coordinates coordinates for all trajectories and 
-!--       times. 
+          ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
+          ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime)       )
+          ALLOCATE( station_h(1:vmea(l)%n_tr_st,1:dim_ntime)   )
+          ALLOCATE( zar(1:vmea(l)%n_tr_st,1:dim_ntime)         )
+          e_utm     = 0.0_wp
+          n_utm     = 0.0_wp
+          station_h = 0.0_wp
+          zar       = 0.0_wp
+
+          ALLOCATE( e_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
+          ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
+!
+!--       Read UTM and height coordinates coordinates for all trajectories and
+!--       times. Note, in case these obtain any missing values, replace them
+!--       with default _FillValues.
+          CALL inquire_fill_value( pids_id,                                    &
+                                   char_eutm // TRIM( dum ),                   &
+                                   nofill,                                     &
+                                   fill_eutm )
+          CALL inquire_fill_value( pids_id,                                    &
+                                   char_nutm // TRIM( dum ),                   &
+                                   nofill,                                     &
+                                   fill_nutm )
+          CALL inquire_fill_value( pids_id,                                    &
+                                   char_zar // TRIM( dum ),                    &
+                                   nofill,                                     &
+                                   fill_zar )
+!
+!--       Further line is just to avoid compiler warnings. nofill might be used
+!--       in future.
+          IF ( nofill == 0  .OR.  nofill /= 0 )  CONTINUE
+!
+!--       Read observation coordinates. Please note, for trajectories the
+!--       observation height is stored directly in z, while for timeSeries
+!--       it is stored in z - station_h, according to UC2-standard.
           IF ( vmea(l)%trajectory )  THEN
-             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ),      &
-                                         vmea_general%id_vm,                   &
-                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
-             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ),      &
-                                         vmea_general%id_vm,                   &
-                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
-             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ),        &
-                                         vmea_general%id_vm,                   &
-                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
+             CALL get_variable( pids_id,                                       &
+                                char_eutm // TRIM( dum ),                      &
+                                e_utm,                                         &
+                                0, dim_ntime-1,                                &
+                                0, vmea(l)%n_tr_st-1 )
+             CALL get_variable( pids_id,                                       &
+                                char_nutm // TRIM( dum ),                      &
+                                n_utm,                                         &
+                                0, dim_ntime-1,                                &
+                                0, vmea(l)%n_tr_st-1 )
+             CALL get_variable( pids_id,                                       &
+                                char_zar // TRIM( dum ),                       &
+                                zar,                                           &
+                                0, dim_ntime-1,                                &
+                                0, vmea(l)%n_tr_st-1 )
           ELSE
-             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
-                                         vmea_general%id_vm )
-             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
-                                         vmea_general%id_vm )
-             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), &
-                                         vmea_general%id_vm )
-          ENDIF         
+             CALL get_variable( pids_id,                                       &
+                                char_eutm // TRIM( dum ),                      &
+                                e_utm(:,1) )
+             CALL get_variable( pids_id,                                       &
+                                char_nutm // TRIM( dum ),                      &
+                                n_utm(:,1) )
+             CALL get_variable( pids_id,                                       &
+                                char_station_h // TRIM( dum ),                 &
+                                station_h(:,1) )
+             CALL get_variable( pids_id,                                       &
+                                char_zar // TRIM( dum ),                       &
+                                zar(:,1) )
+          ENDIF
+
+          e_utm = MERGE( e_utm, vmea(l)%fillout, e_utm /= fill_eutm )
+          n_utm = MERGE( n_utm, vmea(l)%fillout, n_utm /= fill_nutm )
+          zar   = MERGE( zar,   vmea(l)%fillout, zar   /= fill_zar  )
+!
+!--       Compute observation height above ground.
+          zar  = zar - station_h
 !
 !--       Based on UTM coordinates, check if the measurement station or parts
-!--       of the trajectory is on subdomain. This case, setup grid index space
-!--       sample these quantities. 
+!--       of the trajectory are on subdomain. This case, setup grid index space
+!--       sample these quantities.
           meas_flag = 0
-          DO  t = 1, vmea(l)%ntraj
-!             
-!--          First, compute relative x- and y-coordinates with respect to the 
-!--          lower-left origin of the model domain, which is the difference 
+          DO  t = 1, vmea(l)%n_tr_st
+!
+!--          First, compute relative x- and y-coordinates with respect to the
+!--          lower-left origin of the model domain, which is the difference
 !--          between UTM coordinates. Note, if the origin is not correct, the
-!--          virtual sites will be misplaced. 
+!--          virtual sites will be misplaced. Further, in case of an rotated
+!--          model domain, the UTM coordinates must be also rotated.
              e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
              n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
@@ -749,21 +1152,31 @@
 !--          trajectory. This is required as several stations and trajectories
 !--          are merged into one file but they do not have the same number of
-!--          points in time, hence, missing values may occur and cannot be 
+!--          points in time, hence, missing values may occur and cannot be
 !--          processed further. This is actually a work-around for the specific
-!--          (UC)2 dataset, but it won't harm in anyway. 
+!--          (UC)2 dataset, but it won't harm anyway.
              vmea(l)%dim_t(t) = 0
              DO  n = 1, dim_ntime
                 IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
                      n_utm(t,n) /= fill_nutm  .AND.                            &
-                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
+                     zar(t,n)   /= fill_zar )  vmea(l)%dim_t(t) = n
              ENDDO
 !
 !--          Compute grid indices relative to origin and check if these are
-!--          on the subdomain. Note, virtual measurements will be taken also 
-!--          at grid points surrounding the station, hence, check also for 
+!--          on the subdomain. Note, virtual measurements will be taken also
+!--          at grid points surrounding the station, hence, check also for
 !--          these grid points.
+!--          The number of surrounding grid points is set according to the
+!--          featureType.
+             IF ( vmea(l)%timseries_profile )  THEN
+                off = off_pr
+             ELSEIF ( vmea(l)%timseries     )  THEN
+                off = off_ts
+             ELSEIF ( vmea(l)%trajectory    )  THEN
+                off = off_tr
+             ENDIF
+
              DO  n = 1, vmea(l)%dim_t(t)
                 is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
-                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
+                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
 !
 !--             Is the observation point on subdomain?
@@ -777,19 +1190,19 @@
 !--                height.
                    ksurf = topo_top_ind(js,is,0)
-                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
+                   ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
 !
 !--                Set mask array at the observation coordinates. Also, flag the
 !--                surrounding coordinate points, but first check whether the
-!--                surrounding coordinate points are on the subdomain. 
-                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 >= ksurf )
-                   ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
-                 
-                   DO  i = is-1, is+1
-                      DO  j = js-1, js+1
+!--                surrounding coordinate points are on the subdomain.
+                   kl = MERGE( ks-off, ksurf, ks-off >= nzb  .AND. ks-off >= ksurf )
+                   ku = MERGE( ks+off, nzt,   ks+off < nzt+1 )
+
+                   DO  i = is-off, is+off
+                      DO  j = js-off, js+off
                          DO  k = kl, ku
-                            meas_flag(k,j,i) = MERGE(                              &
-                                             IBSET( meas_flag(k,j,i), 0 ),         &
-                                             0,                                    &
-                                             BTEST( wall_flags_total_0(k,j,i), 0 ) &
+                            meas_flag(k,j,i) = MERGE(                           &
+                                          IBSET( meas_flag(k,j,i), 0 ),         &
+                                          0,                                    &
+                                          BTEST( wall_flags_total_0(k,j,i), 0 ) &
                                                     )
                          ENDDO
@@ -798,14 +1211,14 @@
                 ENDIF
              ENDDO
-             
+
           ENDDO
 !
-!--       Based on the flag array count the number of of sampling coordinates.
-!--       Please note, sampling coordinates in atmosphere and soil may be 
-!--       different, as within the soil all levels will be measured.           
+!--       Based on the flag array count the number of sampling coordinates.
+!--       Please note, sampling coordinates in atmosphere and soil may be
+!--       different, as within the soil all levels will be measured.
 !--       Hence, count individually. Start with atmoshere.
           ns = 0
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxl-off, nxr+off
+             DO  j = nys-off, nyn+off
                 DO  k = nzb, nzt+1
                    ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
@@ -813,20 +1226,20 @@
              ENDDO
           ENDDO
-          
+
 !
 !--       Store number of observation points on subdomain and allocate index
 !--       arrays as well as array containing height information.
           vmea(l)%ns = ns
-          
+
           ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
           ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
           ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
-          ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )      
-!
-!--       Based on the flag array store the grid indices which correspond to 
-!--       the observation coordinates. 
+          ALLOCATE( vmea(l)%zar(1:vmea(l)%ns) )
+!
+!--       Based on the flag array store the grid indices which correspond to
+!--       the observation coordinates.
           ns = 0
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxl-off, nxr+off
+             DO  j = nys-off, nyn+off
                 DO  k = nzb, nzt+1
                    IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
@@ -835,5 +1248,5 @@
                       vmea(l)%j(ns) = j
                       vmea(l)%k(ns) = k
-                      vmea(l)%z_ag(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
+                      vmea(l)%zar(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
                    ENDIF
                 ENDDO
@@ -841,6 +1254,8 @@
           ENDDO
 !
-!--       Same for the soil. Based on the flag array, count the number of 
-!--       sampling coordinates in soil. Sample at all soil levels in this case. 
+!--       Same for the soil. Based on the flag array, count the number of
+!--       sampling coordinates in soil. Sample at all soil levels in this case.
+!--       Please note, soil variables can only be sampled on subdomains, not
+!--       on ghost layers.
           IF ( vmea(l)%soil_sampling )  THEN
              DO  i = nxl, nxr
@@ -850,17 +1265,17 @@
                            surf_lsm_h%end_index(j,i) )  THEN
                          vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
-                                                      nzt_soil - nzb_soil + 1 
+                                                      nzt_soil - nzb_soil + 1
                       ENDIF
                       IF ( surf_usm_h%start_index(j,i) <=                      &
                            surf_usm_h%end_index(j,i) )  THEN
                          vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
-                                                      nzt_wall - nzb_wall + 1 
+                                                      nzt_wall - nzb_wall + 1
                       ENDIF
                    ENDIF
                 ENDDO
              ENDDO
-          ENDIF         
-!
-!--       Allocate index arrays as well as array containing height information 
+          ENDIF
+!
+!--       Allocate index arrays as well as array containing height information
 !--       for soil.
           IF ( vmea(l)%soil_sampling )  THEN
@@ -868,6 +1283,6 @@
              ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
              ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
-             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
-          ENDIF     
+             ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil)  )
+          ENDIF
 !
 !--       For soil, store the grid indices.
@@ -885,8 +1300,8 @@
                             vmea(l)%j_soil(ns) = j
                             vmea(l)%k_soil(ns) = k
-                            vmea(l)%depth(ns)  = zs(k)
+                            vmea(l)%depth(ns)  = - zs(k)
                          ENDDO
                       ENDIF
-                      
+
                       IF ( surf_usm_h%start_index(j,i) <=                      &
                            surf_usm_h%end_index(j,i) )  THEN
@@ -897,5 +1312,5 @@
                             vmea(l)%j_soil(ns) = j
                             vmea(l)%k_soil(ns) = k
-                            vmea(l)%depth(ns)  = surf_usm_h%zw(k,m)
+                            vmea(l)%depth(ns)  = - surf_usm_h%zw(k,m)
                          ENDDO
                       ENDIF
@@ -907,5 +1322,5 @@
 !--       Allocate array to save the sampled values.
           ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
-          
+
           IF ( vmea(l)%soil_sampling )                                         &
              ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
@@ -924,19 +1339,21 @@
           IF ( ALLOCATED( n_utm_tmp ) )  DEALLOCATE( n_utm_tmp )
           IF ( ALLOCATED( n_utm )     )  DEALLOCATE( n_utm )
-          IF ( ALLOCATED( z_ag  )     )  DEALLOCATE( z_ag  )
-          IF ( ALLOCATED( z_ag  )     )  DEALLOCATE( vmea(l)%dim_t )
+          IF ( ALLOCATED( zar  )      )  DEALLOCATE( vmea(l)%dim_t )
+          IF ( ALLOCATED( zar  )      )  DEALLOCATE( zar  )
+          IF ( ALLOCATED( station_h ) )  DEALLOCATE( station_h )
+
        ENDIF
     ENDDO
 !
-!-- Close input file for virtual measurements. Therefore, just call 
-!-- the read attribute routine with the "close" option. 
-    CALL netcdf_data_input_att( vmea_general%nvm, char_numstations,            &
-                                vmea_general%id_vm, '',                        &
-                                global_attribute, 'close', '' )
-!
-!-- Sum-up the number of observation coordiates, for atmosphere first. 
+!-- Dellocate flag array
+    DEALLOCATE( meas_flag )
+!
+!-- Close input file for virtual measurements. 
+    CALL close_input_file( pids_id )
+!
+!-- Sum-up the number of observation coordiates, for atmosphere first.
 !-- This is actually only required for data output.
     ALLOCATE( ns_all(1:vmea_general%nvm) )
-    ns_all = 0    
+    ns_all = 0
 #if defined( __parallel )
     CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
@@ -948,5 +1365,5 @@
 !
 !-- Now for soil
-    ns_all = 0    
+    ns_all = 0
 #if defined( __parallel )
     CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
@@ -956,177 +1373,1159 @@
 #endif
     vmea(:)%ns_soil_tot = ns_all(:)
-    
+
     DEALLOCATE( ns_all )
-!                               
-!-- Dellocate flag array
-    DEALLOCATE( meas_flag )
-!
-!-- Initialize binary data output of virtual measurements.
-!-- Open binary output file.
-    CALL check_open( 27 )
-!
-!-- Output header information.
-    CALL vm_data_output
-        
+!
+!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to
+!-- be defined, so that each processor knows where to write the data.
+#if defined( __netcdf4_parallel )
+    ALLOCATE( ns_atmos(0:numprocs-1,1:vmea_general%nvm) )
+    ALLOCATE( ns_soil(0:numprocs-1,1:vmea_general%nvm)  )
+    ns_atmos = 0
+    ns_soil  = 0
+
+    DO  l = 1, vmea_general%nvm
+       ns_atmos(myid,l) = vmea(l)%ns
+       ns_soil(myid,l)  = vmea(l)%ns_soil
+    ENDDO
+
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm,   &
+                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm,    &
+                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
+#else
+    ns_atmos(0,:) = vmea(:)%ns
+    ns_soil(0,:)  = vmea(:)%ns_soil
+#endif
+
+!
+!-- Determine the start coordinate in NetCDF file for the local arrays.
+!-- Note, start coordinates are initialized with zero for sake of simplicity
+!-- in summation. However, in NetCDF the start coordinates must be >= 1,
+!-- so that a one needs to be added at the end. 
+    DO  l = 1, vmea_general%nvm
+       DO  n  = 0, myid - 1
+          vmea(l)%start_coord_a = vmea(l)%start_coord_a + ns_atmos(n,l)
+          vmea(l)%start_coord_s = vmea(l)%start_coord_s + ns_soil(n,l)
+       ENDDO
+!
+!--    Start coordinate in NetCDF starts always at one not at 0.
+       vmea(l)%start_coord_a = vmea(l)%start_coord_a + 1
+       vmea(l)%start_coord_s = vmea(l)%start_coord_s + 1
+!
+!--    Determine the local end coordinate
+       vmea(l)%end_coord_a = vmea(l)%start_coord_a + vmea(l)%ns - 1
+       vmea(l)%end_coord_s = vmea(l)%start_coord_s + vmea(l)%ns_soil - 1
+    ENDDO
+
+    DEALLOCATE( ns_atmos )
+    DEALLOCATE( ns_soil  )
+
+#endif
+
+#endif
+
   END SUBROUTINE vm_init
-  
-  
+
+
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Binary data output.
+!> Initialize output using data-output module
 !------------------------------------------------------------------------------!
-  SUBROUTINE vm_data_output
-    
-     USE pegrid
-    
-     IMPLICIT NONE
-         
-     INTEGER(iwp) ::  i  !< running index over IO blocks    
-     INTEGER(iwp) ::  l  !< running index over all stations
-     INTEGER(iwp) ::  n  !< running index over all measured variables at a station
-!
-!--  Header output on each PE
-     IF ( init )  THEN
-
-        DO  i = 0, io_blocks-1
-           IF ( i == io_group )  THEN
-              WRITE ( 27 )  'number of measurements            '
-              WRITE ( 27 )  vmea_general%nvm
-
-              DO  l = 1, vmea_general%nvm
-                 WRITE ( 27 )  'site                              '
-                 WRITE ( 27 )  vmea(l)%site
-                 WRITE ( 27 )  'file                              '
-                 WRITE ( 27 )  vmea(l)%filename_original
-                 WRITE ( 27 )  'feature_type                      '
-                 WRITE ( 27 )  vmea(l)%feature_type
-                 WRITE ( 27 )  'origin_x_obs                      '
-                 WRITE ( 27 )  vmea(l)%origin_x_obs
-                 WRITE ( 27 )  'origin_y_obs                      '
-                 WRITE ( 27 )  vmea(l)%origin_y_obs
-                 WRITE ( 27 )  'total number of observation points'                               
-                 WRITE ( 27 )  vmea(l)%ns_tot
-                 WRITE ( 27 )  'number of measured variables      '
-                 WRITE ( 27 )  vmea(l)%nmeas
-                 WRITE ( 27 )  'variables                         '
-                 WRITE ( 27 )  vmea(l)%measured_vars_name(:)
-                 WRITE ( 27 )  'number of observation points      '
-                 WRITE ( 27 )  vmea(l)%ns
-                 WRITE ( 27 )  'E_UTM                             '
-                 WRITE ( 27 )  init_model%origin_x                             &
+ SUBROUTINE vm_init_output
+
+    CHARACTER(LEN=100) ::  variable_name  !< name of output variable
+
+    INTEGER(iwp) ::  l              !< loop index
+    INTEGER(iwp) ::  n              !< loop index
+    INTEGER      ::  return_value   !< returned status value of called function
+
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ndim !< dummy to write dimension
+
+    REAL(wp) ::  dum_lat !< transformed geographical coordinate (latitude)
+    REAL(wp) ::  dum_lon !< transformed geographical coordinate (longitude)
+
+!
+!-- Determine the number of output timesteps. Set a maximum value of 80000
+!-- timesteps.
+    ntimesteps = MIN( CEILING(                                                 &
+                  ( end_time - MAX( vm_time_start, time_since_reference_point )&
+                  ) / MAX( 0.1_wp, dt_virtual_measurement ) ), ntimesteps_max )
+!
+!-- Create directory where output files will be stored.
+    CALL local_system( 'mkdir -p VM_OUTPUT' // TRIM( coupling_char ) )
+!
+!-- Loop over all sites.
+    DO  l = 1, vmea_general%nvm
+!
+!--    Skip if no observations will be taken for this site.
+       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
+!
+!--    Define output file.
+       WRITE( vmea(l)%nc_filename, '(A,I4.4)')  'VM_OUTPUT' //                 &
+                                                TRIM( coupling_char ) // '/' //&
+                                                'site', l
+
+
+       return_value = dom_def_file( vmea(l)%nc_filename, 'netcdf4-parallel' )
+!
+!--    Define global attributes.
+!--    Before, transform UTM into geographical coordinates.
+       CALL convert_utm_to_geographic( crs_list,                               &
+                                       vmea(l)%origin_x_obs,                   &
+                                       vmea(l)%origin_y_obs,                   &
+                                       dum_lon,                                &
+                                       dum_lat )
+
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'site',                    &
+                                   value = TRIM( vmea(l)%site ) )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'title',                   &
+                                   value = 'Virtual measurement output')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'source',                  &
+                                   value = 'PALM-4U')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'institution',             &
+                                   value = input_file_atts%institution )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'acronym',                 &
+                                   value = input_file_atts%acronym )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'author',                  &
+                                   value = input_file_atts%author )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'contact_person',          &
+                                   value = input_file_atts%contact_person )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'iop',                     &
+                                   value = input_file_atts%campaign )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'campaign',                &
+                                   value = 'PALM-4U' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_time ',            &
+                                   value = origin_date_time)
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'location',                &
+                                   value = input_file_atts%location )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_x',                &
+                                   value = vmea(l)%origin_x_obs )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_y',                &
+                                   value = vmea(l)%origin_y_obs )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_lon',              &
+                                   value = dum_lon )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_lat',              &
+                                   value = dum_lat )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'origin_z',                &
+                                   value = 0.0 )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'rotation_angle',          &
+                                   value = input_file_atts%rotation_angle )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'featureType',             &
+                                   value = TRIM( vmea(l)%feature_type_out ) )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'data_content',            &
+                                   value = TRIM( vmea(l)%data_content ) )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'creation_time',           &
+                                   value = input_file_atts%creation_time )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'version',                 &
+                                   value = 1 ) !input_file_atts%version )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'creation_time',           &
+                                   value = TRIM( vmea(l)%site ) )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'Conventions',             &
+                                   value = input_file_atts%conventions )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'dependencies',            &
+                                   value = input_file_atts%dependencies )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'history',                 &
+                                   value = input_file_atts%history )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'references',              &
+                                   value = input_file_atts%references )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'comment',                 &
+                                   value = input_file_atts%comment )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'keywords',                &
+                                   value = input_file_atts%keywords )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   attribute_name = 'licence',                 &
+                                   value = '[UC]2 Open Licence; see [UC]2 ' // &
+                                           'data policy available at ' //      &
+                                           'www.uc2-program.org/uc2_data_policy.pdf' )
+!
+!--    Define dimensions.
+!--    station
+       ALLOCATE( ndim(1:vmea(l)%ns_tot) )
+       DO  n = 1, vmea(l)%ns_tot
+          ndim(n) = n
+       ENDDO
+       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
+                                   dimension_name = 'station',                 &
+                                   output_type = 'int32',                      &
+                                   bounds = (/1_iwp, vmea(l)%ns_tot/),         &
+                                   values_int32 = ndim )
+       DEALLOCATE( ndim )
+!
+!--    ntime
+       ALLOCATE( ndim(1:ntimesteps) )
+       DO  n = 1, ntimesteps
+          ndim(n) = n
+       ENDDO
+
+       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
+                                   dimension_name = 'ntime',                   &
+                                   output_type = 'int32',                      &
+                                   bounds = (/1_iwp, ntimesteps/),             &
+                                   values_int32 = ndim )
+       DEALLOCATE( ndim )
+!
+!--    nv
+       ALLOCATE( ndim(1:2) )
+       DO  n = 1, 2
+          ndim(n) = n
+       ENDDO
+
+       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
+                                   dimension_name = 'nv',                      &
+                                   output_type = 'int32',                      &
+                                   bounds = (/1_iwp, 2_iwp/),                  &
+                                   values_int32 = ndim )
+       DEALLOCATE( ndim )
+!
+!--    maximum name length
+       ALLOCATE( ndim(1:maximum_name_length) )
+       DO  n = 1, maximum_name_length
+          ndim(n) = n
+       ENDDO
+
+       return_value = dom_def_dim( vmea(l)%nc_filename,                        &
+                                   dimension_name = 'max_name_len',            &
+                                   output_type = 'int32',                      &
+                                   bounds = (/1_iwp, 32_iwp/),                 &
+                                   values_int32 = ndim )
+       DEALLOCATE( ndim )
+!
+!--    Define coordinate variables.
+!--    time
+       variable_name = 'time'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/ 'station  ',           &
+                                                        'ntime    '/),         &
+                                   output_type = 'real32' )
+!
+!--    station_name. DOM needs to be enabled to define CHARACTER variables.
+!        variable_name = 'station_name'
+!        return_value = dom_def_var( vmea(l)%nc_filename,                        &
+!                                    variable_name = variable_name,              &
+!                                    dimension_names = (/ 'max_name_len',        &
+!                                                         'station     '/),      &
+!                                    output_type = 'char' )
+!
+!--    vrs (vertical reference system)
+       variable_name = 'vrs'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/ 'station' /),          &
+                                   output_type = 'int8' )
+!
+!--    crs (coordinate reference system)
+       variable_name = 'crs'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/ 'station' /),          &
+                                   output_type = 'int8' )
+!
+!--    z
+       variable_name = 'z'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    station_h
+       variable_name = 'station_h'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    x
+       variable_name = 'x'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    y
+       variable_name = 'y'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    E-UTM
+       variable_name = 'E_UTM'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    N-UTM
+       variable_name = 'N_UTM'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    latitude
+       variable_name = 'lat'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    longitude
+       variable_name = 'lon'
+       return_value = dom_def_var( vmea(l)%nc_filename,                        &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'station'/),            &
+                                   output_type = 'real32' )
+!
+!--    Set attributes for the coordinate variables. Note, not all coordinates
+!--    have the same number of attributes.
+!--    Units
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = char_unit,                 &
+                                   value = 'seconds since ' // origin_date_time )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'z',                        &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_h',                &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'x',                        &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'y',                        &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'E_UTM',                    &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'N_UTM',                    &
+                                   attribute_name = char_unit,                 &
+                                   value = 'm' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lat',                      &
+                                   attribute_name = char_unit,                 &
+                                   value = 'degrees_north' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lon',                      &
+                                   attribute_name = char_unit,                 &
+                                   value = 'degrees_east' )
+!
+!--    Long name
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_name',             &
+                                   attribute_name = char_long,                 &
+                                   value = 'station name')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = char_long,                 &
+                                   value = 'time')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'z',                        &
+                                   attribute_name = char_long,                 &
+                                   value = 'height above origin' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_h',                &
+                                   attribute_name = char_long,                 &
+                                   value = 'surface altitude' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'x',                        &
+                                   attribute_name = char_long,                 &
+                                   value = 'distance to origin in x-direction' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'y',                        &
+                                   attribute_name = char_long,                 &
+                                   value = 'distance to origin in y-direction' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'E_UTM',                    &
+                                   attribute_name = char_long,                 &
+                                   value = 'easting' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'N_UTM',                    &
+                                   attribute_name = char_long,                 &
+                                   value = 'northing' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lat',                      &
+                                   attribute_name = char_long,                 &
+                                   value = 'latitude' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lon',                      &
+                                   attribute_name = char_long,                 &
+                                   value = 'longitude' )
+!
+!--    Standard name
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_name',             &
+                                   attribute_name = char_standard,             &
+                                   value = 'platform_name')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = char_standard,             &
+                                   value = 'time')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'z',                        &
+                                   attribute_name = char_standard,             &
+                                   value = 'height_above_mean_sea_level' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_h',                &
+                                   attribute_name = char_standard,             &
+                                   value = 'surface_altitude' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'E_UTM',                    &
+                                   attribute_name = char_standard,             &
+                                   value = 'projection_x_coordinate' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'N_UTM',                    &
+                                   attribute_name = char_standard,             &
+                                   value = 'projection_y_coordinate' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lat',                      &
+                                   attribute_name = char_standard,             &
+                                   value = 'latitude' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'lon',                      &
+                                   attribute_name = char_standard,             &
+                                   value = 'longitude' )
+!
+!--    Axis
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = 'axis',                    &
+                                   value = 'T')
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'z',                        &
+                                   attribute_name = 'axis',                    &
+                                   value = 'Z' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'x',                        &
+                                   attribute_name = 'axis',                    &
+                                   value = 'X' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'y',                        &
+                                   attribute_name = 'axis',                    &
+                                   value = 'Y' )
+!
+!--    Set further individual attributes for the coordinate variables.
+!--    For station name
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'station_name',             &
+                                   attribute_name = 'cf_role',                 &
+                                   value = 'timeseries_id' )
+!
+!--    For time
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = 'calendar',                &
+                                   value = 'proleptic_gregorian' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'time',                     &
+                                   attribute_name = 'bounds',                  &
+                                   value = 'time_bounds' )
+!
+!--    For vertical reference system
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'vrs',                      &
+                                   attribute_name = char_long,                 &
+                                   value = 'vertical reference system' )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'vrs',                      &
+                                   attribute_name = 'system_name',             &
+                                   value = 'DHHN2016' )
+!
+!--    For z
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'z',                        &
+                                   attribute_name = 'positive',                &
+                                   value = 'up' )
+!
+!--    For coordinate reference system
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'epsg_code',               &
+                                   value = coord_ref_sys%epsg_code )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'false_easting',           &
+                                   value = coord_ref_sys%false_easting )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'false_northing',          &
+                                   value = coord_ref_sys%false_northing )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'grid_mapping_name',       &
+                                   value = coord_ref_sys%grid_mapping_name )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'inverse_flattening',      &
+                                   value = coord_ref_sys%inverse_flattening )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'latitude_of_projection_origin',&
+                                   value = coord_ref_sys%latitude_of_projection_origin )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = char_long,                 &
+                                   value = coord_ref_sys%long_name )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'longitude_of_central_meridian', &
+                                   value = coord_ref_sys%longitude_of_central_meridian )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'longitude_of_prime_meridian', &
+                                   value = coord_ref_sys%longitude_of_prime_meridian )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'scale_factor_at_central_meridian', &
+                                   value = coord_ref_sys%scale_factor_at_central_meridian )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = 'semi_major_axis',         &
+                                   value = coord_ref_sys%semi_major_axis )
+       return_value = dom_def_att( vmea(l)%nc_filename,                        &
+                                   variable_name = 'crs',                      &
+                                   attribute_name = char_unit,                 &
+                                   value = coord_ref_sys%units )
+!
+!--    In case of sampled soil quantities, define further dimensions and
+!--    coordinates.
+       IF ( vmea(l)%soil_sampling )  THEN
+!
+!--       station for soil
+          ALLOCATE( ndim(1:vmea(l)%ns_soil_tot) )
+          DO  n = 1, vmea(l)%ns_soil_tot
+             ndim(n) = n
+          ENDDO
+
+          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
+                                      dimension_name = 'station_soil',         &
+                                      output_type = 'int32',                   &
+                                      bounds = (/1_iwp,vmea(l)%ns_soil_tot/),  &
+                                      values_int32 = ndim )
+          DEALLOCATE( ndim )
+!
+!--       ntime for soil
+          ALLOCATE( ndim(1:ntimesteps) )
+          DO  n = 1, ntimesteps
+             ndim(n) = n
+          ENDDO
+
+          return_value = dom_def_dim( vmea(l)%nc_filename,                     &
+                                      dimension_name = 'ntime_soil',           &
+                                      output_type = 'int32',                   &
+                                      bounds = (/1_iwp,ntimesteps/),           &
+                                      values_int32 = ndim )
+          DEALLOCATE( ndim )
+!
+!--       time for soil
+          variable_name = 'time_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil',      &
+                                                          'ntime_soil  '/),    &
+                                      output_type = 'real32' )
+!
+!--       station_name for soil
+          variable_name = 'station_name_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/ 'max_name_len',     &
+                                                           'station_soil'/),   &
+                                      output_type = 'char' )
+!
+!--       z
+          variable_name = 'z_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       station_h for soil
+          variable_name = 'station_h_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       x soil
+          variable_name = 'x_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!-        y soil
+          variable_name = 'y_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       E-UTM soil
+          variable_name = 'E_UTM_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       N-UTM soil
+          variable_name = 'N_UTM_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       latitude soil
+          variable_name = 'lat_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       longitude soil
+          variable_name = 'lon_soil'
+          return_value = dom_def_var( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      dimension_names = (/'station_soil'/),    &
+                                      output_type = 'real32' )
+!
+!--       Set attributes for the coordinate variables. Note, not all coordinates
+!--       have the same number of attributes.
+!--       Units
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = char_unit,              &
+                                      value = 'seconds since ' // origin_date_time )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'z_soil',                &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_h_soil',        &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'x_soil',                &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'y_soil',                &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'E_UTM_soil',            &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'N_UTM_soil',            &
+                                      attribute_name = char_unit,              &
+                                      value = 'm' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lat_soil',              &
+                                      attribute_name = char_unit,              &
+                                      value = 'degrees_north' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lon_soil',              &
+                                      attribute_name = char_unit,              &
+                                      value = 'degrees_east' )
+!
+!--       Long name
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_name_soil',     &
+                                      attribute_name = char_long,              &
+                                      value = 'station name')
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = char_long,              &
+                                      value = 'time')
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'z_soil',                &
+                                      attribute_name = char_long,              &
+                                      value = 'height above origin' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_h_soil',        &
+                                      attribute_name = char_long,              &
+                                      value = 'surface altitude' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'x_soil',                &
+                                      attribute_name = char_long,              &
+                                      value = 'distance to origin in x-direction' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'y_soil',                &
+                                      attribute_name = char_long,              &
+                                      value = 'distance to origin in y-direction' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'E_UTM_soil',            &
+                                      attribute_name = char_long,              &
+                                      value = 'easting' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'N_UTM_soil',            &
+                                      attribute_name = char_long,              &
+                                      value = 'northing' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lat_soil',              &
+                                      attribute_name = char_long,              &
+                                      value = 'latitude' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lon_soil',              &
+                                      attribute_name = char_long,              &
+                                      value = 'longitude' )
+!
+!--       Standard name
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_name_soil',     &
+                                      attribute_name = char_standard,          &
+                                      value = 'platform_name')
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = char_standard,          &
+                                      value = 'time')
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'z_soil',                &
+                                      attribute_name = char_standard,          &
+                                      value = 'height_above_mean_sea_level' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_h_soil',        &
+                                      attribute_name = char_standard,          &
+                                      value = 'surface_altitude' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'E_UTM_soil',            &
+                                      attribute_name = char_standard,          &
+                                      value = 'projection_x_coordinate' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'N_UTM_soil',            &
+                                      attribute_name = char_standard,          &
+                                      value = 'projection_y_coordinate' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lat_soil',              &
+                                      attribute_name = char_standard,          &
+                                      value = 'latitude' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'lon_soil',              &
+                                      attribute_name = char_standard,          &
+                                      value = 'longitude' )
+!
+!--       Axis
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = 'axis',                 &
+                                      value = 'T')
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'z_soil',                &
+                                      attribute_name = 'axis',                 &
+                                      value = 'Z' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'x_soil',                &
+                                      attribute_name = 'axis',                 &
+                                      value = 'X' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'y_soil',                &
+                                      attribute_name = 'axis',                 &
+                                      value = 'Y' )
+!
+!--       Set further individual attributes for the coordinate variables.
+!--       For station name soil
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'station_name_soil',     &
+                                      attribute_name = 'cf_role',              &
+                                      value = 'timeseries_id' )
+!
+!--       For time soil
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = 'calendar',             &
+                                      value = 'proleptic_gregorian' )
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'time_soil',             &
+                                      attribute_name = 'bounds',               &
+                                      value = 'time_bounds' )
+!
+!--       For z soil
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = 'z_soil',                &
+                                      attribute_name = 'positive',             &
+                                      value = 'up' )
+       ENDIF
+!
+!--    Define variables that shall be sampled.
+       DO  n = 1, vmea(l)%nmeas
+          variable_name = TRIM( vmea(l)%var_atts(n)%name )
+!
+!--       In order to link the correct dimension names, atmosphere and soil
+!--       variables need to be distinguished.
+          IF ( vmea(l)%soil_sampling  .AND.                                    &
+               ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
+
+             return_value = dom_def_var( vmea(l)%nc_filename,                  &
+                                         variable_name = variable_name,        &
+                                         dimension_names = (/'station_soil',   &
+                                                             'ntime_soil  '/), &
+                                         output_type = 'real32' )
+          ELSE
+
+             return_value = dom_def_var( vmea(l)%nc_filename,                  &
+                                         variable_name = variable_name,        &
+                                         dimension_names = (/'station',        &
+                                                             'ntime  '/),      &
+                                         output_type = 'real32' )
+          ENDIF
+!
+!--       Set variable attributes. Please note, for some variables not all
+!--       attributes are defined, e.g. standard_name for the horizontal wind
+!--       components.
+          CALL vm_set_attributes( vmea(l)%var_atts(n) )
+
+          IF ( vmea(l)%var_atts(n)%long_name /= 'none' )  THEN
+             return_value = dom_def_att( vmea(l)%nc_filename,                  &
+                                         variable_name = variable_name,        &
+                                         attribute_name = char_long,           &
+                                         value = TRIM( vmea(l)%var_atts(n)%long_name ) )
+          ENDIF
+          IF ( vmea(l)%var_atts(n)%standard_name /= 'none' )  THEN
+             return_value = dom_def_att( vmea(l)%nc_filename,                  &
+                                         variable_name = variable_name,        &
+                                         attribute_name = char_standard,       &
+                                         value = TRIM( vmea(l)%var_atts(n)%standard_name ) )
+          ENDIF
+          IF ( vmea(l)%var_atts(n)%units /= 'none' )  THEN
+             return_value = dom_def_att( vmea(l)%nc_filename,                  &
+                                         variable_name = variable_name,        &
+                                         attribute_name = char_unit,           &
+                                         value = TRIM( vmea(l)%var_atts(n)%units ) )
+          ENDIF
+
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      attribute_name = 'grid_mapping',         &
+                                      value = TRIM( vmea(l)%var_atts(n)%grid_mapping ) )
+
+          return_value = dom_def_att( vmea(l)%nc_filename,                     &
+                                      variable_name = variable_name,           &
+                                      attribute_name = 'coordinates',          &
+                                      value = TRIM( vmea(l)%var_atts(n)%coordinates ) )
+
+!           return_value = dom_def_att( vmea(l)%nc_filename,                     &
+!                                       variable_name = variable_name,           &
+!                                       attribute_name = char_fill,              &
+!                                       value = vmea(l)%var_atts(n)%fill_value )
+
+       ENDDO  ! loop over variables per site
+
+    ENDDO  ! loop over sites
+
+
+ END SUBROUTINE vm_init_output
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parallel NetCDF output via data-output module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE vm_data_output
+
+    CHARACTER(LEN=100) ::  variable_name !< name of output variable
+
+    INTEGER(iwp)       ::  l             !< loop index
+    INTEGER(iwp)       ::  n             !< loop index
+    INTEGER            ::  return_value  !< returned status value of called function
+
+    INTEGER(iwp)       ::  t_ind         !< time index
+
+    LOGICAL            ::  proceed_run   !< flag to indicate that the maximum number of
+                                         !< timesteps in the file will be exceeded
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE           ::  oro_rel                  !< relative altitude of model surface
+    REAL(wp), DIMENSION(:), POINTER               ::  output_values_1d_pointer !< pointer for 1d output array
+    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET   ::  output_values_1d_target  !< target for 1d output array
+    REAL(wp), DIMENSION(:,:), POINTER             ::  output_values_2d_pointer !< pointer for 2d output array
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET ::  output_values_2d_target  !< target for 2d output array
+
+!
+!-- At the first call of this routine write the spatial coordinates.
+    IF ( .NOT. initial_write_coordinates )  THEN
+!
+!--    Write spatial coordinates.
+       DO  l = 1, vmea_general%nvm
+!
+!--       Skip if no observations were taken.
+          IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
+
+          ALLOCATE( output_values_1d_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
+!
+!--       Output of Easting coordinate. Before output, recalculate EUTM.
+          output_values_1d_target = init_model%origin_x                        &
                     + REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
                     * COS( init_model%rotation_angle * pi / 180.0_wp )         &
                     + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
                     * SIN( init_model%rotation_angle * pi / 180.0_wp )
-                 WRITE ( 27 )  'N_UTM                             '
-                 WRITE ( 27 )  init_model%origin_y                             &
+
+          output_values_1d_pointer => output_values_1d_target
+
+          return_value =                                                       &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'E_UTM',                                      &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
+!
+!--       Output of Northing coordinate. Before output, recalculate NUTM.
+          output_values_1d_target = init_model%origin_y                        &
                     - REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
                     * SIN( init_model%rotation_angle * pi / 180.0_wp )         &
                     + REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
                     * COS( init_model%rotation_angle * pi / 180.0_wp )
-                 WRITE ( 27 )  'Z_AG                              '
-                 WRITE ( 27 )  vmea(l)%z_ag(1:vmea(l)%ns)
-                 WRITE ( 27 )  'soil sampling                     '
-                 WRITE ( 27 )  MERGE( 'yes                               ',    &
-                                      'no                                ',    &
-                                      vmea(l)%soil_sampling )
- 
-                 IF ( vmea(l)%soil_sampling )  THEN                 
-                    WRITE ( 27 )  'total number of soil points       '                               
-                    WRITE ( 27 )  vmea(l)%ns_soil_tot
-                    WRITE ( 27 )  'number of soil points             '
-                    WRITE ( 27 )  vmea(l)%ns_soil
-                    WRITE ( 27 )  'E_UTM soil                        '
-                    WRITE ( 27 )  init_model%origin_x                          &
-                       + REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
-                               KIND = wp ) * dx                                &
-                       * COS( init_model%rotation_angle * pi / 180.0_wp )      &
-                       + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
-                               KIND = wp ) * dy                                &
-                       * SIN( init_model%rotation_angle * pi / 180.0_wp )
-                    WRITE ( 27 )  'N_UTM soil                        '
-                    WRITE ( 27 )  init_model%origin_y                          &
-                       - REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
-                               KIND = wp ) * dx                                &
-                       * SIN( init_model%rotation_angle * pi / 180.0_wp )      &
-                       + REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,     &
-                               KIND = wp ) * dy                                &
-                       * COS( init_model%rotation_angle * pi / 180.0_wp )
-                    WRITE ( 27 )  'DEPTH                             '
-                    WRITE ( 27 )  vmea(l)%depth(1:vmea(l)%ns_soil)
-                 ENDIF
-              ENDDO
-
-           ENDIF
-        ENDDO
-        
+
+          output_values_1d_pointer => output_values_1d_target
+          return_value =                                                       &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'N_UTM',                                      &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
+!
+!--       Output of relative height coordinate.
+!--       Before this is output, first define the relative orographie height
+!--       and add this to z.
+          ALLOCATE( oro_rel(1:vmea(l)%ns) )
+          DO  n = 1, vmea(l)%ns
+             oro_rel(n) = zw(topo_top_ind(vmea(l)%j(n),vmea(l)%i(n),3))
+          ENDDO
+
+          output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
+          output_values_1d_pointer => output_values_1d_target
+          return_value =                                                       &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'z',                                          &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
+!
+!--       Write surface altitude for the station. Note, since z is already
+!--       a relative observation height, station_h must be zero, in order
+!--       to obtain the observation level.
+          output_values_1d_target = oro_rel(:)
+          output_values_1d_pointer => output_values_1d_target
+          return_value =                                                       &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'station_h',                                  &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_a/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_a  /) )
+
+          DEALLOCATE( oro_rel )
+          DEALLOCATE( output_values_1d_target )
+!
+!--       In case of sampled soil quantities, output also the respective 
+!--       coordinate arrays.
+          IF ( vmea(l)%soil_sampling )  THEN
+             ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
+!
+!--          Output of Easting coordinate. Before output, recalculate EUTM.
+             output_values_1d_target = init_model%origin_x                     &
+               + REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
+               * COS( init_model%rotation_angle * pi / 180.0_wp )              &
+               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
+               * SIN( init_model%rotation_angle * pi / 180.0_wp )
+             output_values_1d_pointer => output_values_1d_target
+             return_value =                                                    &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'E_UTM_soil',                                 &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
+!
+!--          Output of Northing coordinate. Before output, recalculate NUTM.
+             output_values_1d_target = init_model%origin_y                     &
+               - REAL( vmea(l)%i(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
+               * SIN( init_model%rotation_angle * pi / 180.0_wp )              &
+               + REAL( vmea(l)%j(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
+               * COS( init_model%rotation_angle * pi / 180.0_wp )
+
+             output_values_1d_pointer => output_values_1d_target
+             return_value =                                                    &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'N_UTM_soil',                                 &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
+!
+!--          Output of relative height coordinate.
+!--          Before this is output, first define the relative orographie height
+!--          and add this to z.
+             ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
+             DO  n = 1, vmea(l)%ns
+                oro_rel(n) = zw(topo_top_ind(vmea(l)%j_soil(n),vmea(l)%i_soil(n),3))
+             ENDDO
+
+             output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
+             output_values_1d_pointer => output_values_1d_target
+             return_value =                                                    &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'z_soil',                                     &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
+!
+!--          Write surface altitude for the station. Note, since z is already
+!--          a relative observation height, station_h must be zero, in order
+!--          to obtain the observation level.
+             output_values_1d_target = oro_rel(:)
+             output_values_1d_pointer => output_values_1d_target
+             return_value =                                                    &
+                  dom_write_var( vmea(l)%nc_filename,                          &
+                                 'station_h_soil',                             &
+                                 values_realwp_1d = output_values_1d_pointer,  &
+                                 bounds_start = (/vmea(l)%start_coord_s/),     &
+                                 bounds_end   = (/vmea(l)%end_coord_s  /) )
+
+             DEALLOCATE( oro_rel )
+             DEALLOCATE( output_values_1d_target )
+!
+!--          Write the stations name
+             
+          ENDIF
+
+       ENDDO  ! loop over sites
+
+       initial_write_coordinates = .TRUE.
+    ENDIF
+!
+!-- Loop over all sites.
+    DO  l = 1, vmea_general%nvm
+!
+!--    Skip if no observations were taken.
+       IF ( vmea(l)%ns_tot == 0  .AND.  vmea(l)%ns_soil_tot == 0 )  CYCLE
+!
+!--    Determine time index in file.
+       t_ind = vmea(l)%file_time_index + 1
+!
+!--    Write output variables. Distinguish between atmosphere and soil variables.
+       DO  n = 1, vmea(l)%nmeas
+          IF ( vmea(l)%soil_sampling  .AND.                                    &
+            ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) )  THEN
+!
+!--          Write time coordinate to file
+             variable_name = 'time_soil'
+             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
+             output_values_2d_target(t_ind,:) = time_since_reference_point
+             output_values_2d_pointer => output_values_2d_target
+
+             return_value = dom_write_var( vmea(l)%nc_filename,                &
+                                           variable_name,                      &
+                                           values_realwp_2d = output_values_2d_pointer, &
+                                           bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
+                                           bounds_end   = (/vmea(l)%end_coord_s, t_ind /) )
+
+             variable_name = TRIM( vmea(l)%var_atts(n)%name )
+             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
+             output_values_2d_pointer => output_values_2d_target
+             return_value =                                                    &
+                      dom_write_var( vmea(l)%nc_filename,                      &
+                                     variable_name,                            &
+                                     values_realwp_2d = output_values_2d_pointer, &
+                                     bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
+                                     bounds_end   = (/vmea(l)%end_coord_s, t_ind  /) )
+             DEALLOCATE( output_values_2d_target )
+          ELSE
+!
+!--          Write time coordinate to file
+             variable_name = 'time'
+             ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
+             output_values_2d_target(t_ind,:) = time_since_reference_point
+             output_values_2d_pointer => output_values_2d_target
+
+             return_value = dom_write_var( vmea(l)%nc_filename,                &
+                                           variable_name,                      &
+                                           values_realwp_2d = output_values_2d_pointer, &
+                                           bounds_start = (/vmea(l)%start_coord_a, t_ind/), &
+                                           bounds_end   = (/vmea(l)%end_coord_a, t_ind/) )
+
+             variable_name = TRIM( vmea(l)%var_atts(n)%name )
+
+             output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
+             output_values_2d_pointer => output_values_2d_target
+             return_value =                                                    &
+                      dom_write_var( vmea(l)%nc_filename,                      &
+                                     variable_name,                            &
+                                     values_realwp_2d = output_values_2d_pointer, &
+                                     bounds_start = (/ vmea(l)%start_coord_a, t_ind /), &
+                                     bounds_end   = (/ vmea(l)%end_coord_a, t_ind /) )
+
+             DEALLOCATE( output_values_2d_target )
+          ENDIF
+       ENDDO
+!
+!--    Update number of written time indices
+       vmea(l)%file_time_index = t_ind
+
+    ENDDO  ! loop over sites
+!
+!-- Check if the time index exceeds its maximum values. This case, terminate the
+!-- run and force a restart. This is simply done by creating the file
+!-- DO_RESTART_NOW.
+    proceed_run = .TRUE.
+    IF ( ANY( vmea(:)%file_time_index >= ntimesteps - 1) )  proceed_run = .FALSE.
+
 #if defined( __parallel )
-        CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( MPI_IN_PLACE,                                          &
+                        proceed_run,                                           &
+                        1,                                                     &
+                        MPI_LOGICAL,                                           &
+                        MPI_LAND,                                              &
+                        MPI_COMM_WORLD,                                        &
+                        ierr )
 #endif
 !
-!--     After header information is written, set control flag to .FALSE.
-        init = .FALSE.
-!
-!--  Data output at each measurement timestep on each PE
-     ELSE
-        DO  i = 0, io_blocks-1
-
-           IF ( i == io_group )  THEN
-              WRITE( 27 )  'output time                       '
-              WRITE( 27 )  time_since_reference_point
-              DO  l = 1, vmea_general%nvm
-!
-!--              Skip binary writing if no observation points are defined on PE
-                 IF ( vmea(l)%ns < 1  .AND.  vmea(l)%ns_soil < 1)  CYCLE                 
-                 DO  n = 1, vmea(l)%nmeas
-                    WRITE( 27 )  vmea(l)%measured_vars_name(n)
-                    IF ( vmea(l)%soil_sampling  .AND.                           &
-                         ANY( TRIM( vmea(l)%measured_vars_name(n))  ==          &
-                              soil_vars ) )  THEN                   
-                       WRITE( 27 )  vmea(l)%measured_vars_soil(:,n)
-                    ELSE
-                       WRITE( 27 )  vmea(l)%measured_vars(:,n)
-                    ENDIF
-                 ENDDO
-           
-              ENDDO
-           ENDIF
-        ENDDO
-#if defined( __parallel )
-        CALL MPI_BARRIER( comm2d, ierr )
-#endif
-     ENDIF
-  
-  END SUBROUTINE vm_data_output  
-  
-  
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Write end-of-file statement as last action.
-!------------------------------------------------------------------------------!
-  SUBROUTINE vm_last_actions
-    
-     USE pegrid
-    
-     IMPLICIT NONE
-         
-     INTEGER(iwp) ::  i         !< running index over IO blocks    
- 
-     DO  i = 0, io_blocks-1
-        IF ( i == io_group )  THEN
-           WRITE( 27 )  'EOF                               '
-        ENDIF
-     ENDDO
-#if defined( __parallel )
-        CALL MPI_BARRIER( comm2d, ierr )
-#endif
-!
-!--  Close binary file
-     CALL close_file( 27 )
-  
-  END SUBROUTINE vm_last_actions 
-  
+!-- Create file DO_RESTART_NOW to force a restart. This file is only created
+!-- by rank 0 of the root domain.
+    IF ( myid == 0  .AND.  .NOT. child_domain  .AND.  .NOT. proceed_run )  THEN
+       OPEN( 999, FILE='DO_RESTART_NOW' )
+       CLOSE( 999 )
+
+       message_string = 'Run will be terminated because virtual ' //           &
+                        'measurement times index exceeds its maximum value ' //&
+                        'in the output files. A restart run is forced.' 
+       CALL message( 'vm_data_output', 'PA0707', 0, 0, 0, 6, 0 )
+    ENDIF
+
+  END SUBROUTINE vm_data_output
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1136,15 +2535,12 @@
   SUBROUTINE vm_sampling
 
-    USE arrays_3d,                                                             &
-        ONLY:  exner, pt, q, u, v, w
-
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation  
+        ONLY:  radiation
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
-    
-     IMPLICIT NONE
-      
+        ONLY:  surf_def_h,                                                     &
+               surf_lsm_h,                                                     &
+               surf_usm_h
+
      INTEGER(iwp) ::  i         !< grid index in x-direction
      INTEGER(iwp) ::  j         !< grid index in y-direction
@@ -1156,23 +2552,27 @@
      INTEGER(iwp) ::  n         !< running index over all measured variables at a station
      INTEGER(iwp) ::  nn        !< running index over the number of chemcal species
-     
+
      LOGICAL ::  match_lsm !< flag indicating natural-type surface
      LOGICAL ::  match_usm !< flag indicating urban-type surface
-!
-!--  Loop over all sites. 
+
+     REAL(wp) ::  e_s      !< saturation water vapor pressure
+     REAL(wp) ::  q_s      !< saturation mixing ratio
+     REAL(wp) ::  q_wv     !< mixing ratio
+!
+!--  Loop over all sites.
      DO  l = 1, vmea_general%nvm
 !
 !--     At the beginning, set _FillValues
         IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
-           vmea(l)%measured_vars      = vmea(l)%fillout 
+           vmea(l)%measured_vars      = vmea(l)%fillout
         IF ( ALLOCATED( vmea(l)%measured_vars_soil ) )                         &
-           vmea(l)%measured_vars_soil = vmea(l)%fillout 
-!
-!--     Loop over all variables measured at this site.  
+           vmea(l)%measured_vars_soil = vmea(l)%fillout
+!
+!--     Loop over all variables measured at this site.
         DO  n = 1, vmea(l)%nmeas
-        
-           SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
-           
-              CASE ( 'theta' )
+
+           SELECT CASE ( TRIM( vmea(l)%var_atts(n)%name ) )
+
+              CASE ( 'theta' ) ! potential temperature
                  IF ( .NOT. neutral )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -1183,6 +2583,6 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'ta' )
+
+              CASE ( 'ta' ) ! absolute temperature
                  IF ( .NOT. neutral )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -1190,11 +2590,12 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
+                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )     &
+                                                  - degc_to_k
                     ENDDO
                  ENDIF
-              
+
               CASE ( 't_va' )
-                 
-              CASE ( 'hus', 'haa' )
+
+              CASE ( 'hus' ) ! mixing ratio
                  IF ( humidity )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -1205,6 +2606,47 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'u', 'ua' )
+
+              CASE ( 'haa' ) ! absolute humidity
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
+                                                    / ( 1.0_wp - q(k,j,i) ) )  &
+                                                  * rho_air(k)
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'pwv' ) ! water vapor partial pressure
+                 IF ( humidity )  THEN
+!                     DO  m = 1, vmea(l)%ns
+!                        k = vmea(l)%k(m)
+!                        j = vmea(l)%j(m)
+!                        i = vmea(l)%i(m)
+!                        vmea(l)%measured_vars(m,n) = ( q(k,j,i)                 &
+!                                                     / ( 1.0_wp - q(k,j,i) ) )  &
+!                                                   * rho_air(k)
+!                     ENDDO
+                 ENDIF
+
+              CASE ( 'hur' ) ! relative humidity
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+!
+!--                    Calculate actual temperature, water vapor saturation
+!--                    pressure, and based on this  the saturation mixing ratio.
+                       e_s  = magnus( exner(k) * pt(k,j,i) )
+                       q_s  = rd_d_rv * e_s / ( hyp(k) - e_s )
+                       q_wv = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
+
+                       vmea(l)%measured_vars(m,n) = q_wv / ( q_s + 1E-10_wp )
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'u', 'ua' ) ! u-component
                  DO  m = 1, vmea(l)%ns
                     k = vmea(l)%k(m)
@@ -1213,6 +2655,6 @@
                     vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
                  ENDDO
-                 
-              CASE ( 'v', 'va' )
+
+              CASE ( 'v', 'va' ) ! v-component
                  DO  m = 1, vmea(l)%ns
                     k = vmea(l)%k(m)
@@ -1221,14 +2663,14 @@
                     vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
                  ENDDO
-                 
-              CASE ( 'w' )
+
+              CASE ( 'w' ) ! w-component
                  DO  m = 1, vmea(l)%ns
-                    k = MAX ( 1, vmea(l)%k(m) )  
+                    k = MAX ( 1, vmea(l)%k(m) )
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
                     vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
                  ENDDO
-                 
-              CASE ( 'wspeed' )
+
+              CASE ( 'wspeed' ) ! horizontal wind speed
                  DO  m = 1, vmea(l)%ns
                     k = vmea(l)%k(m)
@@ -1240,17 +2682,18 @@
                                                      )
                  ENDDO
-                 
-              CASE ( 'wdir' )
+
+              CASE ( 'wdir' ) ! wind direction
                  DO  m = 1, vmea(l)%ns
                     k = vmea(l)%k(m)
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    
-                    vmea(l)%measured_vars(m,n) = ATAN2(                        &
-                                       - 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ),   &
-                                       - 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )    &
-                                                      ) * 180.0_wp / pi
+
+                    vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi      &
+                                               * ATAN2(                        &
+                                          0.5_wp * ( v(k,j,i) + v(k,j+1,i) ),  &
+                                          0.5_wp * ( u(k,j,i) + u(k,j,i+1) )   &
+                                                      )
                  ENDDO
-                 
+
               CASE ( 'utheta' )
                  DO  m = 1, vmea(l)%ns
@@ -1262,5 +2705,5 @@
                                                    pt(k,j,i)
                  ENDDO
-                 
+
               CASE ( 'vtheta' )
                  DO  m = 1, vmea(l)%ns
@@ -1272,5 +2715,5 @@
                                                    pt(k,j,i)
                  ENDDO
-                 
+
               CASE ( 'wtheta' )
                  DO  m = 1, vmea(l)%ns
@@ -1282,5 +2725,41 @@
                                                    pt(k,j,i)
                  ENDDO
-                 
+
+              CASE ( 'uqv' )
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
+                                                    ( u(k,j,i) + u(k,j,i+1) ) *&
+                                                      q(k,j,i)
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'vqv' )
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
+                                                    ( v(k,j,i) + v(k,j+1,i) ) *&
+                                                      q(k,j,i)
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'wqv' )
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = MAX ( 1, vmea(l)%k(m) )
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(m,n) = 0.5_wp *                   &
+                                                    ( w(k-1,j,i) + w(k,j,i) ) *&
+                                                      q(k,j,i)
+                    ENDDO
+                 ENDIF
+
               CASE ( 'uw' )
                  DO  m = 1, vmea(l)%ns
@@ -1292,5 +2771,5 @@
                                                  ( u(k,j,i)   + u(k,j,i+1) )
                  ENDDO
-                 
+
               CASE ( 'vw' )
                  DO  m = 1, vmea(l)%ns
@@ -1302,8 +2781,8 @@
                                                  ( v(k,j,i)   + v(k,j+1,i) )
                  ENDDO
-                 
+
               CASE ( 'uv' )
                  DO  m = 1, vmea(l)%ns
-                    k = MAX ( 1, vmea(l)%k(m) )
+                    k = vmea(l)%k(m)
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
@@ -1313,25 +2792,29 @@
                  ENDDO
 !
-!--           List of variables may need extension.
-              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2',    & 
-                     'tro3' )                     
+!--           Chemistry variables. List of variables may need extension.
+!--           Note, gas species in PALM are in ppm and no distinction is made
+!--           between mole-fraction and concentration quantities (all are
+!--           output in ppm so far).
+              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2',            &
+                     'mcno', 'mcno2', 'tro3' )
                  IF ( air_chemistry )  THEN
 !
-!--                 First, search for the measured variable in the chem_vars 
-!--                 list, in order to get the internal name of the variable. 
+!--                 First, search for the measured variable in the chem_vars
+!--                 list, in order to get the internal name of the variable.
                     DO  nn = 1, UBOUND( chem_vars, 2 )
-                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
+                       IF ( TRIM( vmea(l)%var_atts(n)%name ) ==                &
                             TRIM( chem_vars(0,nn) ) )  ind_chem = nn
                     ENDDO
 !
-!--                 Run loop over all chemical species, if the measured 
+!--                 Run loop over all chemical species, if the measured
 !--                 variable matches the interal name, sample the variable.
-                    DO  nn = 1, nvar                    
+!--                 Note, nvar as a chemistry-module variable.
+                    DO  nn = 1, nvar
                        IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
-                            TRIM( chem_species(nn)%name ) )  THEN                            
-                          DO  m = 1, vmea(l)%ns             
+                            TRIM( chem_species(nn)%name ) )  THEN
+                          DO  m = 1, vmea(l)%ns
                              k = vmea(l)%k(m)
                              j = vmea(l)%j(m)
-                             i = vmea(l)%i(m)                    
+                             i = vmea(l)%i(m)
                              vmea(l)%measured_vars(m,n) =                      &
                                                    chem_species(nn)%conc(k,j,i)
@@ -1340,15 +2823,15 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'us' )
+
+              CASE ( 'us' ) ! friction velocity
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not 
-!--                 on ghost points. Hence, limit the indices. 
-                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                    j = MERGE( j           , nyn, j            > nyn )
-                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                    i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
@@ -1364,15 +2847,15 @@
                     ENDDO
                  ENDDO
-                 
-              CASE ( 'ts' )
+
+              CASE ( 'thetas' ) ! scaling parameter temperature
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not 
-!--                 on ghost points. Hence, limit the indices. 
-                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                    j = MERGE( j           , nyn, j            > nyn )
-                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                    i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
@@ -1388,15 +2871,15 @@
                     ENDDO
                  ENDDO
-                 
-              CASE ( 'hfls' )
+
+              CASE ( 'hfls' ) ! surface latent heat flux
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not 
-!--                 on ghost points. Hence, limit the indices. 
-                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                    j = MERGE( j           , nyn, j            > nyn )
-                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                    i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
@@ -1412,15 +2895,15 @@
                     ENDDO
                  ENDDO
-                 
-              CASE ( 'hfss' )
+
+              CASE ( 'hfss' ) ! surface sensible heat flux
                  DO  m = 1, vmea(l)%ns
 !
-!--                 Surface data is only available on inner subdomains, not 
-!--                 on ghost points. Hence, limit the indices. 
-                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                    j = MERGE( j           , nyn, j            > nyn )
-                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                    i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
@@ -1436,16 +2919,48 @@
                     ENDDO
                  ENDDO
-                 
-              CASE ( 'rnds' )
+
+              CASE ( 'hfdg' ) ! ground heat flux
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
+                    DO  mm = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ghf(mm)
+                    ENDDO
+                 ENDDO
+
+              CASE ( 'lwcs' )  ! liquid water of soil layer
+!                  DO  m = 1, vmea(l)%ns
+! !
+! !--                 Surface data is only available on inner subdomains, not
+! !--                 on ghost points. Hence, limit the indices.
+!                     j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+!                     j = MERGE( j           , nyn, j            < nyn )
+!                     i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+!                     i = MERGE( i           , nxr, i            < nxr )
+! 
+!                     DO  mm = surf_lsm_h%start_index(j,i),                      &
+!                              surf_lsm_h%end_index(j,i)
+!                        vmea(l)%measured_vars(m,n) = ?
+!                     ENDDO
+!                  ENDDO
+
+              CASE ( 'rnds' ) ! surface net radiation
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not 
-!--                    on ghost points. Hence, limit the indices. 
-                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                       j = MERGE( j           , nyn, j            > nyn )
-                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                       i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                    Surface data is only available on inner subdomains, not
+!--                    on ghost points. Hence, limit the indices.
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                       j = MERGE( j           , nyn, j            < nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                       i = MERGE( i           , nxr, i            < nxr )
+
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
@@ -1458,16 +2973,16 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'rsus' )
+
+              CASE ( 'rsus' ) ! surface shortwave out
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not 
-!--                    on ghost points. Hence, limit the indices. 
-                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                       j = MERGE( j           , nyn, j            > nyn )
-                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                       i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                    Surface data is only available on inner subdomains, not
+!--                    on ghost points. Hence, limit the indices.
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                       j = MERGE( j           , nyn, j            < nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                       i = MERGE( i           , nxr, i            < nxr )
+
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
@@ -1480,16 +2995,16 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'rsds' )
+
+              CASE ( 'rsds' ) ! surface shortwave in
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not 
-!--                    on ghost points. Hence, limit the indices. 
-                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                       j = MERGE( j           , nyn, j            > nyn )
-                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                       i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                    Surface data is only available on inner subdomains, not
+!--                    on ghost points. Hence, limit the indices.
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                       j = MERGE( j           , nyn, j            < nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                       i = MERGE( i           , nxr, i            < nxr )
+
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
@@ -1502,16 +3017,16 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'rlus' )
+
+              CASE ( 'rlus' ) ! surface longwave out
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not 
-!--                    on ghost points. Hence, limit the indices. 
-                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                       j = MERGE( j           , nyn, j            > nyn )
-                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                       i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                    Surface data is only available on inner subdomains, not
+!--                    on ghost points. Hence, limit the indices.
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                       j = MERGE( j           , nyn, j            < nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                       i = MERGE( i           , nxr, i            < nxr )
+
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
@@ -1524,16 +3039,16 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'rlds' )
+
+              CASE ( 'rlds' ) ! surface longwave in
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
 !
-!--                    Surface data is only available on inner subdomains, not 
-!--                    on ghost points. Hence, limit the indices. 
-                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
-                       j = MERGE( j           , nyn, j            > nyn )
-                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
-                       i = MERGE( i           , nxr, i            > nxr )
-                    
+!--                    Surface data is only available on inner subdomains, not
+!--                    on ghost points. Hence, limit the indices.
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                       j = MERGE( j           , nyn, j            < nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                       i = MERGE( i           , nxr, i            < nxr )
+
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
@@ -1546,75 +3061,115 @@
                     ENDDO
                  ENDIF
-                 
-              CASE ( 'rsd' )
+
+              CASE ( 'rsd' ) ! shortwave in
                  IF ( radiation )  THEN
-                    DO  m = 1, vmea(l)%ns
-                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
-                       j = vmea(l)%j(m)
-                       i = vmea(l)%i(m)
-                    
-                       vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
-                    ENDDO
+                    IF ( radiation_scheme /= 'rrtmg' )  THEN
+                       DO  m = 1, vmea(l)%ns
+                          k = 0
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
+                       ENDDO
+                    ELSE
+                       DO  m = 1, vmea(l)%ns
+                          k = vmea(l)%k(m)
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
+                       ENDDO
+                    ENDIF
                  ENDIF
-                 
-              CASE ( 'rsu' )
+
+              CASE ( 'rsu' ) ! shortwave out
                  IF ( radiation )  THEN
-                    DO  m = 1, vmea(l)%ns
-                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
-                       j = vmea(l)%j(m)
-                       i = vmea(l)%i(m)
-                    
-                       vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
-                    ENDDO
+                    IF ( radiation_scheme /= 'rrtmg' )  THEN
+                       DO  m = 1, vmea(l)%ns
+                          k = 0
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
+                       ENDDO
+                    ELSE
+                       DO  m = 1, vmea(l)%ns
+                          k = vmea(l)%k(m)
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
+                       ENDDO
+                    ENDIF
                  ENDIF
-                 
-              CASE ( 'rlu' )
+
+              CASE ( 'rlu' ) ! longwave out
                  IF ( radiation )  THEN
-                    DO  m = 1, vmea(l)%ns
-                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
-                       j = vmea(l)%j(m)
-                       i = vmea(l)%i(m)
-                    
-                       vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
-                    ENDDO
+                    IF ( radiation_scheme /= 'rrtmg' )  THEN
+                       DO  m = 1, vmea(l)%ns
+                          k = 0
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
+                       ENDDO
+                    ELSE
+                       DO  m = 1, vmea(l)%ns
+                          k = vmea(l)%k(m)
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
+                       ENDDO
+                    ENDIF
                  ENDIF
-                 
-              CASE ( 'rld' )
+
+              CASE ( 'rld' ) ! longwave in
                  IF ( radiation )  THEN
-                    DO  m = 1, vmea(l)%ns
-                       k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' ) 
-                       j = vmea(l)%j(m)
-                       i = vmea(l)%i(m)
-                    
-                       vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
-                    ENDDO
+                    IF ( radiation_scheme /= 'rrtmg' )  THEN
+                       DO  m = 1, vmea(l)%ns
+                          k = 0
+!
+!--                       Surface data is only available on inner subdomains,
+!--                       not on ghost points. Hence, limit the indices.
+                          j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                          j = MERGE( j           , nyn, j            < nyn )
+                          i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                          i = MERGE( i           , nxr, i            < nxr )
+
+                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
+                       ENDDO
+                    ELSE
+                       DO  m = 1, vmea(l)%ns
+                          k = vmea(l)%k(m)
+                          j = vmea(l)%j(m)
+                          i = vmea(l)%i(m)
+                          vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
+                       ENDDO
+                    ENDIF
                  ENDIF
-                 
-              CASE ( 'rsddif' )
+
+              CASE ( 'rsddif' ) ! shortwave in, diffuse part
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                    
+
                        vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
                     ENDDO
                  ENDIF
-                 
-              CASE ( 't_soil' )
+
+              CASE ( 't_soil' ) ! soil and wall temperature
                  DO  m = 1, vmea(l)%ns_soil
-                    i = vmea(l)%i_soil(m)
-                    j = vmea(l)%j_soil(m)
+                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
+                    j = MERGE( j                , nyn, j                 < nyn )
+                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
+                    i = MERGE( i                , nxr, i                 < nxr )
                     k = vmea(l)%k_soil(m)
-                    
+
                     match_lsm = surf_lsm_h%start_index(j,i) <=                 &
                                 surf_lsm_h%end_index(j,i)
                     match_usm = surf_usm_h%start_index(j,i) <=                 &
                                 surf_usm_h%end_index(j,i)
-                                
+
                     IF ( match_lsm )  THEN
                        mm = surf_lsm_h%start_index(j,i)
                        vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
                     ENDIF
-                    
+
                     IF ( match_usm )  THEN
                        mm = surf_usm_h%start_index(j,i)
@@ -1622,19 +3177,84 @@
                     ENDIF
                  ENDDO
-                 
-              CASE ( 'm_soil' )
+
+              CASE ( 'm_soil' ) ! soil moisture
                  DO  m = 1, vmea(l)%ns_soil
-                    i = vmea(l)%i_soil(m)
-                    j = vmea(l)%j_soil(m)
+                    j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
+                    j = MERGE( j                , nyn, j                 < nyn )
+                    i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
+                    i = MERGE( i                , nxr, i                 < nxr )
                     k = vmea(l)%k_soil(m)
-                    
+
                     match_lsm = surf_lsm_h%start_index(j,i) <=                 &
                                 surf_lsm_h%end_index(j,i)
-                                
+
                     IF ( match_lsm )  THEN
                        mm = surf_lsm_h%start_index(j,i)
                        vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
                     ENDIF
-                    
+
+                 ENDDO
+
+              CASE ( 'ts' ) ! surface temperature
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
+                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
+                             surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
+                    ENDDO
+                    DO  mm = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%pt_surface(mm)
+                    ENDDO
+                    DO  mm = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%pt_surface(mm)
+                    ENDDO
+                 ENDDO
+
+              CASE ( 'lwp' ) ! liquid water path
+                 IF ( ASSOCIATED( ql ) )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+
+                       vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i)       &
+                                                       * dzw(1:nzt+1) )        &
+                                                  * rho_surface
+                    ENDDO
+                 ENDIF
+
+              CASE ( 'ps' ) ! surface pressure
+                 vmea(l)%measured_vars(:,n) = surface_pressure
+
+              CASE ( 'pswrtg' ) ! platform speed above ground
+                 vmea(l)%measured_vars(:,n) = 0.0_wp
+
+              CASE ( 'pswrta' ) ! platform speed in air
+                 vmea(l)%measured_vars(:,n) = 0.0_wp
+
+              CASE ( 't_lw' ) ! water temperature
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not
+!--                 on ghost points. Hence, limit the indices.
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
+                    j = MERGE( j           , nyn, j            < nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
+                    i = MERGE( i           , nxr, i            < nxr )
+
+                    DO  mm = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                       IF ( surf_lsm_h%water_surface(m) )                      &
+                          vmea(l)%measured_vars(m,n) = t_soil_h%var_2d(nzt,m)
+                    ENDDO
+
                  ENDDO
 !
@@ -1650,7 +3270,7 @@
 
      ENDDO
-          
+
   END SUBROUTINE vm_sampling
- 
+
 
  END MODULE virtual_measurement_mod
