Index: palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- palm/trunk/SCRIPTS/.palm.iofiles	(revision 3703)
+++ palm/trunk/SCRIPTS/.palm.iofiles	(revision 3704)
@@ -60,4 +60,6 @@
 SURFACE_DATA_AV_BIN*       out:tr   *         $base_data/$run_identifier/OUTPUT         _av_surf   bin
 #
+VIRTUAL_MEAS_BIN*          out   *            $base_data_output/$run_identifier/OUTPUT  _vmeas     bin
+#
 DATA_1D_PR_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _pr        nc
 DATA_1D_SP_NETCDF*         out:tr   *         $base_data/$run_identifier/OUTPUT         _sp        nc
Index: palm/trunk/SOURCE/check_open.f90
===================================================================
--- palm/trunk/SOURCE/check_open.f90	(revision 3703)
+++ palm/trunk/SOURCE/check_open.f90	(revision 3704)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Open binary files for virtual measurements
 ! 
 ! Former revisions:
@@ -470,5 +470,10 @@
           OPEN ( 26, FILE='SURFACE_DATA_AV_BIN'//TRIM( coupling_char )//myid_char,  &
                  FORM='UNFORMATTED', POSITION='APPEND' )
-       
+                 
+       CASE ( 27 )
+!
+!--       Binary files for virtual measurement data
+          OPEN ( 27, FILE='VIRTUAL_MEAS_BIN'//TRIM( coupling_char )//myid_char,&
+                 FORM='UNFORMATTED', POSITION='APPEND' )
 
        CASE ( 30 )
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 3703)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 3704)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Add last_actions for virtual measurements
 ! 
 ! Former revisions:
@@ -290,6 +290,7 @@
 
    USE virtual_measurement_mod,                                                &
-       ONLY:  vm_parin,                                                        &
-              vm_init
+       ONLY:  vm_init,                                                         &
+              vm_last_actions,                                                 &
+              vm_parin
 
    USE wind_turbine_model_mod,                                                 &
@@ -1227,5 +1228,5 @@
 SUBROUTINE module_interface_last_actions
 
-
+   IF ( virtual_measurement )  CALL vm_last_actions
    IF ( user_module_enabled )  CALL user_last_actions
 
Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3703)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 3704)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Formatting adjustments
 ! 
 ! Former revisions:
@@ -623,60 +623,64 @@
        ENDIF
 #else
-    u_ref  = u_ref_l
-    v_ref  = v_ref_l
-    IF ( humidity )       q_ref  = q_ref_l
-    IF ( .NOT. neutral )  pt_ref = pt_ref_l
+       u_ref  = u_ref_l
+       v_ref  = v_ref_l
+       IF ( humidity )       q_ref  = q_ref_l
+       IF ( .NOT. neutral )  pt_ref = pt_ref_l
 #endif
 !
-!-- Average data. Note, reference profiles up to nzt are derived from lateral
-!-- boundaries, at the model top it is derived from the top boundary. Thus,
-!-- number of input data is different from nzb:nzt compared to nzt+1. 
-!-- Derived from lateral boundaries. 
-    u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),      &
-                                            KIND = wp )  
-    v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),      &
-                                            KIND = wp )
-    IF ( humidity )                                                            &
-       q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
-                                                 KIND = wp )
-    IF ( .NOT. neutral )                                                       &
-       pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
-                                           KIND = wp )
-!
-!-- Derived from top boundary.    
-    u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
-    v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
-    IF ( humidity )                                                            &
-       q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
-    IF ( .NOT. neutral )                                                       &
-       pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
-!
-!-- Write onto init profiles, which are used for damping
-    u_init = u_ref
-    v_init = v_ref
-    IF ( humidity      )  q_init  = q_ref
-    IF ( .NOT. neutral )  pt_init = pt_ref
-!
-!-- Set bottom boundary condition
-    IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
-    IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
-!
-!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
-!-- Therefore, calculate boundary-layer depth first. 
-    CALL calc_zi
-    CALL adjust_sponge_layer 
-    
-!
-!-- Update geostrophic wind components from dynamic input file. 
-    DO  k = nzb+1, nzt
-       ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),      &
-                                    fac_dt )
-       vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),      &
-                                    fac_dt )
-    ENDDO
-    ug(nzt+1) = ug(nzt)
-    vg(nzt+1) = vg(nzt)
-    
-    CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
+!--    Average data. Note, reference profiles up to nzt are derived from lateral
+!--    boundaries, at the model top it is derived from the top boundary. Thus,
+!--    number of input data is different from nzb:nzt compared to nzt+1. 
+!--    Derived from lateral boundaries. 
+       u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx     ),   &
+                                               KIND = wp )  
+       v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny   + nx + 1   ),   &
+                                               KIND = wp )
+       IF ( humidity )                                                         &
+          q_ref(nzb:nzt) = q_ref(nzb:nzt)   / REAL( 2.0_wp *                   &
+                                                          ( ny + 1 + nx + 1 ), &
+                                                    KIND = wp )
+       IF ( .NOT. neutral )                                                    &
+          pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp *                   &
+                                                          ( ny + 1 + nx + 1 ), &
+                                              KIND = wp )
+!
+!--    Derived from top boundary.    
+       u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx     ), KIND = wp ) 
+       v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny     ) * ( nx + 1 ), KIND = wp )
+       IF ( humidity )                                                         &
+          q_ref(nzt+1) = q_ref(nzt+1)   / REAL( ( ny + 1 ) * ( nx + 1 ),       &
+                                                KIND = wp )
+       IF ( .NOT. neutral )                                                    &
+          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
+                                                KIND = wp )
+!
+!--    Write onto init profiles, which are used for damping
+       u_init = u_ref
+       v_init = v_ref
+       IF ( humidity      )  q_init  = q_ref
+       IF ( .NOT. neutral )  pt_init = pt_ref
+!
+!--    Set bottom boundary condition
+       IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
+       IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
+!
+!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
+!--    Therefore, calculate boundary-layer depth first. 
+       CALL calc_zi
+       CALL adjust_sponge_layer 
+    
+!
+!--    Update geostrophic wind components from dynamic input file. 
+       DO  k = nzb+1, nzt
+          ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k),   &
+                                       fac_dt )
+          vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k),   &
+                                       fac_dt )
+       ENDDO
+       ug(nzt+1) = ug(nzt)
+       vg(nzt+1) = vg(nzt)
+    
+       CALL  cpu_log( log_point(58), 'offline nesting', 'stop' )
 
     END SUBROUTINE nesting_offl_bc
Index: palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3703)
+++ palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3704)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Interface for attribute input of 8-bit and 32-bit integer 
 ! 
 ! Former revisions:
@@ -777,5 +777,6 @@
     
     INTERFACE netcdf_data_input_att
-       MODULE PROCEDURE netcdf_data_input_att_int
+       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
@@ -1174,9 +1175,54 @@
 ! Description:
 ! ------------
-!> Read a global integer attribute 
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_att_int( val, search_string, id_mod,          &
-                                          input_file, global, openclose,       &
-                                          variable_name )
+!> 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
@@ -1214,5 +1260,5 @@
 #endif            
 
-    END SUBROUTINE netcdf_data_input_att_int
+    END SUBROUTINE netcdf_data_input_att_int32
     
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3703)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3704)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Make variables that are sampled in virtual measurement module public 
 ! 
 ! Former revisions:
@@ -1276,5 +1276,6 @@
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
-           skyvf, skyvft, radiation_interactions_on, average_radiation
+           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
+           rad_sw_in_diff, rad_sw_in_dir
 
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3703)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3704)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Data output for virtual measurements added
 ! 
 ! Former revisions:
@@ -632,5 +632,5 @@
                
     USE virtual_measurement_mod,                                               &
-        ONLY:  vm_sampling, vm_time_start
+        ONLY:  vm_data_output, vm_sampling, vm_time_start
 
     IMPLICIT NONE
@@ -1612,5 +1612,8 @@
 !--    Take virtual measurements
        IF ( virtual_measurement  .AND.                                         &
-            vm_time_start <= time_since_reference_point )  CALL vm_sampling
+            vm_time_start <= time_since_reference_point )  THEN
+          CALL vm_sampling
+          CALL vm_data_output
+       ENDIF
 
 !
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3703)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3704)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! make nzb_wall public, required for virtual-measurements
 ! 
 ! Former revisions:
@@ -1150,6 +1150,6 @@
 !-- Public parameters, constants and initial values
     PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod,           &
-           usm_green_heat_model, building_pars,  &
-           nzt_wall, t_wall_h, t_wall_v,         &
+           usm_green_heat_model, building_pars,                                &
+           nzb_wall, nzt_wall, t_wall_h, t_wall_v,                             &
            t_window_h, t_window_v, building_type
 
@@ -7469,5 +7469,5 @@
         IMPLICIT NONE
 
-        INTEGER(iwp)                          :: i, j, k, l, m   !< running indices
+        INTEGER(iwp) :: i, j, k, l, m   !< running indices
         
         INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3703)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3704)
@@ -20,12 +20,11 @@
 ! Current revisions:
 ! -----------------
-! 
+! - initialization revised
+! - binary data output
+! - list of allowed variables extended
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
-! unused variables removed
-! 
-! 3522 2018-11-13 12:14:36Z suehring
 ! Sampling of variables
 ! 
@@ -44,5 +43,5 @@
 !> 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' measurement positions. 
+!> 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, 
@@ -53,6 +52,4 @@
 !>
 !> @todo list_of_allowed variables needs careful checking
-!> @todo output (binary or NetCDF) needs to be implemented
-!> @todo clean-up anything from current test modus
 !> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
 !>       the (UC)2 standard 
@@ -71,6 +68,6 @@
         
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, dz, humidity, neutral, message_string,           &
-               virtual_measurement
+        ONLY:  air_chemistry, dz, humidity, io_blocks, io_group, neutral,      &
+               message_string, time_since_reference_point, virtual_measurement
 
     USE cpulog,                                                                &
@@ -81,40 +78,75 @@
 
     USE indices,                                                               &
-        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny
+        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0
 
     USE kinds
+    
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  init_model
+        
+    USE pegrid
+    
+    USE surface_mod,                                                           &
+        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
+        
+    USE urban_surface_mod,                                                     &
+        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_mea
     
-       CHARACTER(LEN=100)  ::  feature_type  !< type of the measurement
-       CHARACTER(LEN=100)  ::  site          !< name of the measurement site
+       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 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
-       INTEGER(iwp) ::  ntraj  !< number of trajectories of a measurement
-       INTEGER(iwp) ::  nvar   !< number of 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) ::  nvar            !< 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
        
-       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       !< 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) ::  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       !< measured variables 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars_soil  !< measured variables 
        
     END TYPE virt_mea
@@ -122,4 +154,8 @@
     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
+    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
@@ -133,7 +169,25 @@
     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',                                          &
+                            'm_soil',                                          &
+                            'lwc   ',                                          &
+                            'lwcs  ',                                          &
+                            'smp   '                       /)
+                            
+    CHARACTER(LEN=10), DIMENSION(0:1,1:8) ::  chem_vars = RESHAPE( (/          &
+                                              'mcpm1     ', 'kc_PM1    ',      &
+                                              'mcpm2p5   ', 'kc_PM2.5  ',      &
+                                              'mcpm25    ', 'kc_PM25   ',      &
+                                              'mcpm10    ', 'kc_PM10   ',      &
+                                              'mfno2     ', 'kc_NO2    ',      &
+                                              'mfno      ', 'kc_NO     ',      &
+                                              'tro3      ', 'kc_O3     ',      &
+                                              'mfco      ', 'kc_CO     '       &
+                                                                   /), (/ 2, 8 /) )
 !
 !-- MS: List requires careful revision!
-    CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
+    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)
@@ -143,5 +197,4 @@
           'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
           'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
-          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
           'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
           'mcco      ',  & ! mass concentration of CO (kg/m3)
@@ -173,8 +226,11 @@
           '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)
@@ -182,20 +238,21 @@
           'rsu       ',  & ! upward shortwave radiative flux (W/m2)
           'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
-          'rnds      '   & ! surface net downward radiative flux (W/m2)
+          'rnds      ',  & ! surface net downward radiative flux (W/m2)
+          't_soil    ',  &
+          'm_soil    ',  &
+          'lwc       ',  &
+          'lwcs      ',  &
+          'smp       '   &
        /)
-                                                             
-    INTEGER(iwp) ::  id_vm                           !< NetCDF file id for virtual measurements
-    INTEGER(iwp) ::  nvm = 0                         !< number of virtual measurements
-!    INTEGER(iwp) ::  observation_coverage_xy = 0     !< horizontal distance from the measurement point where observations should be taken in the surrounding
-!    INTEGER(iwp) ::  observation_coverage_z  = 0     !< vertical distance from the measurement point where observations should be taken in the surrounding
-    
+                                                            
+    
+    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
-    LOGICAL ::  global_attribute = .TRUE.         !< flag indicating a global attribute
     
     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
-    
     
     INTERFACE vm_check_parameters
@@ -203,8 +260,16 @@
     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_parin
        MODULE PROCEDURE vm_parin
@@ -221,9 +286,10 @@
 !
 !-- Public interfaces
-    PUBLIC  vm_check_parameters, vm_init, vm_parin, vm_sampling
+    PUBLIC  vm_check_parameters, vm_data_output, vm_init, vm_last_actions,     &
+            vm_parin, vm_sampling
 
 !
 !-- Public variables
-    PUBLIC  vmea, vm_time_start
+    PUBLIC  vmea, vmea_general, vm_time_start
 
  CONTAINS
@@ -251,5 +317,5 @@
 !-- ToDo: Revise this later and remove this requirement.
     IF ( virtual_measurement  .AND.  .NOT. input_pids_static )  THEN
-       message_string = 'If virtual measurements are taken a static input ' // &
+       message_string = 'If virtual measurements are taken, a static input ' //&
                         'file is mandatory.'
        CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
@@ -319,5 +385,6 @@
  
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  input_file_vm, netcdf_data_input_get_dimension_length,          &
+        ONLY:  init_model, input_file_vm,                                      &
+               netcdf_data_input_get_dimension_length,                         &
                netcdf_data_input_att, netcdf_data_input_var
                
@@ -328,17 +395,18 @@
     
     CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
+    CHARACTER(LEN=5)    ::  dummy_read                !< 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_eutm  !< dimension size of UTM easting coordinate
-!    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
+    INTEGER(iwp) ::  dim_eutm  !< dimension size of UTM easting coordinate
+    INTEGER(iwp) ::  dim_nutm  !< dimension size of UTM northing coordinate
     INTEGER(iwp) ::  dim_ntime !< dimension size of time coordinate 
-!    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
-!    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
-!    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
+    INTEGER(iwp) ::  dim_zag   !< dimension size of height coordinate
+    INTEGER(iwp) ::  i         !< grid index of virtual observation point in x-direction
+    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
     INTEGER(iwp) ::  is        !< grid index of real observation point of the respective station in x-direction
-!    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
+    INTEGER(iwp) ::  j         !< grid index of observation point in x-direction
     INTEGER(iwp) ::  js        !< grid index of real observation point of the respective station in y-direction
-!    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
+    INTEGER(iwp) ::  k         !< grid index of observation point in x-direction
     INTEGER(iwp) ::  kl        !< lower vertical index of surrounding grid points of an observation coordinate
     INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
@@ -352,4 +420,9 @@
     INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
     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(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
@@ -366,24 +439,28 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  z_ag  !< height coordinate relative to origin_z, temporary variable
 !
-!-- Obtain number of virtual measurement stations
-    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm,   &
+!-- 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', '' )
                                 
-!     write(9,*) "num stationi", nvm
-!     flush(9)
-!
-!-- ALLOCATE data structure
-    ALLOCATE( vmea(1:nvm) )
-!
-!-- Allocate flag array
+!
+!-- Allocate data structure which encompass all required information, such as
+!-- 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) )
     meas_flag = 0
- 
-!
-!-- Read station coordinates and further attributes. 
-!-- Note all coordinates are in UTM coordinates.
-    DO  l = 1, nvm
-!
-!--    Determine suffix which contains the ID
+!
+!-- Loop over all sites.
+    DO  l = 1, vmea_general%nvm
+!
+!--    Determine suffix which contains the ID, ordered according to the number
+!--    of measurements. 
        IF( l < 10 )  THEN
           WRITE( dum, '(I1)')  l
@@ -397,18 +474,36 @@
           WRITE( dum, '(I5)')  l
        ENDIF
-       
-       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx            &
-                                   // TRIM( dum ), id_vm, '', global_attribute,&
+!
+!--    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,   &
                                    '', '' )
-       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy            &
-                                   // TRIM( dum ), id_vm, '', global_attribute,&
+!
+!--    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,   &
                                    '', '' )
-       CALL netcdf_data_input_att( vmea(l)%site,         char_site             &
-                                   // TRIM( dum ), id_vm, '', global_attribute,&
-                                   '', '' )
-       CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature          &
-                                   // TRIM( dum ), id_vm, '', global_attribute,&
-                                   '', '' )
-       
+!
+!--    Set flag for soil-sampling.
+       IF ( soil_dum == 1 )  vmea(l)%soil_sampling = .TRUE.
 !
 !---   Set logicals depending on the type of the measurement
@@ -419,7 +514,7 @@
        ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
           vmea(l)%trajectory        = .TRUE.
+!
+!--   Give error message in case the type matches non of the pre-defined types.
        ELSE
-!
-!--       Give error message
           message_string = 'Attribue featureType = ' //                        &
                            TRIM( vmea(l)%feature_type ) //                     &
@@ -428,16 +523,16 @@
        ENDIF
 !
-!--    Read string with all measured variables at this station
+!--    Read string with all measured variables at this site
        measured_variables_file = ''
        CALL netcdf_data_input_var( measured_variables_file,                    &
-                                   char_mv // TRIM( dum ), id_vm )
-!
-!--    Count the number of measured variables which match with the variables
-!--    which are allowed to be measured in PALM. 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)%nvar = 0
+                                   char_mv // TRIM( dum ), vmea_general%id_vm )
+!
+!--    Count the number of measured variables. Only count variables that match
+!--    with the allowed 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)%nvar   = 0
        DO ll = 1, SIZE( measured_variables_file )
           IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
@@ -465,5 +560,5 @@
        ENDDO
 !
-!--    Allocate array for the measured variables names for the station l.
+!--    Allocate array for the measured variables names for the respective site.
        ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
 
@@ -474,52 +569,56 @@
 !--    In case of chemistry, check if species is considered in the modelled
 !--    chemistry mechanism.
-       IF ( air_chemistry )  THEN
-          DO  ll = 1, vmea(l)%nvar
-             chem_include = .FALSE.
-             DO  n = 1, nspec
-                IF ( TRIM( vmea(l)%measured_vars_name(ll) ) ==                 &
-                     TRIM( chem_species(n)%name ) )  chem_include = .TRUE.
-             ENDDO
-             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
-!
-!--    For the actual measurement ID read the UTM coordinates. Based on these,
-!--    define the index space on each subdomain where measurements should be
-!--    taken. Note, the entire coordinate arrays will not be stored on data
-!--    type as this would exceed memory requirements, particularly for 
-!--    trajectory measurements. If no variable will be virtually measured, 
-!--    skip the reading. 
+!        IF ( air_chemistry )  THEN
+!           DO  ll = 1, vmea(l)%nvar
+!              chem_include = .FALSE.
+!              DO  n = 1, nspec
+!                 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,
+!--    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. 
        IF ( vmea(l)%nvar > 0 )  THEN
 !
 !--       For stationary measurements UTM coordinates are just one value and 
 !--       its dimension is "station", while for mobile measurements UTM 
-!--       coordinates are arrays. First, inquire dimension length for 
-!--       UTM coordinates.
+!--       coordinates are arrays depending on the number of trajectories and
+!--       time, according to (UC)2 standard. First, inquire dimension length 
+!--       of the UTM coordinates.
           IF ( vmea(l)%trajectory )  THEN
 !
 !--          For non-stationary measurements read the number of trajectories
-             CALL netcdf_data_input_get_dimension_length( id_vm,              &
+!--          and the number of time coordinates.
+             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
                                                           vmea(l)%ntraj,      &
                                                           "traj" //           &
                                                           TRIM( dum ) )
-             CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime,   &
+             CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
+                                                          dim_ntime,          &
                                                           "ntime" //          &
                                                           TRIM( dum ) )
 !
-!--       For stationary measurements the dimension for UTM coordinates is 1
+!--       For stationary measurements the dimension for UTM and time 
+!--       coordinates is 1. 
           ELSE
              vmea(l)%ntraj  = 1
              dim_ntime = 1
           ENDIF
-          
 !
 !-        Allocate array which defines individual time frame for each 
-!--       trajectory or station
+!--       trajectory or station.
           ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
 !
@@ -530,13 +629,16 @@
           ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
 !
-!--       Read _FillValue attributes
+!--       Read _FillValue attributes of the coordinate dimensions.
           CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
-                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      vmea_general%id_vm, '',                  &
+                                      .NOT. global_attribute, '',              &
                                       char_eutm // TRIM( dum ) )
           CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
-                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      vmea_general%id_vm, '',                  &
+                                      .NOT. global_attribute, '',              &
                                       char_nutm // TRIM( dum ) )
           CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
-                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      vmea_general%id_vm, '',                  &
+                                      .NOT. global_attribute, '',              &
                                       char_zag  // TRIM( dum ) )
 !
@@ -544,29 +646,21 @@
 !--       times. 
           IF ( vmea(l)%trajectory )  THEN
-             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm,  &
+             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 ), id_vm,  &
+             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 ), id_vm,  &
+             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ),        &
+                                         vmea_general%id_vm,                   &
                                          0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
           ELSE
-             CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
-             CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
-             CALL netcdf_data_input_var( z_ag(1,:),  char_zag  // TRIM( dum ), id_vm )
-          ENDIF
-!
-!-- For testing:
-          e_utm = e_utm - e_utm(1,1)
-          n_utm = n_utm - n_utm(1,1)
-          
-!             
-!--       First, compute relative x- and y-coordinates with respect to the 
-!--       lower-left origin of the model domain, which is the difference 
-!--       betwen UTM coordinates.  
-!           e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t))            &
-!                                       - init_model%origin_x
-!           n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t))            &
-!                                       - init_model%origin_y
-
+             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         
 !
 !--       Based on UTM coordinates, check if the measurement station or parts
@@ -575,4 +669,11 @@
           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 
+!--          betwen UTM coordinates. Note, if the origin is not correct, the
+!--          virtual sites will be misplaced. 
+             e_utm(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
+             n_utm(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
 !
 !--          Determine the individual time coordinate length for each station and
@@ -580,5 +681,6 @@
 !--          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 
-!--          processed further.
+!--          processed further. This is actually a work-around for the specific
+!--          (UC)2 dataset, but it won't harm in anyway. 
              vmea(l)%dim_t(t) = 0
              DO  n = 1, dim_ntime
@@ -587,5 +689,4 @@
                      z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
              ENDDO
-
 !
 !--          Compute grid indices relative to origin and check if these are
@@ -612,9 +713,18 @@
 !--                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 )
-                 
-                   meas_flag(kl:ku,js-1:js+1,is-1:is+1) =                      &
-                               IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 ) 
+                   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
+                         DO  k = kl, ku
+                            meas_flag(k,j,i) = MERGE(                          &
+                                             IBSET( meas_flag(k,j,i), 0 ),     &
+                                             0,                                &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                                    )
+                         ENDDO
+                      ENDDO
+                   ENDDO
                 ENDIF
              ENDDO
@@ -622,54 +732,121 @@
           ENDDO
 !
-!--       Based on the flag array count the number of observation points.
+!--       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.           
+!--       Hence, count individually. Start with atmoshere.
           ns = 0
-          DO  is = nxl-1, nxr+1
-             DO  js = nys-1, nyn+1
-                DO  ks = nzb, nzt+1
-                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                DO  k = nzb, nzt+1
+                   ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
                 ENDDO
              ENDDO
           ENDDO
+          
 !
 !--       Store number of observation points on subdomain and allocate index
-!--       arrays.
+!--       arrays as well as array containing height information.
           vmea(l)%ns = ns
-          ns = 0
           
           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. 
-          DO  is = nxl-1, nxr+1
-             DO  js = nys-1, nyn+1
-                DO  ks = nzb, nzt+1
-                   IF ( BTEST( meas_flag(ks,js,is), 0 ) )  THEN
+          ns = 0
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                DO  k = nzb, nzt+1
+                   IF ( BTEST( meas_flag(k,j,i), 0 ) )  THEN
                       ns = ns + 1
-                      vmea(l)%i(ns) = is
-                      vmea(l)%j(ns) = js
-                      vmea(l)%k(ns) = ks
+                      vmea(l)%i(ns) = i
+                      vmea(l)%j(ns) = j
+                      vmea(l)%k(ns) = k
+                      vmea(l)%z_ag(ns)  = zu(k) -                              &
+                                   zw(get_topography_top_index_ji( j, i, 's' ))
                    ENDIF
                 ENDDO
              ENDDO
           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. 
+          IF ( vmea(l)%soil_sampling )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
+                      IF ( surf_lsm_h%start_index(j,i) <=                      &
+                           surf_lsm_h%end_index(j,i) )  THEN
+                         vmea(l)%ns_soil = vmea(l)%ns_soil +                   &
+                                                      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 
+                      ENDIF
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDIF         
+!
+!--       Allocate index arrays as well as array containing height information 
+!--       for soil.
+          IF ( vmea(l)%soil_sampling )  THEN
+             ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
+             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     
+!
+!--       For soil, store the grid indices.
+          ns = 0
+          IF ( vmea(l)%soil_sampling )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) )  THEN
+                      IF ( surf_lsm_h%start_index(j,i) <=                      &
+                           surf_lsm_h%end_index(j,i) )  THEN
+                         m = surf_lsm_h%start_index(j,i)
+                         DO  k = nzb_soil, nzt_soil
+                            ns = ns + 1
+                            vmea(l)%i_soil(ns) = i
+                            vmea(l)%j_soil(ns) = j
+                            vmea(l)%k_soil(ns) = 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
+                         m = surf_usm_h%start_index(j,i)
+                         DO  k = nzb_wall, nzt_wall
+                            ns = ns + 1
+                            vmea(l)%i_soil(ns) = i
+                            vmea(l)%j_soil(ns) = j
+                            vmea(l)%k_soil(ns) = k
+                            vmea(l)%depth(ns)  = surf_usm_h%zw(k,m)
+                         ENDDO
+                      ENDIF
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDIF
+!
+!--       Allocate array to save the sampled values.
+          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) )
           
-!           write(9,*) l, "NS: ", ns
-!           IF ( ns > 0 )   THEN
-!              write(9,*) "i ", vmea(l)%i
-!              write(9,*) "j ", vmea(l)%j
-!              write(9,*) "k ", vmea(l)%k
-!              write(9,*)
-!           ENDIF
-!
-!--       Allocate array to save the sampled values.
-!--       Todo: Is it better to allocate for all variables at a station
-!--       and store all the values before writing, or sample the variables
-!--       directly in the data output?
-          ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
-!
-!--       Initialize with _FillValue
-          vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
+          IF ( vmea(l)%soil_sampling )                                         &
+             ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,           &
+                                                  1:vmea(l)%nvar) )
+!
+!--       Initialize with _FillValues
+          vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) = vmea(l)%fillout
+          IF ( vmea(l)%soil_sampling )                                         &
+             vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nvar) =    &
+                                                                vmea(l)%fillout
 !
 !--       Deallocate temporary coordinate arrays
@@ -677,16 +854,47 @@
           IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
           IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
+          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( vmea(l)%dim_t )
        ENDIF
     ENDDO
-!     flush(9)
-    
 !
 !-- Close input file for virtual measurements. Therefore, just call 
 !-- the read attribute routine with the "close" option. 
-    CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
+    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. 
+!-- This is actually only required for data output.
+    ALLOCATE( ns_all(1:vmea_general%nvm) )
+    ns_all = 0    
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER,  &
+                        MPI_SUM, comm2d, ierr )
+#else
+    ns_all(:) = vmea(:)%ns
+#endif
+    vmea(:)%ns_tot = ns_all(:)
+!
+!-- Now for soil
+    ns_all = 0    
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm,          &
+                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
+#else
+    ns_all(:) = vmea(:)%ns_soil
+#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
         
   END SUBROUTINE vm_init
@@ -696,4 +904,150 @@
 ! Description:
 ! ------------
+!> Binary data output.
+!------------------------------------------------------------------------------!
+  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)%nvar
+                 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 +                           &
+                        REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx
+                 WRITE ( 27 )  'N_UTM                             '
+                 WRITE ( 27 )  init_model%origin_y +                           &
+                        REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy
+                 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
+                    print*, "vmea(l)%ns_soil_tot", 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
+                    WRITE ( 27 )  'N_UTM soil                        '
+                    WRITE ( 27 )  init_model%origin_y +                        &
+                           REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp,   &
+                                 KIND = wp ) * dy
+                    WRITE ( 27 )  'DEPTH                             '
+                    WRITE ( 27 )  vmea(l)%depth(1:vmea(l)%ns_soil)
+                 ENDIF
+              ENDDO
+
+           ENDIF
+        ENDDO
+        
+#if defined( __parallel )
+        CALL MPI_BARRIER( comm2d, 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)%nvar
+                    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    
+     INTEGER(iwp) ::  l         !< running index over all stations
+     INTEGER(iwp) ::  n         !< running index over all measured variables at a station
+ 
+     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 
+  
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Sampling of the actual quantities along the observation coordinates
 !------------------------------------------------------------------------------!
@@ -714,21 +1068,29 @@
      IMPLICIT NONE
      
-!     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
+     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
       
-     INTEGER(iwp) ::  i  !< grid index in x-direction
-     INTEGER(iwp) ::  j  !< grid index in y-direction
-     INTEGER(iwp) ::  k  !< grid index in z-direction
-     INTEGER(iwp) ::  l  !< running index over the number of stations
-     INTEGER(iwp) ::  m  !< running index over all virtual observation coordinates
-     INTEGER(iwp) ::  mm !< index of surface element which corresponds to the virtual observation coordinate
-     INTEGER(iwp) ::  n  !< running index over all measured variables at a station
-     INTEGER(iwp) ::  nn !< running index over the number of chemcal species
-!
-!--  Loop over all stations. For each possible variable loop over all 
-!--  observation points
-     DO  l = 1, nvm
-!
-!--     Loop over all measured variables at this station. Please note, 
-!--     velocity components are interpolated onto scalar grid.  
+     INTEGER(iwp) ::  i         !< grid index in x-direction
+     INTEGER(iwp) ::  j         !< grid index in y-direction
+     INTEGER(iwp) ::  k         !< grid index in z-direction
+     INTEGER(iwp) ::  ind_chem  !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
+     INTEGER(iwp) ::  l         !< running index over the number of stations
+     INTEGER(iwp) ::  m         !< running index over all virtual observation coordinates
+     INTEGER(iwp) ::  mm        !< index of surface element which corresponds to the virtual observation coordinate
+     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. 
+     DO  l = 1, vmea_general%nvm
+!
+!--     At the beginning, set _FillValues
+        IF ( ALLOCATED( vmea(l)%measured_vars      ) )                         &
+           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.  
         DO  n = 1, vmea(l)%nvar
         
@@ -741,9 +1103,9 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(n,m) = pt(k,j,i)
+                       vmea(l)%measured_vars(m,n) = pt(k,j,i)
                     ENDDO
                  ENDIF
                  
-              CASE ( 'ta', 't_va' )
+              CASE ( 'ta' )
                  IF ( .NOT. neutral )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -751,7 +1113,9 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
+                       vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
                     ENDDO
                  ENDIF
+              
+              CASE ( 't_va' )
                  
               CASE ( 'hus', 'haa' )
@@ -761,5 +1125,5 @@
                        j = vmea(l)%j(m)
                        i = vmea(l)%i(m)
-                       vmea(l)%measured_vars(n,m) = q(k,j,i)
+                       vmea(l)%measured_vars(m,n) = q(k,j,i)
                     ENDDO
                  ENDIF
@@ -770,5 +1134,5 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
                  ENDDO
                  
@@ -778,13 +1142,13 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
                  ENDDO
                  
               CASE ( 'w' )
                  DO  m = 1, vmea(l)%ns
-                    k = 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(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
+                    vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
                  ENDDO
                  
@@ -794,5 +1158,5 @@
                     j = vmea(l)%j(m)
                     i = vmea(l)%i(m)
-                    vmea(l)%measured_vars(n,m) = SQRT(                         &
+                    vmea(l)%measured_vars(m,n) = SQRT(                         &
                                    ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
                                    ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2   &
@@ -806,17 +1170,86 @@
                     i = vmea(l)%i(m)
                     
-                    vmea(l)%measured_vars(n,m) = ATAN2(                        &
+                    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
                  ENDDO
-!
-!--           MS: list of variables needs extension.
-              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda',  &
-                     'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno',        &
-                     'mfno2', 'mfso2', 'mfh20', 'tr03' )
-                 IF ( air_chemistry )  THEN                 
+                 
+              CASE ( 'utheta' )
+                 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) ) *   &
+                                                   pt(k,j,i)
+                 ENDDO
+                 
+              CASE ( 'vtheta' )
+                 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) ) *   &
+                                                   pt(k,j,i)
+                 ENDDO
+                 
+              CASE ( 'wtheta' )
+                 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) ) *   &
+                                                   pt(k,j,i)
+                 ENDDO
+                 
+              CASE ( 'uw' )
+                 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.25_wp *                     &
+                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
+                                                 ( u(k,j,i)   + u(k,j,i+1) )
+                 ENDDO
+                 
+              CASE ( 'vw' )
+                 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.25_wp *                     &
+                                                 ( w(k-1,j,i) + w(k,j,i) ) *   &
+                                                 ( 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) )
+                    j = vmea(l)%j(m)
+                    i = vmea(l)%i(m)
+                    vmea(l)%measured_vars(m,n) = 0.25_wp *                     &
+                                                 ( u(k,j,i)   + u(k,j,i+1) ) * &
+                                                 ( v(k,j,i)   + v(k,j+1,i) )
+                 ENDDO
+!
+!--           List of variables may need extension.
+              CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2',    & 
+                     '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. 
+                    DO  nn = 1, UBOUND( chem_vars, 2 )
+                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
+                            TRIM( chem_vars(0,nn) ) )  ind_chem = nn
+                    ENDDO
+!
+!--                 Run loop over all chemical species, if the measured 
+!--                 variable matches the interal name, sample the variable.
                     DO  nn = 1, nspec                    
-                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
+                       IF ( TRIM( chem_vars(1,ind_chem) ) ==                   &
                             TRIM( chem_species(nn)%name ) )  THEN                            
                           DO  m = 1, vmea(l)%ns             
@@ -824,5 +1257,5 @@
                              j = vmea(l)%j(m)
                              i = vmea(l)%i(m)                    
-                             vmea(l)%measured_vars(n,m) =                      &
+                             vmea(l)%measured_vars(m,n) =                      &
                                                    chem_species(nn)%conc(k,j,i)
                           ENDDO
@@ -843,13 +1276,13 @@
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
                     ENDDO
                     DO  mm = surf_lsm_h%start_index(j,i),                      &
                              surf_lsm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
                     ENDDO
                     DO  mm = surf_usm_h%start_index(j,i),                      &
                              surf_usm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
                     ENDDO
                  ENDDO
@@ -867,13 +1300,13 @@
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
                     ENDDO
                     DO  mm = surf_lsm_h%start_index(j,i),                      &
                              surf_lsm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
                     ENDDO
                     DO  mm = surf_usm_h%start_index(j,i),                      &
                              surf_usm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
                     ENDDO
                  ENDDO
@@ -891,13 +1324,13 @@
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
                     ENDDO
                     DO  mm = surf_lsm_h%start_index(j,i),                      &
                              surf_lsm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
                     ENDDO
                     DO  mm = surf_usm_h%start_index(j,i),                      &
                              surf_usm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
                     ENDDO
                  ENDDO
@@ -915,13 +1348,13 @@
                     DO  mm = surf_def_h(0)%start_index(j,i),                   &
                              surf_def_h(0)%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
+                       vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
                     ENDDO
                     DO  mm = surf_lsm_h%start_index(j,i),                      &
                              surf_lsm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
+                       vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
                     ENDDO
                     DO  mm = surf_usm_h%start_index(j,i),                      &
                              surf_usm_h%end_index(j,i)
-                       vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
+                       vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
                     ENDDO
                  ENDDO
@@ -940,14 +1373,14 @@
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
+                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
                        ENDDO
                        DO  mm = surf_usm_h%start_index(j,i),                   &
                                 surf_usm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
+                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
                        ENDDO
                     ENDDO
                  ENDIF
                  
-              CASE ( 'rsus', 'rsu' )
+              CASE ( 'rsus' )
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -962,14 +1395,14 @@
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
+                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
                        ENDDO
                        DO  mm = surf_usm_h%start_index(j,i),                   &
                                 surf_usm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
+                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
                        ENDDO
                     ENDDO
                  ENDIF
                  
-              CASE ( 'rsds', 'rsd' )
+              CASE ( 'rsds' )
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -984,14 +1417,14 @@
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
+                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
                        ENDDO
                        DO  mm = surf_usm_h%start_index(j,i),                   &
                                 surf_usm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
+                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
                        ENDDO
                     ENDDO
                  ENDIF
                  
-              CASE ( 'rlus', 'rlu' )
+              CASE ( 'rlus' )
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -1006,14 +1439,14 @@
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
+                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
                        ENDDO
                        DO  mm = surf_usm_h%start_index(j,i),                   &
                                 surf_usm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
+                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
                        ENDDO
                     ENDDO
                  ENDIF
                  
-              CASE ( 'rlds', 'rld' )
+              CASE ( 'rlds' )
                  IF ( radiation )  THEN
                     DO  m = 1, vmea(l)%ns
@@ -1028,15 +1461,111 @@
                        DO  mm = surf_lsm_h%start_index(j,i),                   &
                                 surf_lsm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
+                          vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
                        ENDDO
                        DO  mm = surf_usm_h%start_index(j,i),                   &
                                 surf_usm_h%end_index(j,i)
-                          vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
+                          vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
                        ENDDO
                     ENDDO
                  ENDIF
+                 
+              CASE ( 'rsd' )
+                 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
+                 ENDIF
+                 
+              CASE ( 'rsu' )
+                 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
+                 ENDIF
+                 
+              CASE ( 'rlu' )
+                 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
+                 ENDIF
+                 
+              CASE ( 'rld' )
+                 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
+                 ENDIF
+                 
+              CASE ( 'rsddif' )
+                 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' )
+                 DO  m = 1, vmea(l)%ns_soil
+                    i = vmea(l)%i_soil(m)
+                    j = vmea(l)%j_soil(m)
+                    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)
+                       vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
+                    ENDIF
+                 ENDDO
+                 
+              CASE ( 'm_soil' )
+                 DO  m = 1, vmea(l)%ns_soil
+                    i = vmea(l)%i_soil(m)
+                    j = vmea(l)%j_soil(m)
+                    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
 !
 !--           More will follow ...
-                 
+
+!
+!--           No match found - just set a fill value
+              CASE DEFAULT
+                 vmea(l)%measured_vars(:,n) = vmea(l)%fillout
            END SELECT
 
@@ -1044,5 +1573,5 @@
 
      ENDDO
-     
+          
   END SUBROUTINE vm_sampling
  
Index: palm/trunk/UTIL/combine_virtual_measurements/combine_virtual_measurements.f90
===================================================================
--- palm/trunk/UTIL/combine_virtual_measurements/combine_virtual_measurements.f90	(revision 3704)
+++ palm/trunk/UTIL/combine_virtual_measurements/combine_virtual_measurements.f90	(revision 3704)
@@ -0,0 +1,820 @@
+!> @file combine_virtual_measurements.f90
+!------------------------------------------------------------------------------!
+! 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-2018  Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! Initial revsion
+! 
+! Former revisions:
+! -----------------
+! $Id$
+!
+! Authors:
+! --------
+! @author Matthias Suehring
+!
+!
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routines merges binary output from virtual measurements taken from 
+!> different subdomains and creates a NetCDF output file according to the (UC)2
+!> data standard. 
+!------------------------------------------------------------------------------!
+ PROGRAM combine_virtual_measurements
+ 
+#if defined( __netcdf )  
+    USE NETCDF
+#endif
+ 
+    IMPLICIT NONE
+    
+    CHARACTER(LEN=34)  ::  char_in                !< dummy string
+    CHARACTER(LEN=4)   ::  file_suffix = '.bin'   !< string which contain the suffix indicating virtual measurement data  
+    CHARACTER(LEN=30)  ::  myid_char              !< combined string indicating binary file
+    CHARACTER(LEN=100) ::  path                   !< path to the binary data 
+    CHARACTER(LEN=100) ::  run                    !< name of the run
+    
+    CHARACTER(LEN=100), DIMENSION(:),   ALLOCATABLE ::  site          !< name of the site
+    CHARACTER(LEN=100), DIMENSION(:),   ALLOCATABLE ::  filename      !< name of the original file
+    CHARACTER(LEN=100), DIMENSION(:),   ALLOCATABLE ::  feature_type  !< string indicating the type of the measurement
+    CHARACTER(LEN=100), DIMENSION(:),   ALLOCATABLE ::  soil_quantity !< string indicating soil measurements
+    CHARACTER(LEN=10),  DIMENSION(:,:), ALLOCATABLE ::  variables     !< list of measured variables 
+    
+    CHARACTER(LEN=6), DIMENSION(5) ::  soil_quantities = (/                   & !< list of measurable soil variables
+                            "t_soil",                                         &
+                            "m_soil",                                         &
+                            "lwc   ",                                         &
+                            "lwcs  ",                                         &
+                            "smp   "                    /)
+    
+    INTEGER, PARAMETER ::  iwp = 4                !< integer precision
+    INTEGER, PARAMETER ::  wp  = 8                !< float precision
+    
+    INTEGER(iwp) ::  cycle_number                 !< cycle number
+    INTEGER(iwp) ::  f                            !< running index over all binary files
+    INTEGER(iwp) ::  file_id_in = 18              !< file unit for input binaray file    
+    INTEGER(iwp) ::  l                            !< running index indicating the actual site
+    INTEGER(iwp) ::  n                            !< running index over all variables measured at a site
+    INTEGER(iwp) ::  num_pe                       !< number of processors used for the run
+    INTEGER(iwp) ::  nvm                          !< number of sites
+    INTEGER(iwp) ::  status_nc                    !< NetCDF error code, return value
+       
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns           !< number of observation coordinates on current subdomain
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_tot       !< total number of observation coordinates for a site 
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_soil      !< number of observation coordinates for soil quantities (on current subdomain)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ns_soil_tot  !< total number of observation coordinates for a site (for the soil)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvar         !< number of sampled variables at a site 
+!
+!-- NetCDF varialbes
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nc_id                 !< NetCDF file ID
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_eutm               !< NetCDF dimension ID for E_UTM
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_nutm               !< NetCDF dimension ID for N_UTM
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_hao                !< NetCDF dimension ID for the height coordinate
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_eutm_soil          !< NetCDF dimension ID for E_UTM for soil quantity
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_nutm_soil          !< NetCDF dimension ID for N_UTM for soil quantity
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_depth              !< NetCDF dimension ID for soil depth
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_time               !< NetCDF dimension ID for time
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_eutm           !< NetCDF variable ID for E_UTM
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_nutm           !< NetCDF variable ID for N_UTM
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_hao            !< NetCDF variable ID for the height coordinate
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_eutm_soil      !< NetCDF variable ID for E_UTM for soil quantity
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_nutm_soil      !< NetCDF variable ID for N_UTM for soil quantity
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_depth          !< NetCDF variable ID for soil depth 
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  id_var_time           !< NetCDF variable ID for the time coordinate
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  start_count_time      !< NetCDF start index for the time dimension
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  start_count_utm       !< NetCDF start index for the UTM dimension 
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  start_count_utm_soil  !< NetCDF start index for the UTM dimension in the soil
+    
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  id_var              !< NetCDF variable IDs for the sampled variables at a site
+    
+    LOGICAL, DIMENSION(:), ALLOCATABLE ::  soil  !< flag indicating sampled soil quantities
+    
+    REAL(wp)                              ::  output_time    !< output time   
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var            !< sampled data
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_soil       !< sampled data of a soil varialbe
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  origin_x_obs   !< site coordinate (x)
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  origin_y_obs   !< site coordinate (y)
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  e_utm          !< E_UTM coordinates where measurements were taken
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  n_utm          !< N_UTM coordinates where measurements were taken
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  z_ag           !< height coordinates where measurements were taken
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  e_utm_soil     !< E_UTM coordinates where measurements were taken (soil)
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  n_utm_soil     !< N_UTM coordinates where measurements were taken (soil)
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  depth          !< soil depth where measurements were taken (soil)
+       
+!
+!-- Read namelist.
+    CALL cvm_parin
+!
+!-- Create filename suffix of the treated binary file.
+    f = 0
+    CALL create_file_string
+!
+!-- Open binary file for processor 0.
+    OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                   &
+           TRIM( myid_char ), FORM = 'UNFORMATTED' )
+!
+!-- Reader global information, such as number of stations, their type, 
+!-- number of observation coordinates for each station on subdomain and 
+!-- total. Read number of sites. 
+    READ( file_id_in )  char_in   
+    READ( file_id_in )  nvm
+!
+!-- Allocate arrays required to describe measurements
+    ALLOCATE( site(1:nvm)          )
+    ALLOCATE( filename(1:nvm)      )
+    ALLOCATE( feature_type(1:nvm)  )
+    ALLOCATE( soil_quantity(1:nvm) )
+    ALLOCATE( nvar(1:nvm)          )
+    ALLOCATE( origin_x_obs(1:nvm)  )
+    ALLOCATE( origin_y_obs(1:nvm)  )    
+    ALLOCATE( ns(1:nvm)            )
+    ALLOCATE( ns_tot(1:nvm)        )
+    ALLOCATE( ns_soil(1:nvm)       )
+    ALLOCATE( ns_soil_tot(1:nvm)   )
+    ALLOCATE( soil(1:nvm)          )
+    
+    ns_soil = 0
+    ns_soil_tot = 0  
+!
+!-- Allocate array with the measured variables at each station
+    ALLOCATE( variables(1:100,1:nvm) )
+!
+!-- Allocate arrays for NetCDF IDs
+    ALLOCATE( nc_id(1:nvm)        )
+    ALLOCATE( id_eutm(1:nvm)      )
+    ALLOCATE( id_nutm(1:nvm)      )
+    ALLOCATE( id_hao(1:nvm)       )
+    ALLOCATE( id_eutm_soil(1:nvm) )
+    ALLOCATE( id_nutm_soil(1:nvm) )
+    ALLOCATE( id_depth(1:nvm)     )
+    ALLOCATE( id_time(1:nvm)      )
+    ALLOCATE( id_var_eutm(1:nvm)  )
+    ALLOCATE( id_var_nutm(1:nvm)  )
+    ALLOCATE( id_var_hao(1:nvm)   ) 
+    ALLOCATE( id_var_eutm_soil(1:nvm) )
+    ALLOCATE( id_var_nutm_soil(1:nvm) )
+    ALLOCATE( id_var_depth(1:nvm) )    
+    ALLOCATE( id_var_time(1:nvm)  )
+    ALLOCATE( id_var(1:50,1:nvm)  )
+    id_var = 0    
+    nc_id  = 0
+!
+!-- Allocate arrays that contain information about the start index in the 
+!-- dimension array, used to write binary data at the correct position in
+!-- the NetCDF file.    
+    ALLOCATE( start_count_utm(1:nvm) )
+    ALLOCATE( start_count_utm_soil(1:nvm) )
+    
+    ALLOCATE( start_count_time(1:nvm) )
+!
+!-- Read further global information from processor 0, such as filenames,
+!-- global attributes, dimension sizes, etc. used to create NetCDF files.
+    DO  l = 1, nvm
+!
+!--    Read sitename
+       READ( file_id_in )  char_in
+       READ( file_id_in )  site(l)   
+!
+!--    Read filename (original name where real-world data is stored) 
+       READ( file_id_in )  char_in
+       READ( file_id_in )  filename(l)
+!
+!--    Read type of the measurement
+       READ( file_id_in )  char_in
+       READ( file_id_in )  feature_type(l)
+!
+!--    Read x-y origin coordinates (in UTM)
+       READ( file_id_in )  char_in
+       READ( file_id_in )  origin_x_obs(l)
+       READ( file_id_in )  char_in
+       READ( file_id_in )  origin_y_obs(l)
+!
+!--    Read total number of observation grid points (dimension size of the 
+!--    virtual measurement)
+       READ( file_id_in )  char_in
+       READ( file_id_in )  ns_tot(l)
+!
+!--    Read number of observed quantities at each station
+       READ( file_id_in )  char_in
+       READ( file_id_in )  nvar(l)
+!
+!--    Read names of observed quantities 
+       READ( file_id_in )  char_in
+       READ( file_id_in )  variables(1:nvar(l),l)
+!
+!--    Further dummy arguments are read (number of observation points 
+!--    on subdomains and its UTM coordinates). 
+       READ( file_id_in )  char_in
+       READ( file_id_in )  ns(l)
+             
+       ALLOCATE( e_utm(1:ns(l)) )
+       ALLOCATE( n_utm(1:ns(l)) )
+       ALLOCATE( z_ag(1:ns(l)) )
+!
+!--    Read the local coordinate arrays
+       READ( file_id_in )  char_in
+       READ( file_id_in )  e_utm 
+      
+       READ( file_id_in )  char_in
+       READ( file_id_in )  n_utm
+
+       READ( file_id_in )  char_in
+       READ( file_id_in )  z_ag
+       
+       DEALLOCATE( e_utm )
+       DEALLOCATE( n_utm )
+       DEALLOCATE( z_ag  )  
+       
+!
+!--    Read flag indicating whether soil data is also present or not
+       READ( file_id_in )  char_in
+       READ( file_id_in )  char_in
+       
+       soil(l) = MERGE( .TRUE., .FALSE., TRIM( char_in ) == "yes" )
+       
+       IF ( soil(l) )  THEN
+       
+          READ( file_id_in )  char_in
+          READ( file_id_in )  ns_soil_tot(l)
+          
+          READ( file_id_in )  char_in
+          READ( file_id_in )  ns_soil(l)
+       
+          ALLOCATE( e_utm_soil(1:ns_soil(l)) )
+          ALLOCATE( n_utm_soil(1:ns_soil(l)) )
+          ALLOCATE( depth(1:ns_soil(l)) )
+!
+!--       Read the local coordinate arrays
+          READ( file_id_in )  char_in
+          READ( file_id_in )  e_utm_soil 
+      
+          READ( file_id_in )  char_in
+          READ( file_id_in )  n_utm_soil
+
+          READ( file_id_in )  char_in
+          READ( file_id_in )  depth
+       
+          DEALLOCATE( e_utm_soil )
+          DEALLOCATE( n_utm_soil )
+          DEALLOCATE( depth      )
+       
+       ENDIF
+!
+!--    Create netcdf file and setup header information
+       CALL netcdf_create_file
+       
+    ENDDO
+!
+!-- Close binary file created by processor 0.
+    CLOSE( file_id_in )
+!
+!-- Initialize UTM coordinate start index
+    start_count_utm = 1
+    start_count_utm_soil = 1
+!
+!-- Read data from all PEs and write into file
+    DO  f = 0, num_pe - 1
+!
+!--    Create filename suffix of the treated binary file.
+       CALL create_file_string
+!
+!--    Open binary file for processor f.
+       OPEN ( file_id_in, FILE = TRIM( path ) // TRIM( run ) //                &
+              TRIM( myid_char ), FORM = 'UNFORMATTED' )
+!
+!--    Initialize time coordinate start index
+       start_count_time = 1
+!
+!--    Reader global information, such as number of stations, their type, 
+!--    number of observation coordinates for each station on subdomain and 
+!--    total. 
+!--    Again, read number of sites. 
+       READ( file_id_in ) char_in
+       READ( file_id_in ) nvm
+             
+       DO  l = 1, nvm
+!
+!--       Read sitename
+          READ( file_id_in )  char_in
+          READ( file_id_in )  site(l)   
+!
+!--       Read filename (original name where real-world data is stored) 
+          READ( file_id_in )  char_in
+          READ( file_id_in )  filename(l)
+!
+!--       Read type of the measurement
+          READ( file_id_in )  char_in
+          READ( file_id_in )  feature_type(l)
+!
+!--       Read x-y origin coordinates (in UTM)
+          READ( file_id_in )  char_in
+          READ( file_id_in )  origin_x_obs(l)
+          READ( file_id_in )  char_in
+          READ( file_id_in )  origin_y_obs(l)
+!
+!--       Read total number of observation grid points (dimension size of the 
+!--       virtual measurement)
+          READ( file_id_in )  char_in
+          READ( file_id_in )  ns_tot(l)
+!
+!--       Read number of observed quantities at each station
+          READ( file_id_in )  char_in
+          READ( file_id_in )  nvar(l)
+!
+!--       Read names of observed quantities 
+          READ( file_id_in )  char_in
+          READ( file_id_in )  variables(1:nvar(l),l)
+!
+!--       Further dummy arguments are read (number of observation points 
+!--       on subdomains and its UTM coordinates). 
+          READ( file_id_in )  char_in
+          READ( file_id_in )  ns(l)
+             
+          ALLOCATE( e_utm(1:ns(l)) )
+          ALLOCATE( n_utm(1:ns(l)) )
+          ALLOCATE( z_ag(1:ns(l)) )
+!
+!--       Read the local coordinate arrays
+          READ( file_id_in )  char_in
+          READ( file_id_in )  e_utm 
+      
+          READ( file_id_in )  char_in
+          READ( file_id_in )  n_utm
+   
+          READ( file_id_in )  char_in
+          READ( file_id_in )  z_ag      
+!
+!--       Read flag indicating whether soil data is also present or not
+          READ( file_id_in )  char_in
+          READ( file_id_in )  char_in
+       
+          soil(l) = MERGE( .TRUE., .FALSE., TRIM( char_in ) == "yes" )
+       
+          IF ( soil(l) )  THEN
+       
+             READ( file_id_in )  char_in
+             READ( file_id_in )  ns_soil_tot(l)
+          
+             READ( file_id_in )  char_in
+             READ( file_id_in )  ns_soil(l)
+       
+             ALLOCATE( e_utm_soil(1:ns_soil(l)) )
+             ALLOCATE( n_utm_soil(1:ns_soil(l)) )
+             ALLOCATE( depth(1:ns_soil(l)) )
+!
+!--          Read the local coordinate arrays
+             READ( file_id_in )  char_in
+             READ( file_id_in )  e_utm_soil 
+      
+             READ( file_id_in )  char_in
+             READ( file_id_in )  n_utm_soil
+   
+             READ( file_id_in )  char_in
+             READ( file_id_in )  depth
+       
+          ENDIF
+!
+!--       Setup the NetCDF dimensions for the UTM coordinates
+          CALL netcdf_define_utm_dimension
+          
+          DEALLOCATE( e_utm )
+          DEALLOCATE( n_utm )
+          DEALLOCATE( z_ag )
+          
+          IF ( soil(l) )  THEN
+             DEALLOCATE( e_utm_soil )
+             DEALLOCATE( n_utm_soil )
+             DEALLOCATE( depth      )
+          ENDIF
+          
+       ENDDO
+!
+!--    Read the actual data, starting with the identification string for the 
+!--    output time
+       READ( file_id_in )  char_in          
+       DO WHILE ( TRIM( char_in ) == 'output time')
+       
+          READ( file_id_in )  output_time 
+!
+!--       Loop over all sites
+          DO  l = 1, nvm
+!
+!--          Cycle loop if no observation coordinates are on local subdomain
+             IF ( ns(l) < 1  .AND.  ns_soil(l) < 1 )  CYCLE
+!
+!--          Setup time coordinate
+             CALL netcdf_define_time_dimension
+!
+!--          Read the actual data, therefore, allocate appropriate array with
+!--          size of the subdomain coordinates. Output data immediately into 
+!--          NetCDF file. 
+             ALLOCATE( var(1:ns(l)) )
+             IF ( soil(l) )  ALLOCATE( var_soil(1:ns_soil(l)) )
+             
+             DO  n = 1, nvar(l)
+                READ( file_id_in )  variables(n,l)
+                IF ( soil(l)  .AND.                                            &
+                     ANY( TRIM( variables(n,l) ) == soil_quantities ) )  THEN
+                   READ( file_id_in )  var_soil
+                ELSE
+                   READ( file_id_in )  var
+                ENDIF
+!
+!--             Write data to NetCDF file 
+                CALL netcdf_data_output
+             ENDDO        
+             
+             DEALLOCATE( var )
+             IF( ALLOCATED(var_soil) )  DEALLOCATE( var_soil )
+!
+!--          Increment NetCDF index of the time coordinate
+             start_count_time(l) = start_count_time(l) + 1
+          ENDDO
+!
+!--       Read next identification string
+          READ( file_id_in )  char_in
+          
+       ENDDO
+!
+!--    After data from processor f is read and output into NetCDF file, 
+!--    the start index of the UTM coordinate array need to be incremented
+       start_count_utm      = start_count_utm      + ns
+       start_count_utm_soil = start_count_utm_soil + ns_soil
+!
+!--    Close binary file for processor f
+       CLOSE( file_id_in )
+    ENDDO
+!
+!-- Close Netcdf files
+    DO  l = 1, nvm
+       CALL netcdf_close_file
+    ENDDO
+    
+ CONTAINS
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine read the namelist file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE cvm_parin
+       
+       IMPLICIT NONE
+       
+       INTEGER(iwp) ::  file_id_parin = 90
+       
+       NAMELIST /vm/  cycle_number, num_pe, path, run
+
+!
+!--    Open namelist file.
+       OPEN( file_id_parin, FILE='vm_parin', STATUS='OLD', FORM='FORMATTED')
+!
+!--    Read namelist.
+       READ ( file_id_parin, vm )
+!
+!--    Close namelist file.
+       CLOSE( file_id_parin )
+       
+    END SUBROUTINE cvm_parin
+      
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine creates the filename string of the treated binary file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE create_file_string
+       
+       IMPLICIT NONE
+       
+       CHARACTER(LEN=4)   ::  char_cycle = '' !< dummy string for cycle number
+       CHARACTER(LEN=10)  ::  char_dum        !< dummy string for processor ID
+       
+!
+!--    Create substring for the cycle number.
+       IF ( cycle_number /= 0 )  THEN
+          IF ( cycle_number < 10 )  THEN
+             WRITE( char_cycle, '(I1)')  cycle_number
+             char_cycle = '.00' // TRIM( char_cycle )
+          ELSEIF ( cycle_number < 100 )  THEN
+             WRITE( char_cycle, '(I2)')  cycle_number
+             char_cycle = '.0' // TRIM( char_cycle )
+          ELSEIF ( cycle_number < 1000 )  THEN
+             WRITE( char_cycle, '(I3)')  cycle_number
+             char_cycle = '.' // TRIM( char_cycle )
+          ENDIF
+       ENDIF
+!
+!--    Create substring for the processor id and combine all substrings.
+       IF ( f < 10 )  THEN
+          WRITE( char_dum, '(I1)')  f
+          myid_char = '_vmeas_00000' // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 100     )  THEN
+          WRITE( char_dum, '(I2)')  f
+          myid_char = '_vmeas_0000'  // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 1000    )  THEN
+          WRITE( char_dum, '(I3)')  f
+          myid_char = '_vmeas_000'   // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 10000   )  THEN
+          WRITE( char_dum, '(I4)')  f
+          myid_char = '_vmeas_00'    // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 100000  )  THEN
+          WRITE( char_dum, '(I5)')  f
+          myid_char = '_vmeas_0'     // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ELSEIF ( f < 1000000 )  THEN
+          WRITE( char_dum, '(I6)')  f
+          myid_char = '_vmeas_'      // TRIM( char_dum ) //                    &
+                      TRIM( char_cycle ) // file_suffix
+       ENDIF
+       
+    END SUBROUTINE create_file_string
+    
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine creates the NetCDF file and defines dimesions and varialbes.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_create_file
+       
+       IMPLICIT NONE
+       
+       
+       CHARACTER(LEN=5)   ::  char_cycle = ''  !< dummy string for cycle number
+       CHARACTER(LEN=200) ::  nc_filename = '' !< NetCDF filename
+       
+!
+!--    Create substring for the cycle number.
+       IF ( cycle_number /= 0 )  THEN
+          IF ( cycle_number < 10 )  THEN
+             WRITE( char_cycle, '(I1)')  cycle_number
+             char_cycle = '.00' // TRIM( char_cycle ) // '.'
+          ELSEIF ( cycle_number < 100 )  THEN
+             WRITE( char_cycle, '(I2)')  cycle_number
+             char_cycle = '.0' // TRIM( char_cycle ) // '.'
+          ELSEIF ( cycle_number < 1000 )  THEN
+             WRITE( char_cycle, '(I3)')  cycle_number
+             char_cycle = '.' // TRIM( char_cycle ) // '.'
+          ENDIF
+       ELSE
+          char_cycle = '.'
+       ENDIF
+#if defined( __netcdf )     
+
+       nc_filename = site(l)(1:LEN_TRIM(site(l))-1) // '_palm4U' //            &
+                     TRIM( char_cycle )  // 'nc'
+!
+!--    Create NetCDF file
+       status_nc = NF90_CREATE( nc_filename(1:LEN_TRIM(nc_filename)),          &
+                                IOR( NF90_CLOBBER, NF90_NETCDF4 ), nc_id(l) )
+       CALL handle_error( "create file" )
+!
+!--    Define attributes
+       status_nc = NF90_PUT_ATT( nc_id(l), NF90_GLOBAL, "featureType",         &
+                                 TRIM( feature_type(l) ) )
+       CALL handle_error( "define attribue featureType" )
+       
+       status_nc = NF90_PUT_ATT( nc_id(l), NF90_GLOBAL, "origin_x",            &
+                                 origin_x_obs(l) )
+       CALL handle_error( "define attribue origin_x" )
+       
+       status_nc = NF90_PUT_ATT( nc_id(l), NF90_GLOBAL, "origin_y",            &
+                                 origin_y_obs(l) )
+       CALL handle_error( "define attribue origin_y" )
+       
+       status_nc = NF90_PUT_ATT( nc_id(l), NF90_GLOBAL, "site",                &
+                                 TRIM( site(l) ) )
+       CALL handle_error( "define attribue site" )
+!
+!--    Define dimensions
+       status_nc = NF90_DEF_DIM( nc_id(l), 'E_UTM', ns_tot(l), id_eutm(l) )
+       CALL handle_error( "define dimension E_UTM" )
+       
+       status_nc = NF90_DEF_DIM( nc_id(l), 'N_UTM', ns_tot(l), id_nutm(l) )
+       CALL handle_error( "define dimension N_UTM" )
+       
+       status_nc = NF90_DEF_DIM( nc_id(l), 'height_above_origin', ns_tot(l),   &
+                                 id_hao(l) )
+       CALL handle_error( "define dimension height_above_origin" )
+       
+       IF ( soil(l) )  THEN
+          status_nc = NF90_DEF_DIM( nc_id(l), 'E_UTM soil', ns_soil_tot(l),    &
+                                    id_eutm_soil(l) )
+          CALL handle_error( "define dimension E_UTM soil" )
+          
+          status_nc = NF90_DEF_DIM( nc_id(l), 'N_UTM soil', ns_soil_tot(l),    &
+                                    id_nutm_soil(l) )
+          CALL handle_error( "define dimension N_UTM soil" )
+          status_nc = NF90_DEF_DIM( nc_id(l), 'depth', ns_soil_tot(l),         &
+                                    id_depth(l) )
+          CALL handle_error( "define dimension depth" )
+       ENDIF
+       
+       status_nc = NF90_DEF_DIM( nc_id(l), 'time', NF90_UNLIMITED, id_time(l) )
+       CALL handle_error( "define dimension time" )      
+
+!
+!--    Define dimension variables
+       status_nc = NF90_DEF_VAR( nc_id(l), 'E_UTM', NF90_DOUBLE,               &
+                                 (/ id_eutm(l) /), id_var_eutm(l) )
+       CALL handle_error( "define variable E_UTM" )
+       
+       status_nc = NF90_DEF_VAR( nc_id(l), 'N_UTM', NF90_DOUBLE,               &
+                                 (/ id_nutm(l) /), id_var_nutm(l) )
+       CALL handle_error( "define variable N_UTM" )
+       status_nc = NF90_DEF_VAR( nc_id(l), 'height_above_origin', NF90_DOUBLE, &
+                                 (/ id_hao(l)  /), id_var_hao(l)  )
+       CALL handle_error( "define variable height_above_origin" )
+       status_nc = NF90_DEF_VAR( nc_id(l), 'time', NF90_DOUBLE,                &
+                                 (/ id_time(l) /), id_var_time(l) )
+       CALL handle_error( "define variable time" )
+       
+       IF ( soil(l) )  THEN
+          status_nc = NF90_DEF_VAR( nc_id(l), 'E_UTM soil', NF90_DOUBLE,       &
+                                    (/ id_eutm_soil(l) /), id_var_eutm_soil(l) )
+          CALL handle_error( "define variable E_UTM soil" )
+          
+          status_nc = NF90_DEF_VAR( nc_id(l), 'N_UTM soil', NF90_DOUBLE,       &
+                                    (/ id_nutm_soil(l) /), id_var_nutm_soil(l) )
+          CALL handle_error( "define variable N_UTM soil" )
+          
+          status_nc = NF90_DEF_VAR( nc_id(l), 'depth', NF90_DOUBLE,            &
+                                    (/ id_depth(l)  /), id_var_depth(l) )
+          CALL handle_error( "define variable depth" )
+       ENDIF
+!
+!--    Define the measured quantities
+       DO  n = 1, nvar(l)
+          IF ( soil(l)  .AND.                                                  &
+               ANY( TRIM( variables(n,l) ) == soil_quantities ) )  THEN
+             status_nc = NF90_DEF_VAR( nc_id(l), TRIM( variables(n,l) ),       &
+                                       NF90_DOUBLE,                            &
+                                       (/ id_time(l), id_eutm_soil(l) /),      &
+                                       id_var(n,l) )
+             CALL handle_error( "define variable " // TRIM( variables(n,l) ) )
+          ELSE
+             status_nc = NF90_DEF_VAR( nc_id(l), TRIM( variables(n,l) ),       &
+                                       NF90_DOUBLE,                            &
+                                       (/ id_time(l), id_eutm(l) /),           &
+                                       id_var(n,l) )
+             CALL handle_error( "define variable " // TRIM( variables(n,l) ) )
+          ENDIF
+       ENDDO
+#endif        
+    END SUBROUTINE netcdf_create_file
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine closes a NetCDF file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_close_file
+       
+       IMPLICIT NONE
+       
+#if defined( __netcdf )        
+       status_nc = NF90_CLOSE( nc_id(l) )
+       CALL handle_error( "close file" )
+#endif   
+       
+    END SUBROUTINE netcdf_close_file
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine defines the UTM dimensions
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_define_utm_dimension
+       
+       IMPLICIT NONE
+       
+!
+!--    Define dimensions
+#if defined( __netcdf ) 
+       status_nc = NF90_PUT_VAR( nc_id(l), id_var_eutm(l), e_utm,              &
+                                 start = (/ start_count_utm(l) /),             & 
+                                 count = (/ ns(l) /) )
+       CALL handle_error( "write variable E_UTM" )
+       
+       status_nc = NF90_PUT_VAR( nc_id(l), id_var_nutm(l), n_utm,              &
+                                 start = (/ start_count_utm(l) /),             &
+                                 count = (/ ns(l) /) )
+       CALL handle_error( "write variable N_UTM" )
+       
+       status_nc = NF90_PUT_VAR( nc_id(l), id_var_hao(l),  z_ag,               &
+                                 start = (/ start_count_utm(l) /),             &
+                                 count = (/ ns(l) /) )
+       CALL handle_error( "write variable height_above_origin" )
+       
+       IF ( soil(l) )  THEN
+          status_nc = NF90_PUT_VAR( nc_id(l), id_var_eutm_soil(l), e_utm_soil, &
+                                 start = (/ start_count_utm_soil(l) /),        & 
+                                 count = (/ ns_soil(l) /) )
+          CALL handle_error( "write variable E_UTM soil" )
+          
+          status_nc = NF90_PUT_VAR( nc_id(l), id_var_nutm_soil(l), n_utm_soil, &
+                                 start = (/ start_count_utm_soil(l) /),        &
+                                 count = (/ ns_soil(l) /) )
+          CALL handle_error( "write variable N_UTM soil" )
+          
+          status_nc = NF90_PUT_VAR( nc_id(l), id_var_depth(l),  depth,         &
+                                 start = (/ start_count_utm_soil(l) /),        &
+                                 count = (/ ns_soil(l) /) )
+          CALL handle_error( "write variable depth" )
+       ENDIF
+!
+!--    End of NetCDF file definition
+       status_nc = NF90_ENDDEF( nc_id(l) )
+#endif       
+    END SUBROUTINE netcdf_define_utm_dimension
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine updates the unlimited time dimension.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_define_time_dimension
+       
+       IMPLICIT NONE
+       
+!
+!--    Define dimensions
+#if defined( __netcdf )   
+       status_nc = NF90_PUT_VAR( nc_id(l), id_var_time(l), (/ output_time /),  &
+                                 start = (/ start_count_time(l) /),            &
+                                 count = (/ 1 /) )
+       CALL handle_error( "write variable time" )
+#endif
+       
+    END SUBROUTINE netcdf_define_time_dimension
+
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine writes the sampled variable to the NetCDF file.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_output
+       
+       IMPLICIT NONE
+       
+
+       IF ( soil(l)  .AND.                                                     &
+            ANY( TRIM( variables(n,l) ) == soil_quantities ) )  THEN
+          status_nc = NF90_PUT_VAR( nc_id(l), id_var(n,l), (/ var_soil /),     &
+                                    start = (/ start_count_time(l),            &
+                                               start_count_utm_soil(l) /),     &
+                                    count = (/ 1, ns_soil(l) /) )
+          CALL handle_error( "write variable " // TRIM( variables(n,l) ) )
+       ELSE
+          status_nc = NF90_PUT_VAR( nc_id(l), id_var(n,l), (/ var /),          &
+                                    start = (/ start_count_time(l),            &
+                                               start_count_utm(l) /),          &
+                                    count = (/ 1, ns(l) /) )
+          CALL handle_error( "write variable " // TRIM( variables(n,l) ) )
+       ENDIF
+       
+    END SUBROUTINE netcdf_data_output
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> NetCDF error handling.
+!------------------------------------------------------------------------------!    
+    SUBROUTINE handle_error( action )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  action !< string indicating the current file action
+       
+#if defined( __netcdf )       
+       IF ( status_nc /= NF90_NOERR )  THEN
+          PRINT*, TRIM( NF90_STRERROR( status_nc ) ) // ' -- ' // action 
+          STOP
+       ENDIF
+#endif
+
+    END SUBROUTINE handle_error
+        
+ END PROGRAM combine_virtual_measurements
Index: palm/trunk/UTIL/combine_virtual_measurements/vm_parin
===================================================================
--- palm/trunk/UTIL/combine_virtual_measurements/vm_parin	(revision 3704)
+++ palm/trunk/UTIL/combine_virtual_measurements/vm_parin	(revision 3704)
@@ -0,0 +1,14 @@
+&vm
+!
+!-- Number of processors of the run.
+    num_pe = 36, 
+!
+!-- Path to the output binary files.
+    path = '~/palm/current_version/JOBS/test_urban/OUTPUT/', 
+!
+!-- Name of the run.
+    run = 'test_berlin',
+!
+!-- Cycle number of the output data which should be processed
+    cycle_number = 0, 
+/  
