Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3469)
+++ palm/trunk/SOURCE/Makefile	(revision 3472)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# Add virtual measurement module
 # 
 # Former revisions:
@@ -696,4 +696,5 @@
 	vertical_nesting_mod.f90 \
 	virtual_flight_mod.f90 \
+	virtual_measurement_mod.f90 \
 	wind_turbine_model_mod.f90 \
 	wrd_write_string.f90 \
@@ -1134,4 +1135,5 @@
 	uv_exposure_model_mod.o \
 	virtual_flight_mod.o \
+	virtual_measurement_mod.o \
 	wind_turbine_model_mod.o
 init_advec.o: \
@@ -1451,4 +1453,5 @@
 	vertical_nesting_mod.o \
 	virtual_flight_mod.o \
+	virtual_measurement_mod.o \
 	wind_turbine_model_mod.o
 plant_canopy_model_mod.o: \
@@ -1750,4 +1753,5 @@
 	vertical_nesting_mod.o \
 	virtual_flight_mod.o \
+	virtual_measurement_mod.o \
 	wind_turbine_model_mod.o
 time_integration_spinup.o: \
@@ -1968,4 +1972,10 @@
 	user_flight.o \
 	user_init_flight.o
+virtual_measurement_mod.o: \
+	cpulog_mod.o \
+	mod_kinds.o \
+	modules.o \
+	netcdf_data_input_mod.o \
+	netcdf_interface_mod.o
 wind_turbine_model_mod.o: \
 	basic_constants_and_equations_mod.o \
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 3469)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 3472)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Add virtual measurement module
 ! 
 ! Former revisions:
@@ -664,4 +664,7 @@
     USE uv_exposure_model_mod,                                                 &
         ONLY:  uvem_init, uvem_init_arrays
+        
+    USE virtual_measurement_mod,                                               &
+        ONLY:  vm_init
 
     USE wind_turbine_model_mod,                                                &
@@ -2443,4 +2446,8 @@
        CALL surface_output_init
     ENDIF
+!
+!-- If virtual measurements should be taken, initialize all relevant 
+!-- arrays and quantities.
+    IF ( virtual_measurement )  CALL vm_init
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3469)
+++ palm/trunk/SOURCE/modules.f90	(revision 3472)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +virtual_measurement
 ! 
 ! Former revisions:
@@ -1398,5 +1398,6 @@
     LOGICAL ::  use_upstream_for_tke = .FALSE.                   !< namelist parameter
     LOGICAL ::  uv_exposure = .FALSE.                            !< switch for uv exposure model
-    LOGICAL ::  virtual_flight = .FALSE.                         !< use virtual flight model?
+    LOGICAL ::  virtual_flight = .FALSE.                         !< use virtual flight model
+    LOGICAL ::  virtual_measurement = .FALSE.                    !< control parameter to switch-on virtual measurements
     LOGICAL ::  wall_adjustment = .TRUE.                         !< namelist parameter
     LOGICAL ::  wind_turbine = .FALSE.                           !< flag for use of wind turbine model
Index: palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3469)
+++ palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3472)
@@ -228,5 +228,5 @@
 !> Modulue contains routines to input data according to Palm input data
 !> standart using dynamic and static input files.
-!> @todo - Review Reading of netcdf files for chemistry
+!> @todo - Chemistry: revise reading of netcdf file and ajdust formatting according to standard!!!
 !> @todo - Order input alphabetically
 !> @todo - Revise error messages and error numbers
@@ -670,8 +670,9 @@
     CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC' !< Name of file which comprises dynamic input data
     CHARACTER(LEN=100) ::  input_file_chem    = 'PIDS_CHEM'    !< Name of file which comprises chemistry input data
-
-    CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:)    :: string_values  !< output of string variables read from netcdf input files
+    CHARACTER(LEN=100) ::  input_file_vm      = 'PIDS_VM'      !< Name of file which comprises virtual measurement data
+
+    CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:)    ::  string_values  !< output of string variables read from netcdf input files
  
-    INTEGER(iwp)                                     :: id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
+    INTEGER(iwp)                                     ::  id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
 
     INTEGER(iwp) ::  nc_stat         !< return value of nf90 function call
@@ -679,6 +680,7 @@
     LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
     LOGICAL ::  input_pids_dynamic = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing dynamic information exists
-    LOGICAL ::  input_pids_chem = .FALSE.      !< Flag indicating whether Palm-input-data-standard file containing chemistry information exists
-
+    LOGICAL ::  input_pids_chem    = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing chemistry information exists
+    LOGICAL ::  input_pids_vm      = .FALSE.   !< Flag indicating whether input file for virtual measurements exist
+    
     LOGICAL ::  collective_read = .FALSE.      !< Enable NetCDF collective read
 
@@ -707,4 +709,8 @@
        MODULE PROCEDURE netcdf_data_input_chemistry_data
     END INTERFACE netcdf_data_input_chemistry_data
+    
+    INTERFACE netcdf_data_input_get_dimension_length                        
+       MODULE PROCEDURE netcdf_data_input_get_dimension_length
+    END INTERFACE netcdf_data_input_get_dimension_length
 
     INTERFACE netcdf_data_input_inquire_file
@@ -715,4 +721,10 @@
        MODULE PROCEDURE netcdf_data_input_init
     END INTERFACE netcdf_data_input_init
+    
+    INTERFACE netcdf_data_input_att
+       MODULE PROCEDURE netcdf_data_input_att_int
+       MODULE PROCEDURE netcdf_data_input_att_real
+       MODULE PROCEDURE netcdf_data_input_att_string
+    END INTERFACE netcdf_data_input_att
 
     INTERFACE netcdf_data_input_init_3d
@@ -732,10 +744,12 @@
     END INTERFACE netcdf_data_input_surface_data
 
-    INTERFACE netcdf_data_input_topo
-       MODULE PROCEDURE netcdf_data_input_topo
-    END INTERFACE netcdf_data_input_topo
+    INTERFACE netcdf_data_input_var
+       MODULE PROCEDURE netcdf_data_input_var_char
+       MODULE PROCEDURE netcdf_data_input_var_real_1d
+       MODULE PROCEDURE netcdf_data_input_var_real_2d
+    END INTERFACE netcdf_data_input_var
 
     INTERFACE get_variable
-       MODULE PROCEDURE get_variable_string
+       MODULE PROCEDURE get_variable_1d_char
        MODULE PROCEDURE get_variable_1d_int
        MODULE PROCEDURE get_variable_1d_real
@@ -749,4 +763,5 @@
        MODULE PROCEDURE get_variable_4d_real
        MODULE PROCEDURE get_variable_5d_to_4d_real
+       MODULE PROCEDURE get_variable_string       
     END INTERFACE get_variable
 
@@ -769,5 +784,6 @@
            coord_ref_sys,                                                      &
            init_3d, init_model, input_file_static, input_pids_static,          &
-           input_pids_dynamic, leaf_area_density_f, nest_offl,                 &
+           input_pids_dynamic, input_pids_vm, input_file_vm,                   &
+           leaf_area_density_f, nest_offl,                                     &
            pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
            root_area_density_lad_f, root_area_density_lsm_f, soil_pars_f,      &
@@ -779,10 +795,12 @@
 !-- Public subroutines
     PUBLIC netcdf_data_input_check_dynamic, netcdf_data_input_check_static,    &
-           netcdf_data_input_chemistry_data, netcdf_data_input_inquire_file,   &
+           netcdf_data_input_chemistry_data,                                   &
+           netcdf_data_input_get_dimension_length,                             &
+           netcdf_data_input_inquire_file,                                     &
            netcdf_data_input_init, netcdf_data_input_init_lsm,                 &
-           netcdf_data_input_init_3d,                                          &
+           netcdf_data_input_init_3d, netcdf_data_input_att,                   &
            netcdf_data_input_interpolate, netcdf_data_input_offline_nesting,   &
            netcdf_data_input_surface_data, netcdf_data_input_topo,             &
-           get_attribute, get_dimension_length, get_variable, open_read_file
+           netcdf_data_input_var, get_attribute, get_variable, open_read_file
 
  CONTAINS
@@ -802,11 +820,12 @@
 
 #if defined ( __netcdf )
-       INQUIRE( FILE = TRIM( input_file_static ) // TRIM( coupling_char ),     &
+       INQUIRE( FILE = TRIM( input_file_static )   // TRIM( coupling_char ),   &
                 EXIST = input_pids_static  )
        INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ),    &
                 EXIST = input_pids_dynamic )
-       INQUIRE( FILE = TRIM( input_file_chem ) // TRIM( coupling_char ),    &
+       INQUIRE( FILE = TRIM( input_file_chem )    // TRIM( coupling_char ),    &
                 EXIST = input_pids_chem )
-
+       INQUIRE( FILE = TRIM( input_file_vm )      // TRIM( coupling_char ),    &
+                EXIST = input_pids_vm )
 #endif
 
@@ -933,4 +952,213 @@
 
     END SUBROUTINE netcdf_data_input_init
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read an array of characters.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_var_char( val, search_string, id_mod )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string     !< name of the variable
+       CHARACTER(LEN=*), DIMENSION(:) ::  val !< variable which should be read
+       
+       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
+
+#if defined ( __netcdf )
+!
+!--    Read variable
+       CALL get_variable( id_mod, search_string, val )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_var_char
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read an 1D array of REAL values.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_var_real_1d( val, search_string, id_mod )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string     !< name of the variable      
+       
+       INTEGER(iwp) ::  id_mod        !< NetCDF id of input file
+       
+       REAL(wp), DIMENSION(:) ::  val !< variable which should be read
+
+#if defined ( __netcdf )
+!
+!--    Read variable
+       CALL get_variable( id_mod, search_string, val )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_var_real_1d
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read an 1D array of REAL values.
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_var_real_2d( val, search_string,              &
+                                              id_mod, d1s, d1e, d2s, d2e )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string     !< name of the variable     
+       
+       INTEGER(iwp) ::  id_mod  !< NetCDF id of input file
+       INTEGER(iwp) ::  d1e     !< end index of first dimension to be read
+       INTEGER(iwp) ::  d2e     !< end index of second dimension to be read
+       INTEGER(iwp) ::  d1s     !< start index of first dimension to be read
+       INTEGER(iwp) ::  d2s     !< start index of second dimension to be read
+       
+       REAL(wp), DIMENSION(:,:) ::  val !< variable which should be read
+
+#if defined ( __netcdf )
+!
+!--    Read character variable
+       CALL get_variable( id_mod, search_string, val, d1s, d1e, d2s, d2e )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_var_real_2d
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read a global string attribute 
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_att_string( val, search_string, id_mod,       &
+                                             input_file, global, openclose,    &
+                                             variable_name )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string !< name of the attribue
+       CHARACTER(LEN=*) ::  val           !< attribute
+       
+       CHARACTER(LEN=*) ::  input_file    !< name of input file
+       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
+       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed  
+       
+       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
+       
+       LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
+
+#if defined ( __netcdf )
+!
+!--    Open file in read-only mode if necessary
+       IF ( openclose == 'open' )  THEN
+          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
+                                  id_mod )
+       ENDIF
+!
+!--    Read global attribute
+       IF ( global )  THEN
+          CALL get_attribute( id_mod, search_string, val, global )
+!
+!--    Read variable attribute
+       ELSE
+          CALL get_attribute( id_mod, search_string, val, global, variable_name )
+       ENDIF
+!
+!--    Close input file 
+       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_att_string
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read a global integer attribute 
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_att_int( val, search_string, id_mod,          &
+                                          input_file, global, openclose,       &
+                                          variable_name )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string !< name of the attribue
+       
+       CHARACTER(LEN=*) ::  input_file    !< name of input file
+       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
+       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
+       
+       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
+       INTEGER(iwp) ::  val      !< value of the attribute
+       
+       LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
+
+#if defined ( __netcdf )
+!
+!--    Open file in read-only mode
+       IF ( openclose == 'open' )  THEN
+          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
+                                  id_mod )
+       ENDIF
+!
+!--    Read global attribute
+       IF ( global )  THEN
+          CALL get_attribute( id_mod, search_string, val, global )
+!
+!--    Read variable attribute
+       ELSE
+          CALL get_attribute( id_mod, search_string, val, global, variable_name )
+       ENDIF
+!
+!--    Finally, close input file
+       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_att_int
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read a global real attribute 
+!------------------------------------------------------------------------------!
+    SUBROUTINE netcdf_data_input_att_real( val, search_string, id_mod,         &
+                                           input_file, global, openclose,      &
+                                           variable_name )
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*) ::  search_string !< name of the attribue
+       
+       CHARACTER(LEN=*) ::  input_file    !< name of input file
+       CHARACTER(LEN=*) ::  openclose     !< string indicating whether NetCDF needs to be opend or closed
+       CHARACTER(LEN=*) ::  variable_name !< string indicating whether NetCDF needs to be opend or closed 
+       
+       INTEGER(iwp) ::  id_mod            !< NetCDF id of input file
+       
+       LOGICAL ::  global                 !< flag indicating a global or a variable's attribute
+       
+       REAL(wp) ::  val                   !< value of the attribute
+
+#if defined ( __netcdf )
+!
+!--    Open file in read-only mode
+       IF ( openclose == 'open' )  THEN
+          CALL open_read_file( TRIM( input_file ) // TRIM( coupling_char ), &
+                                  id_mod )
+       ENDIF
+!
+!--    Read global attribute
+       IF ( global )  THEN
+          CALL get_attribute( id_mod, search_string, val, global )
+!
+!--    Read variable attribute
+       ELSE
+          CALL get_attribute( id_mod, search_string, val, global, variable_name )
+       ENDIF
+!
+!--    Finally, close input file
+       IF ( openclose == 'close' )  CALL close_input_file( id_mod )
+#endif            
+
+    END SUBROUTINE netcdf_data_input_att_real
 
 !------------------------------------------------------------------------------!
@@ -1022,5 +1250,5 @@
           !                                  the other dimensions depend on the mode of the emissions or on the presence of specific components
           !nspecies
-          CALL get_dimension_length( id_emis, emt_att%nspec, 'nspecies' )
+          CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nspec, 'nspecies' )
 
  
@@ -1051,5 +1279,5 @@
  
              !number of categories
-             CALL get_dimension_length( id_emis, emt_att%ncat, 'ncat' )
+             CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%ncat, 'ncat' )
 
              !-- Read EMISSION CATEGORIES INDEX
@@ -1064,5 +1292,5 @@
                 IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
                    !Allocate Array
-                   CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
+                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
                    ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
                    !Read-in Variable
@@ -1081,5 +1309,5 @@
                 !-- EMISSION_PM_NAME (1-DIMENSIONAL)
                 IF (TRIM(emt_att%species_name(ispec))=="PM" .OR. TRIM(emt_att%species_name(ispec))=="pm") THEN
-                   CALL get_dimension_length( id_emis, emt_att%npm, 'npm' )
+                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%npm, 'npm' )
                    ALLOCATE(emt_att%pm_name(1:emt_att%npm))
                    !Read-in Variable
@@ -1118,5 +1346,5 @@
              IF (TRIM(time_fac_type)=="HOUR" .OR. TRIM(time_fac_type)=="hour") THEN
                 !-- Allocate Array
-                CALL get_dimension_length( id_emis, emt_att%nhoursyear, 'nhoursyear' )
+                CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nhoursyear, 'nhoursyear' )
                 ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear))
                 !Read-in Variable
@@ -1127,5 +1355,5 @@
              ELSE IF (TRIM(time_fac_type) == "MDH" .OR. TRIM(time_fac_type) == "mdh") THEN
                 !-- Allocate Array
-                CALL get_dimension_length( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
+                CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
                 ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour))
                 !-- Read-in Variable
@@ -1176,5 +1404,5 @@
                 IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
                    !Allocate Array
-                   CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
+                   CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
                    ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
                    !Read-in Variable
@@ -1194,5 +1422,5 @@
              !-- EMISSION_VALUES (4-DIMENSIONAL)
              !Calculate temporal dimension length
-             CALL get_dimension_length( id_emis, emt_att%dt_emission, 'time' )   
+             CALL netcdf_data_input_get_dimension_length( id_emis, emt_att%dt_emission, 'time' )   
          
 
@@ -1294,6 +1522,7 @@
 !
 !--          Inquire number of vertical vegetation layer
-             CALL get_dimension_length( id_surf, leaf_area_density_f%nz,       &
-                                        'zlad' )
+             CALL netcdf_data_input_get_dimension_length( id_surf,             &
+                                                 leaf_area_density_f%nz,       &
+                                                 'zlad' )
 !
 !--          Allocate variable for leaf-area density
@@ -1318,6 +1547,7 @@
 !
 !--          Inquire number of vertical vegetation layer
-             CALL get_dimension_length( id_surf, basal_area_density_f%nz,      &
-                                        'zlad' )
+             CALL netcdf_data_input_get_dimension_length( id_surf,             &
+                                                 basal_area_density_f%nz,      &
+                                                 'zlad' )
 !
 !--          Allocate variable
@@ -1341,7 +1571,7 @@
 !
 !--          Inquire number of vertical soil layers
-             CALL get_dimension_length( id_surf,                               &
-                                        root_area_density_lad_f%nz,            &
-                                        'zsoil' )
+             CALL netcdf_data_input_get_dimension_length( id_surf,             &
+                                                   root_area_density_lad_f%nz, &
+                                                  'zsoil' )
 !
 !--          Allocate variable
@@ -1427,5 +1657,6 @@
 !
 !--          Obtain number of soil layers from file.
-             CALL get_dimension_length( id_surf, nz_soil, 'zsoil' )
+             CALL netcdf_data_input_get_dimension_length( id_surf, nz_soil,    &
+                                                          'zsoil' )
 
              ALLOCATE ( soil_type_f%var_3d(0:nz_soil,nys:nyn,nxl:nxr) )
@@ -1479,7 +1710,7 @@
 !
 !--       Inquire number of surface fractions
-          CALL get_dimension_length( id_surf,                                  &
-                                     surface_fraction_f%nf,                    &
-                                     'nsurface_fraction' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       surface_fraction_f%nf,  &
+                                                       'nsurface_fraction' )
 !
 !--       Allocate dimension array and input array for surface fractions
@@ -1508,7 +1739,7 @@
 !
 !--       Inquire number of building parameters
-          CALL get_dimension_length( id_surf,                                  &
-                                     building_pars_f%np,                       &
-                                     'nbuilding_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       building_pars_f%np,     &
+                                                       'nbuilding_pars' )
 !
 !--       Allocate dimension array and input array for building parameters
@@ -1551,6 +1782,7 @@
 !
 !--       Inquire number of albedo parameters
-          CALL get_dimension_length( id_surf, albedo_pars_f%np,                &
-                                     'nalbedo_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       albedo_pars_f%np,       &
+                                                       'nalbedo_pars' )
 !
 !--       Allocate dimension array and input array for albedo parameters
@@ -1578,6 +1810,7 @@
 !
 !--       Inquire number of pavement parameters
-          CALL get_dimension_length( id_surf, pavement_pars_f%np,              &
-                                     'npavement_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       pavement_pars_f%np,     &
+                                                       'npavement_pars' )
 !
 !--       Allocate dimension array and input array for pavement parameters
@@ -1606,12 +1839,12 @@
 !
 !--       Inquire number of parameters
-          CALL get_dimension_length( id_surf,                                  &
-                                     pavement_subsurface_pars_f%np,            &
-                                     'npavement_subsurface_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                pavement_subsurface_pars_f%np, &
+                                               'npavement_subsurface_pars' )
 !
 !--       Inquire number of soil layers
-          CALL get_dimension_length( id_surf,                                  &
-                                     pavement_subsurface_pars_f%nz,            &
-                                     'zsoil' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                pavement_subsurface_pars_f%nz, &
+                                                'zsoil' )
 !
 !--       Allocate dimension array and input array for pavement parameters
@@ -1646,6 +1879,7 @@
 !
 !--       Inquire number of vegetation parameters
-          CALL get_dimension_length( id_surf, vegetation_pars_f%np,            &
-                                     'nvegetation_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       vegetation_pars_f%np,   &
+                                                       'nvegetation_pars' )
 !
 !--       Allocate dimension array and input array for surface fractions
@@ -1679,7 +1913,7 @@
 !
 !--       Inquire number of soil parameters
-          CALL get_dimension_length( id_surf,                                  &
-                                     soil_pars_f%np,                           &
-                                     'nsoil_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       soil_pars_f%np,         &
+                                                       'nsoil_pars' )
 !
 !--       Read parameters array
@@ -1691,5 +1925,7 @@
 !--       soil layers, allocate memory and read the respective dimension
           IF ( soil_pars_f%lod == 2 )  THEN
-             CALL get_dimension_length( id_surf, soil_pars_f%nz, 'zsoil' )
+             CALL netcdf_data_input_get_dimension_length( id_surf,             &
+                                                          soil_pars_f%nz,      &
+                                                          'zsoil' )
 
              ALLOCATE( soil_pars_f%layers(0:soil_pars_f%nz-1) )
@@ -1730,7 +1966,7 @@
 !
 !--       Inquire number of water parameters
-          CALL get_dimension_length( id_surf,                                  &
-                                     water_pars_f%np,                          &
-                                     'nwater_pars' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                       water_pars_f%np,        &
+                                                       'nwater_pars' )
 !
 !--       Allocate dimension array and input array for water parameters
@@ -1756,6 +1992,7 @@
 !
 !--       Obtain number of soil layers from file and allocate variable
-          CALL get_dimension_length( id_surf, root_area_density_lsm_f%nz,      &
-                                     'zsoil' )
+          CALL netcdf_data_input_get_dimension_length( id_surf,                &
+                                                   root_area_density_lsm_f%nz, &
+                                                   'zsoil' )
           ALLOCATE( root_area_density_lsm_f%var                                &
                                         (0:root_area_density_lsm_f%nz-1,       &
@@ -2309,6 +2546,6 @@
 !
 !--       Read x, y - dimensions. Only required for consistency checks.
-          CALL get_dimension_length( id_topo, dim_static%nx, 'x' )
-          CALL get_dimension_length( id_topo, dim_static%ny, 'y' )
+          CALL netcdf_data_input_get_dimension_length( id_topo, dim_static%nx, 'x' )
+          CALL netcdf_data_input_get_dimension_length( id_topo, dim_static%ny, 'y' )
           ALLOCATE( dim_static%x(0:dim_static%nx-1) )
           ALLOCATE( dim_static%y(0:dim_static%ny-1) )
@@ -2371,5 +2608,6 @@
                                  .FALSE., 'buildings_3d' )
 
-             CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
+             CALL netcdf_data_input_get_dimension_length( id_topo,             &
+                                                          buildings_f%nz, 'z' )
 !
 !--          Read 3D buildings
@@ -2609,13 +2847,13 @@
 !
 !--    Read vertical dimension of scalar und w grid. 
-       CALL get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
-       CALL get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
 !
 !--    Read also the horizontal dimensions. These are used just used fo
 !--    checking the compatibility with the PALM grid before reading.
-       CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
-       CALL get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
-       CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
-       CALL get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
 
 !
@@ -3065,5 +3303,6 @@
 !--    Read vertical dimension for soil depth. 
        IF ( check_existence( var_names, 'zsoil' ) )                            &
-          CALL get_dimension_length( id_dynamic, init_3d%nzs, 'zsoil' )
+          CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nzs,&
+                                                       'zsoil' )
 !
 !--    Read also the horizontal dimensions required for soil initialization.
@@ -3071,6 +3310,8 @@
 !--    these data is already available, but will be read again for the sake
 !--    of clearness. 
-       CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
-       CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%nx,    &
+                                                    'x'  )
+       CALL netcdf_data_input_get_dimension_length( id_dynamic, init_3d%ny,    &
+                                                    'y'  )
 !
 !--    Check for correct horizontal and vertical dimension. Please note,
@@ -3207,5 +3448,6 @@
 !
 !--       Read time dimension, allocate memory and finally read time array
-          CALL get_dimension_length( id_dynamic, nest_offl%nt, 'time' )
+          CALL netcdf_data_input_get_dimension_length( id_dynamic,             &
+                                                       nest_offl%nt, 'time' )
 
           IF ( check_existence( nest_offl%var_names, 'time' ) )  THEN
@@ -3215,6 +3457,8 @@
 !
 !--       Read vertical dimension of scalar und w grid
-          CALL get_dimension_length( id_dynamic, nest_offl%nzu, 'z' )
-          CALL get_dimension_length( id_dynamic, nest_offl%nzw, 'zw' )
+          CALL netcdf_data_input_get_dimension_length( id_dynamic,             &
+                                                       nest_offl%nzu, 'z' )
+          CALL netcdf_data_input_get_dimension_length( id_dynamic,             &
+                                                       nest_offl%nzw, 'zw' )
 
           IF ( check_existence( nest_offl%var_names, 'z' ) )  THEN
@@ -4538,5 +4782,6 @@
 !> Get dimension array for a given dimension
 !------------------------------------------------------------------------------!
-     SUBROUTINE get_dimension_length( id, dim_len, variable_name )
+     SUBROUTINE netcdf_data_input_get_dimension_length( id, dim_len,           &
+                                                        variable_name )
 #if defined( __netcdf )
 
@@ -4555,87 +4800,138 @@
 !--    First, inquire dimension ID
        nc_stat = NF90_INQ_DIMID( id, TRIM( variable_name ), id_dim )
-       CALL handle_error( 'get_dimension_length', 526, variable_name )
+       CALL handle_error( 'netcdf_data_input_get_dimension_length', 526,       &
+                          variable_name )
 !
 !--    Inquire dimension length
        nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, dum, LEN = dim_len )
-       CALL handle_error( 'get_dimension_length', 526, variable_name )
+       CALL handle_error( 'netcdf_data_input_get_dimension_length', 526,       &
+                          variable_name )
 
 #endif
-    END SUBROUTINE get_dimension_length
+    END SUBROUTINE netcdf_data_input_get_dimension_length
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Routine for reading-in a character string from the chem emissions netcdf input file.  
-!------------------------------------------------------------------------------!
- 
- SUBROUTINE get_variable_string( id, variable_name, var_string, names_number)
+!> Routine for reading-in a character string from the chem emissions netcdf 
+!> input file.  
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_variable_string( id, variable_name, var_string, names_number)
 #if defined( __netcdf )
 
-    USE indices
-    USE pegrid
-
-    IMPLICIT NONE
-
-    CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  :: var_string
-
-    CHARACTER(LEN=*)                                              :: variable_name          !> variable name 
-
-    CHARACTER (LEN=1), ALLOCATABLE, DIMENSION(:,:)                :: tmp_var_string         !> variable to be read
-
-
-    INTEGER(iwp), INTENT(IN)                                      :: id                     !> file id
-
-    INTEGER(iwp), INTENT(IN)                                      :: names_number           !> number of names
-
-    INTEGER(iwp)                                                  :: id_var                 !> variable id
-
-    INTEGER(iwp)                                                  :: i,j                    !> index to go through the length of the dimensions
-
-    INTEGER(iwp)                                                  :: max_string_length=25   !> this is both the maximum length of a name, but also  
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  :: var_string
+
+       CHARACTER(LEN=*)                                              :: variable_name          !> variable name 
+
+       CHARACTER (LEN=1), ALLOCATABLE, DIMENSION(:,:)                :: tmp_var_string         !> variable to be read
+
+
+       INTEGER(iwp), INTENT(IN)                                      :: id                     !> file id
+
+       INTEGER(iwp), INTENT(IN)                                      :: names_number           !> number of names
+
+       INTEGER(iwp)                                                  :: id_var                 !> variable id
+
+       INTEGER(iwp)                                                  :: i,j                    !> index to go through the length of the dimensions
+
+       INTEGER(iwp)                                                  :: max_string_length=25   !> this is both the maximum length of a name, but also  
                                                                                             ! the number of the components of the first dimensions
                                                                                             ! (rows)
 
 
-    ALLOCATE(tmp_var_string(max_string_length,names_number))
-
-    ALLOCATE(var_string(names_number))
+       ALLOCATE(tmp_var_string(max_string_length,names_number))
+
+       ALLOCATE(var_string(names_number))
 
     !-- Inquire variable id
-    nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 
 
     !-- Get variable
     !-- Start cycle over the emission species
-    DO i = 1, names_number 
+       DO i = 1, names_number 
        !-- read the first letter of each component
-       nc_stat = NF90_GET_VAR( id, id_var, var_string(i), start = (/ 1,i /), &
+          nc_stat = NF90_GET_VAR( id, id_var, var_string(i), start = (/ 1,i /), &
                                  count = (/ 1,1 /) )
-       CALL handle_error( 'get_variable_string', 701 )
+          CALL handle_error( 'get_variable_string', 701 )
 
        !-- Start cycle over charachters
-       DO j = 1, max_string_length
+          DO j = 1, max_string_length
                        
           !-- read the rest of the components of the name
-          nc_stat = NF90_GET_VAR( id, id_var, tmp_var_string(j,i), start = (/ j,i /),&
+             nc_stat = NF90_GET_VAR( id, id_var, tmp_var_string(j,i), start = (/ j,i /),&
                                      count = (/ 1,1 /) )
-          CALL handle_error( 'get_variable_string', 702 )
-
-          IF ( iachar(tmp_var_string(j,i) ) == 0 ) THEN
-               tmp_var_string(j,i)=''
-          ENDIF
-
-          IF ( j>1 ) THEN
-          !-- Concatenate first letter of the name and the others
-             var_string(i)=TRIM(var_string(i)) // TRIM(tmp_var_string(j,i))
-
-          ENDIF 
-       ENDDO 
-    ENDDO
+             CALL handle_error( 'get_variable_string', 702 )
+
+             IF ( iachar(tmp_var_string(j,i) ) == 0 ) THEN
+                  tmp_var_string(j,i)=''
+             ENDIF
+
+             IF ( j>1 ) THEN
+             !-- Concatenate first letter of the name and the others
+                var_string(i)=TRIM(var_string(i)) // TRIM(tmp_var_string(j,i))
+
+             ENDIF 
+          ENDDO 
+       ENDDO
 
 #endif
- END SUBROUTINE get_variable_string
-
-
+    END SUBROUTINE get_variable_string
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads a character variable in a 1D array
+!------------------------------------------------------------------------------!
+     SUBROUTINE get_variable_1d_char( id, variable_name, var )
+
+       USE pegrid
+
+       IMPLICIT NONE
+
+       CHARACTER(LEN=*)            ::  variable_name          !< variable name
+       CHARACTER(LEN=*), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
+
+       INTEGER(iwp)                ::  i                !< running index over variable dimension
+       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
+       INTEGER(iwp)                ::  id_var           !< dimension id
+       
+       INTEGER(iwp), DIMENSION(2)  ::  dimid            !< dimension IDs
+       INTEGER(iwp), DIMENSION(2)  ::  dimsize          !< dimension size
+
+#if defined( __netcdf )
+
+!
+!--    First, inquire variable ID
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
+!
+!--    Inquire dimension IDs
+       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, dimids = dimid(1:2) )
+       CALL handle_error( 'get_variable_1d_char', 527, variable_name )
+!
+!--    Read dimesnion length
+       nc_stat = NF90_INQUIRE_DIMENSION( id, dimid(1), LEN = dimsize(1) )
+       nc_stat = NF90_INQUIRE_DIMENSION( id, dimid(2), LEN = dimsize(2) )
+       
+!
+!--    Read character array. Note, each element is read individually, in order
+!--    to better separate single strings. 
+       DO  i = 1, dimsize(2)
+          nc_stat = NF90_GET_VAR( id, id_var, var(i),                          &
+                                  start = (/ 1, i /),                          &
+                                  count = (/ dimsize(1), 1 /) )
+          CALL handle_error( 'get_variable_1d_char', 527, variable_name )
+       ENDDO     
+                          
+#endif
+    END SUBROUTINE get_variable_1d_char
+
+    
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4803,15 +5099,14 @@
        ELSE
 !
-!--    Allocate temporary variable according to memory access on file. 
-       ALLOCATE( tmp(is:ie,js:je) )
-!
-!--    Get variable
-          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
-                               start = (/ is+1,      js+1 /),                  &
-                               count = (/ ie-is + 1, je-js+1 /) )   
-                               
+!--       Allocate temporary variable according to memory access on file. 
+          ALLOCATE( tmp(is:ie,js:je) )
+!
+!--       Get variable
+          nc_stat = NF90_GET_VAR( id, id_var, tmp,                             &
+                                  start = (/ is+1,      js+1 /),               &
+                                  count = (/ ie-is + 1, je-js+1 /) )   
           CALL handle_error( 'get_variable_2d_real', 530, variable_name )
 !
-!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+!--       Resort data. Please note, dimension subscripts of var all start at 1. 
           DO  i = is, ie 
              DO  j = js, je 
@@ -5089,85 +5384,4 @@
 #endif
     END SUBROUTINE get_variable_3d_real
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Reads a 3D float array from file.
-!------------------------------------------------------------------------------!
-!     SUBROUTINE get_variable_3d_real_v( id, variable_name, is, ie, js, je, var )
-! 
-!        USE indices
-!        USE pegrid
-! 
-!        IMPLICIT NONE
-! 
-!        CHARACTER(LEN=*)              ::  variable_name   !< variable name
-! 
-!        INTEGER(iwp), INTENT(IN)      ::  is,ie           !< index range along x direction
-!        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
-!        INTEGER(iwp)                  ::  id_var          !< variable id
-!        INTEGER(iwp), INTENT(IN)      ::  js,je           !< index range along y direction
-!        INTEGER(iwp)                  ::  n3              !< number of data-points along 3rd dimension
-! 
-!        INTEGER(iwp)                  ::  i,j,k
-!        INTEGER(iwp), DIMENSION(3)    ::  id_dim
-! 
-!        REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var         !< variable to be read
-! #if defined( __netcdf )
-! !
-! !--    Inside the ...static NetCDF files, the array is stored as float.
-! !--    Therefore single precision is sufficiant for the temporary array
-! 
-!        REAL(sp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_var     !< temporary array to read NetCDF data in i,j,k direction
-! 
-! !kk    Please check, if it is time consuming to do the inquire every time
-! !
-! !--    Inquire variable id
-!        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-! !
-! !--    Get length of third dimension, required for the count parameter.
-! !--    Therefore, first inquired dimension ids
-!        nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
-!        nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
-! 
-! !
-! !--    Check for collective read-operation and set respective NetCDF flags if
-! !--    required. 
-!        IF ( collective_read )  THEN
-! #if defined( __netcdf4_parallel )
-!           nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
-! #endif
-!        ENDIF
-! 
-! !
-! !--    Allocate temporary array ro read NetCDF data in i,j,k direction
-! 
-!        ALLOCATE(tmp_var(is:ie,js:je,n3))
-! !
-! !--    Get variable
-! !--    Read complete local 3-D array in oone call
-! 
-!        nc_stat = NF90_GET_VAR( id, id_var, tmp_var,                            &
-!                                start = (/ is+1, js+1, 1 /),                    &
-!                                count = (/ ie-is+1, je-js+1, n3 /) )
-! 
-!        CALL handle_error( 'get_variable_3d_real', 532 )
-! 
-! !
-! !--    Resort data in k,j,i direction
-! 
-!        DO i=is,ie
-!           DO j=js,je
-!              DO K=1,n3
-!                 var (k,j-js+1,i-is+1) = tmp_var(i,j,k)
-!              END DO
-!           END DO
-!        END DO
-! 
-!        DEALLOCATE(tmp_var)
-! 
-! #endif
-!     END SUBROUTINE get_variable_3d_real_v
-
 
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3469)
+++ palm/trunk/SOURCE/parin.f90	(revision 3472)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Add virtual measurement module
 ! 
 ! Former revisions:
@@ -558,4 +558,7 @@
     USE vertical_nesting_mod,                                                  &
         ONLY:  vnest_start_time
+        
+    USE virtual_measurement_mod,                                               &
+        ONLY:  vm_parin
 
     USE wind_turbine_model_mod,                                                &
@@ -926,4 +929,5 @@
           CALL im_parin
           CALL salsa_parin
+          CALL vm_parin
 !
 !--       Read user-defined variables
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3469)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3472)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! new module for virtual measurements introduced
 ! 
 ! Former revisions:
@@ -458,5 +458,6 @@
                use_initial_profile_as_reference,                               &
                use_single_reference_value, uv_exposure, u_gtrans, v_gtrans,    &
-               virtual_flight, wind_turbine, ws_scheme_mom, ws_scheme_sca
+               virtual_flight, virtual_measurement, wind_turbine,              &
+               ws_scheme_mom, ws_scheme_sca
 
     USE cpulog,                                                                &
@@ -581,4 +582,7 @@
                vnest_deallocate, vnest_init, vnest_init_fine,                  &
                vnest_start_time
+               
+    USE virtual_measurement_mod,                                               &
+        ONLY:  vm_sampling, vm_time_start
 
     IMPLICIT NONE
@@ -1458,4 +1462,8 @@
           CALL data_output_flight
        ENDIF
+!
+!--    Take virtual measurements
+       IF ( virtual_measurement  .AND.                                         &
+            vm_time_start <= time_since_reference_point )  CALL vm_sampling
 
 !
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3472)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3472)
@@ -0,0 +1,769 @@
+!> @virtual_measurement_mod.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 2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! Initial revision
+! 
+! Former revisions:
+! -----------------
+! $Id$
+!
+! Authors:
+! --------
+! @author Matthias Suehring
+! 
+!
+!
+! Description:
+! ------------
+!> 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. 
+!> More precisely, coordinates and measured quanties will be read from a
+!> NetCDF file which contains all required information. In the model, 
+!> the same quantities (as long as all the required components are switched-on)
+!> will be sampled at the respective positions and output into an extra file,
+!> which allows for straight-forward comparison of model results with 
+!> observations. 
+!------------------------------------------------------------------------------!
+ MODULE virtual_measurement_mod
+
+
+    USE arrays_3d,                                                             &
+        ONLY:  q, pt, u, v, w, zu, zw
+
+    USE control_parameters,                                                    &
+        ONLY:  dz, message_string, virtual_measurement
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point
+        
+    USE grid_variables,                                                        &
+        ONLY:  dx, dy
+
+    USE indices,                                                               &
+        ONLY:  nzb, nzt, nxl, nxr, nys, nyn
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    TYPE virt_mea
+    
+       CHARACTER(LEN=100)  ::  feature_type  !< type of the measurement
+       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    !< 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), DIMENSION(:), ALLOCATABLE ::  dim_t !< number observations individual for each trajectory or station that are no _FillValues
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ngp   !< number of grid points where observations for a site took place, 
+                                                         !<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
+           
+       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
+       
+       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 ::  xmea   !< measurement x-position in absolute UTM coordinates  
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ymea   !< measurement y-position in absolute UTM coordinates  
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zmea   !< measurement z-position in height above ground level 
+       
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  measured_vars  !< measured variables 
+       
+    END TYPE virt_mea
+
+    CHARACTER(LEN=5)  ::  char_eutm = "E_UTM"                      !< dimension name for UTM coordinate easting
+    CHARACTER(LEN=11) ::  char_feature = "featureType"             !< attribute name for feature type
+    CHARACTER(LEN=10) ::  char_fillvalue = "_FillValue"            !< variable attribute name for _FillValue 
+    CHARACTER(LEN=18) ::  char_mv = "measured_variables"           !< variable name for the array with the measured variable names
+    CHARACTER(LEN=5)  ::  char_nutm = "N_UTM"                      !< dimension name for UTM coordinate northing
+    CHARACTER(LEN=18) ::  char_numstations = "number_of_stations"  !< attribute name for number of stations
+    CHARACTER(LEN=8)  ::  char_origx = "origin_x"                  !< attribute name for station coordinate in x
+    CHARACTER(LEN=8)  ::  char_origy = "origin_y"                  !< attribute name for station coordinate in y
+    CHARACTER(LEN=4)  ::  char_site = "site"                       !< attribute name for site name
+    CHARACTER(LEN=19) ::  char_zag = "height_above_ground"         !< attribute name for height above ground variable
+    CHARACTER(LEN=10) ::  type_ts   = 'timeSeries'                 !< name of stationary point measurements
+    CHARACTER(LEN=10) ::  type_traj = 'trajectory'                 !< name of line measurements
+    CHARACTER(LEN=17) ::  type_tspr = 'timeSeriesProfile'          !< name of stationary profile measurements
+    
+    CHARACTER(LEN=10), DIMENSION(1:53), PARAMETER ::  list_allowed_variables = & !< variables that can be sampled in PALM
+       (/ 'hfls      ',  & ! surface latent heat flux (W/m2)
+          'hfss      ',  & ! surface sensible heat flux (W/m2)
+          'hur       ',  & ! relative humidity (-)
+          'hus       ',  & ! specific humidity (g/kg)
+          'haa       ',  & ! absolute atmospheric humidity (kg/m3)
+          'mcpm1     ',  & ! mass concentration of PM1 (kg/m3)
+          'mcpm2p5   ',  & ! mass concentration of PM2.5 (kg/m3)
+          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
+          'mcpm10    ',  & ! mass concentration of PM10 (kg/m3)
+          'mcco      ',  & ! mass concentration of CO (kg/m3)
+          'mcco2     ',  & ! mass concentration of CO2 (kg/m3)
+          'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
+          'ncaa      ',  & ! number concentation of particles (1/m3)
+          'mfco2     ',  & ! mole fraction of CO (mol/mol)
+          'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
+          'mfch4     ',  & ! mole fraction of methane (mol/mol)
+          'mfnh3     ',  & ! mole fraction of amonia (mol/mol)
+          'mfno      ',  & ! mole fraction of nitrogen monoxide (mol/mol)
+          'mfno2     ',  & ! mole fraction of nitrogen dioxide (mol/mol)
+          'mfso2     ',  & ! mole fraction of sulfur dioxide (mol/mol)
+          'mfh20     ',  & ! mole fraction of water (mol/mol)
+          'plev      ',  & ! ? air pressure - hydrostaic + perturbation?
+          'rlds      ',  & ! surface downward longwave flux  (W/m2)
+          'rlus      ',  & ! surface upward longwave flux (W/m2)
+          'rsds      ',  & ! surface downward shortwave flux (W/m2)
+          'rsus      ',  & ! surface upward shortwave flux (W/m2)
+          'ta        ',  & ! air temperature (degree C)
+          't_va      ',  & ! virtual accoustic temperature (K)
+          'theta     ',  & ! potential temperature (K)
+          'tro3      ',  & ! mole fraction of ozone air (mol/mol)
+          'ts        ',  & ! scaling parameter of temperature (K)
+          'wspeed    ',  & ! ? wind speed - horizontal?
+          'wdir      ',  & ! wind direction
+          'us        ',  & ! friction velocity
+          'msoil     ',  & ! ? soil moisture - which depth?  
+          'tsoil     ',  & ! ? soil temperature - which depth?                                                                
+          'u         ',  & ! u-component
+          'ua        ',  & ! eastward wind (is there any difference to u?)
+          'uw        ',  & ! ? vertical momentum flux - total ?
+          'utheta    ',  & ! ? horizontal heat flux - total ? 
+          'uv        ',  & ! upward-northward horizontal momentum flux
+          'v         ',  & ! v-component
+          'va        ',  & ! northward wind (is there any difference to v?)
+          'vw        ',  & ! ? vertical momentum flux - total ?
+          'vtheta    ',  & ! ? horizontal heat flux - total ? 
+          'w         ',  & ! w-component
+          'wtheta    ',  & ! ? vertical heat flux - total ? 
+          'rld       ',  & ! downward longwave radiative flux (W/m2)
+          'rlu       ',  & ! upnward longwave radiative flux (W/m2)
+          'rsd       ',  & ! downward shortwave radiative flux (W/m2)
+          'rsu       ',  & ! upward shortwave radiative flux (W/m2)
+          'rsddif    ',  & ! downward shortwave diffuse radiative flux (W/m2)
+          'rnds      '   & ! surface net downward radiative flux (W/m2)
+       /)
+                                                             
+    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 ::  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_mea ), DIMENSION(:), ALLOCATABLE ::  vmea !< virtual measurement data structure
+    
+    
+    INTERFACE vm_check_parameters
+       MODULE PROCEDURE vm_check_parameters
+    END INTERFACE vm_check_parameters
+    
+    INTERFACE vm_init
+       MODULE PROCEDURE vm_init
+    END INTERFACE vm_init
+    
+    INTERFACE vm_parin
+       MODULE PROCEDURE vm_parin
+    END INTERFACE vm_parin
+    
+    INTERFACE vm_sampling
+       MODULE PROCEDURE vm_sampling
+    END INTERFACE vm_sampling
+
+    SAVE
+
+    PRIVATE
+
+!
+!-- Public interfaces
+    PUBLIC  vm_check_parameters, vm_init, vm_parin, vm_sampling
+
+!
+!-- Public variables
+    PUBLIC  vmea, vm_time_start
+
+ CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters for virtual measurement module
+!------------------------------------------------------------------------------!
+ SUBROUTINE vm_check_parameters
+
+    USE control_parameters,                                                    &
+        ONLY:  message_string, virtual_measurement
+ 
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  input_pids_static
+        
+    IMPLICIT NONE
+    
+!
+!-- In case virtual measurements are taken, a static input file is required.
+!-- This is because UTM coordinates for the PALM domain origin are required
+!-- for correct mapping of the measurements. 
+!-- 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 ' // &
+                        'file is mandatory.'
+       CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
+    ENDIF
+ 
+ END SUBROUTINE vm_check_parameters
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read namelist for the virtual measurement module
+!------------------------------------------------------------------------------!
+ SUBROUTINE vm_parin
+ 
+    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
+ 
+    NAMELIST /virtual_measurement_parameters/  use_virtual_measurement,        &
+                                               vm_time_start
+
+    line = ' '
+
+!
+!-- Try to find stg package
+    REWIND ( 11 )
+    line = ' '
+    DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
+       READ ( 11, '(A)', END=20 )  line
+    ENDDO
+    BACKSPACE ( 11 )
+
+!
+!-- Read namelist
+    READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
+
+!
+!-- Set flag that indicates that the virtual measurement module is switched on
+    IF ( use_virtual_measurement )  virtual_measurement = .TRUE.
+    
+    GOTO 20
+
+ 10 BACKSPACE( 11 )
+    READ( 11 , '(A)') line
+    CALL parin_fail_message( 'virtual_measurement_parameters', line )
+
+ 20 CONTINUE
+ 
+ END SUBROUTINE vm_parin
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize virtual measurements: read coordiante arrays and measured 
+!> variables, set indicies indicating the measurement points, read further
+!> attributes, etc..
+!------------------------------------------------------------------------------!
+ SUBROUTINE vm_init
+
+    USE arrays_3d,                                                             &
+        ONLY:  zu, zw
+ 
+    USE control_parameters,                                                    &
+        ONLY:  message_string
+        
+    USE grid_variables,                                                        &
+        ONLY:  ddx, ddy, dx, dy
+        
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys
+ 
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  init_model, input_file_vm,                                      &
+               netcdf_data_input_get_dimension_length,                         &
+               netcdf_data_input_att, netcdf_data_input_var
+               
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index_ji
+        
+    IMPLICIT NONE
+    
+    CHARACTER(LEN=5)    ::  dum                !< dummy string indicate station id
+    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables_file = '' !< array with all measured variables read from NetCDF
+    CHARACTER(LEN=10), DIMENSION(50) ::  measured_variables      = '' !< dummy array with all measured variables that are allowed
+    
+    LOGICAL ::  on_pe !< flag indicating that the respective measurement coordinate is on subdomain
+    
+    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) ::  icov      !< index range where observations should be taken in x-direction
+    INTEGER(iwp) ::  ii        !< running index over all coordinate points of a measurement
+    INTEGER(iwp) ::  i_prev    !< grid index along x for UTM coordinate at previous observation time step 
+    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) ::  jcov      !< index range where observations should be taken in y-direction
+    INTEGER(iwp) ::  j_prev    !< grid index along y for UTM coordinate at previous observation time step 
+    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) ::  kcov      !< index range where observations should be taken in z-direction
+    INTEGER(iwp) ::  ks        !< grid index of real observation point of the respective station in z-direction
+    INTEGER(iwp) ::  k_prev    !< grid index along z for UTM coordinate at previous observation time step 
+    INTEGER(iwp) ::  ksurf     !< topography top index
+    INTEGER(iwp) ::  l         !< running index over all stations
+    INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
+    INTEGER(iwp) ::  ll        !< running index over all measured variables in file
+    INTEGER(iwp) ::  lll       !< running index over all allowed variables
+    INTEGER(iwp) ::  n         !< running index over trajectory coordinates
+    INTEGER(iwp) ::  ns        !< counter variable for number of observation points on subdomain
+    INTEGER(iwp) ::  t         !< running index over number of trajectories
+    
+    REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
+    REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
+    REAL(wp)     ::  fill_zag  !< _FillValue for height coordinate
+    
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  e_utm !< easting UTM coordinate, temporary variable
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_utm !< northing UTM coordinate, temporary variable,
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  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,   &
+                                global_attribute, 'open', '' )
+!
+!-- ALLOCATE data structure
+    ALLOCATE( vmea(1:nvm) )
+    
+!    print*, "nvm", nvm
+!
+!-- Read station coordinates and further attributes. 
+!-- Note all coordinates are in UTM coordinates.
+    DO  l = 1, nvm
+!
+!--    Determine suffix which contains the ID
+       IF( l < 10 )  THEN
+          WRITE( dum, '(I1)')  l
+       ELSEIF( l < 100 )  THEN
+          WRITE( dum, '(I2)')  l
+       ELSEIF( l < 1000 )  THEN
+          WRITE( dum, '(I3)')  l
+       ELSEIF( l < 10000 )  THEN
+          WRITE( dum, '(I4)')  l
+       ELSEIF( l < 100000 )  THEN
+          WRITE( dum, '(I5)')  l
+       ENDIF
+       
+       CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx            &
+                                   // TRIM( dum ), id_vm, '', global_attribute,&
+                                   '', '' )
+       CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy            &
+                                   // TRIM( dum ), 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 logicals depending on the type of the measurement
+       IF ( INDEX( vmea(l)%feature_type, type_tspr     ) /= 0 )  THEN
+          vmea(l)%timseries_profile = .TRUE.
+       ELSEIF ( INDEX( vmea(l)%feature_type, type_ts   ) /= 0 )  THEN
+          vmea(l)%timseries         = .TRUE.
+       ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 )  THEN
+          vmea(l)%trajectory        = .TRUE.
+       ELSE
+!
+!--       Give error message
+          message_string = 'Attribue featureType = ' //                        &
+                           TRIM( vmea(l)%feature_type ) //                     &
+                           ' is not allowed.' 
+          CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    Read string with all measured variables at this station
+       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
+       DO ll = 1, SIZE( measured_variables_file )
+          IF ( measured_variables_file(ll)(1:1) /= CHAR(0)  .AND.              &
+               measured_variables_file(ll)(1:1) /= ' ')  THEN
+!
+!--          Obtain character length of the character
+             len_char = 1
+             DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
+                 .AND.  measured_variables_file(ll)(len_char:len_char) /= ' ' )
+                len_char = len_char + 1
+             ENDDO
+             len_char = len_char - 1
+!
+!--          Now, compare the measured variable with the list of allowed 
+!--          variables. 
+             DO  lll= 1, SIZE( list_allowed_variables )
+                IF ( measured_variables_file(ll)(1:len_char) ==                &
+                     TRIM( list_allowed_variables(lll) ) )  THEN
+                   vmea(l)%nvar = vmea(l)%nvar + 1
+                   measured_variables(vmea(l)%nvar) =                          &
+                                       measured_variables_file(ll)(1:len_char)
+                ENDIF
+             ENDDO
+          ENDIF
+       ENDDO
+!
+!--    Allocate array for the measured variables names for the station l.
+       ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
+
+       DO  ll = 1, vmea(l)%nvar
+          vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
+       ENDDO
+       
+!        print*, "numvars", vmea(l)%nvar, vmea(l)%measured_vars_name(1:vmea(l)%nvar)
+!
+!--    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 ( 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.
+          IF ( vmea(l)%trajectory )  THEN
+!
+!--          For non-stationary measurements read the number of trajectories
+             CALL netcdf_data_input_get_dimension_length( id_vm,              &
+                                                          vmea(l)%ntraj,      &
+                                                          "traj" //           &
+                                                          TRIM( dum ) )
+             CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime,   &
+                                                          "ntime" //          &
+                                                          TRIM( dum ) )
+!
+!--       For stationary measurements the dimension for UTM coordinates is 1
+          ELSE
+             vmea(l)%ntraj  = 1
+             dim_ntime = 1
+          ENDIF
+          
+!
+!-        Allocate array which defines individual time frame for each 
+!--       trajectory or station
+          ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
+          ALLOCATE( vmea(l)%ngp(1:vmea(l)%ntraj)   )
+!
+!--       Allocate temporary arrays for UTM and height coordinates. Note, 
+!--       on file UTM coordinates might be 1D or 2D variables
+          ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
+          ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
+          ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime)  )
+!
+!--       Read _FillValue attributes
+          CALL netcdf_data_input_att( fill_eutm, char_fillvalue,               &
+                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      char_eutm // TRIM( dum ) )
+          CALL netcdf_data_input_att( fill_nutm, char_fillvalue,               &
+                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      char_nutm // TRIM( dum ) )
+          CALL netcdf_data_input_att( fill_zag, char_fillvalue,                &
+                                      id_vm, '', .NOT. global_attribute, '',   &
+                                      char_zag  // TRIM( dum ) )
+!
+!--       Read UTM and height coordinates coordinates for all trajectories and 
+!--       times. 
+          IF ( vmea(l)%trajectory )  THEN
+             CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm,  &
+                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
+             CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm,  &
+                                         0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
+             CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), 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
+!
+!--       Based on UTM coordinates, check if the measurement station or parts
+!--       of the trajectory is on subdomain. This case, setup grid index space
+!--       sample these quantities. 
+          ns = 0
+          DO  t = 1, vmea(l)%ntraj
+!
+!--          Determine the individual time coordinate length for each station and
+!--          trajectory. This is required as several stations and trajectories
+!--          are merged into one file but they do not have the same number of
+!--          points in time, hence, missing values may occur and cannot be 
+!--          processed further.
+             vmea(l)%dim_t(t) = 0
+             DO  n = 1, dim_ntime
+                IF ( e_utm(t,n) /= fill_eutm  .AND.                            &
+                     n_utm(t,n) /= fill_nutm  .AND.                            &
+                     z_ag(t,n)  /= fill_zag )  vmea(l)%dim_t(t) = n
+             ENDDO
+!             
+!--          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.
+!              write(9,*) l, "before Eutm", e_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_x
+!              write(9,*) l, "before Nutm", n_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_y
+             
+             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
+!              write(9,*) l, "Eutm", e_utm(t,1:vmea(l)%dim_t(t)) 
+!              write(9,*) l, "Nutm", n_utm(t,1:vmea(l)%dim_t(t)) 
+!              write(9,*)
+!
+!--          Compute grid indices relative to origin and check if these are
+!--          on the subdomain. Note, virtual measurements will be taken also 
+!--          at grid points surrounding the station, hence, check also for 
+!--          these grid points.
+             vmea(l)%ngp(t) = 0
+             k_prev = -999
+             j_prev = -999
+             i_prev = -999
+             DO  n = 1, vmea(l)%dim_t(t)
+                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
+                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
+!
+!--             Is the observation point on subdomain?
+                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
+                          js >= nys  .AND.  js <= nyn )
+!
+!--             If the measurement is on subdomain, determine vertical index
+!--             which refers to the observation height above ground level.
+                ks = k_prev
+                IF ( on_pe )  THEN
+!                    write(9,*) " I am onpe" 
+                   ksurf = get_topography_top_index_ji( js, is, 's' )
+                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
+                ENDIF
+!
+!--             Count the number of observation points in index space on 
+!--             subdomain. Only increment if grid indices are different from
+!--             the previous one. 
+                IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
+                                   ks /= k_prev )  THEN
+                   ns             = ns             + 1
+                   vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1
+                ENDIF
+
+!--             Store arrays for next iteration - avoid double counting 
+                i_prev = is
+                j_prev = js
+                k_prev = ks
+             ENDDO
+             
+          ENDDO
+
+!
+!--       Store number of observation points on subdomain and allocate index
+!--       arrays.
+          vmea(l)%ns = ns
+          
+          ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
+          ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
+          ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
+          
+!           print*, "Num ns: ", vmea(l)%ns, "per traj", vmea(l)%ngp(:)
+!
+!--       Repeat the prior loop and save the grid indices relevant for 
+!--       sampling.
+          ns = 0
+          DO  t = 1, vmea(l)%ntraj
+!
+!--          Compute grid indices relative to origin and check if these are
+!--          on the subdomain. Note, virtual measurements will be taken also 
+!--          at grid points surrounding the station, hence, check also for 
+!--          these grid points.
+             k_prev = -999
+             j_prev = -999
+             i_prev = -999
+             DO  n = 1, vmea(l)%dim_t(t)
+                is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
+                js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )             
+!
+!--             Is the observation point on subdomain?
+                on_pe = ( is >= nxl  .AND.  is <= nxr  .AND.                   &
+                          js >= nys  .AND.  js <= nyn )
+!
+!--             If the measurement is on subdomain, determine vertical index
+!--             which refers to the observation height above ground level.
+                ks = k_prev
+                IF ( on_pe )  THEN
+                   ksurf = get_topography_top_index_ji( js, is, 's' )
+                   ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
+                ENDIF
+!
+!--             Count the number of observation points in index space on 
+!--             subdomain. Only increment if grid indices are different from
+!--             the previous one. 
+                IF ( on_pe  .AND.  is /= i_prev  .AND.  js /= j_prev  .AND.    &
+                                   ks /= k_prev )  THEN
+                   ns             = ns + 1
+                   vmea(l)%i(ns)  = is
+                   vmea(l)%j(ns)  = js
+                   vmea(l)%k(ns)  = ks
+!                    write(9,*) TRIM( vmea(l)%feature_type), l, "ns", ns, "ijk", vmea(l)%i(ns), vmea(l)%j(ns), vmea(l)%k(ns)
+                ENDIF
+!
+!--             Store arrays for next iteration - avoid double counting 
+                i_prev = is
+                j_prev = js
+                k_prev = ks
+             ENDDO
+             
+          ENDDO
+!
+!--       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
+!
+!--       Deallocate temporary coordinate arrays
+          IF ( ALLOCATED( e_utm ) )  DEALLOCATE( e_utm )
+          IF ( ALLOCATED( n_utm ) )  DEALLOCATE( n_utm )
+          IF ( ALLOCATED( z_ag  ) )  DEALLOCATE( z_ag  )
+       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, '',              &
+                                global_attribute, 'close', '' )
+  END SUBROUTINE vm_init
+  
+  
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sampling of the actual quantities along the observation coordinates
+!------------------------------------------------------------------------------!
+  SUBROUTINE vm_sampling
+
+    USE arrays_3d !,                                                             &
+!         ONLY:  pt
+
+    USE surface_mod
+    
+     IMPLICIT NONE
+     
+     CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
+      
+     INTEGER(iwp) ::  l !< 
+     INTEGER(iwp) ::  m !<
+     INTEGER(iwp) ::  var !< 
+     
+     INTEGER(iwp) ::  mm, j, i
+
+!      write(9,*) "sampling"
+!      flush(9)
+!
+!--  Loop over all stations. For each possible variable loop over all 
+!--  observation points
+     DO  l = 1, nvm
+!
+!--     Loop over all measured variables. Please note, for the moment
+!--     the same indices for scalar and velocity components are used.
+!--     ToDo: Revise this later.
+!         DO  m = 1, vmea(l)%ns
+!            j = vmea(l)%j(m)
+!            i = vmea(l)%i(m)
+!            
+!            IF ( i >= nxl  .AND.  i <= nxr  .AND.                               &
+!                 j >= nys  .AND.  j <= nyn )  THEN
+!               IF ( surf_def_h(0)%start_index(j,i) <= &
+!                    surf_def_h(0)%end_index(j,i) )  THEN
+!                 
+!                  write(9,*) "sampled"
+!                  flush(9)
+!                  mm = surf_def_h(0)%end_index(j,i)
+!               
+!                  surf_def_h(0)%pt_surface(mm) = -99.0
+!               ENDIF
+!            ENDIF
+!         ENDDO
+!         DO  var = 1, vmea(l)%nvar
+!            trimvar = TRIM( vmea(l)%measured_vars_name(var) )
+!            
+!            IF ( TRIM( trimvar ) == 'theta' )  THEN           
+!               DO  m = 1, vmea(l)%ns
+!                  vmea(l)%measured_vars(var,m) = pt(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
+!               ENDDO
+!            ENDIF           
+!            IF ( TRIM( trimvar ) == 'w' )  THEN           
+!               DO  m = 1, vmea(l)%ns
+!                  vmea(l)%measured_vars(var,m) = w(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
+!               ENDDO
+!            ENDIF
+!            IF ( TRIM( trimvar ) == 'v' )  THEN           
+!               DO  m = 1, vmea(l)%ns
+!                  vmea(l)%measured_vars(var,m) = v(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
+!               ENDDO
+!            ENDIF
+!            IF ( TRIM( trimvar ) == 'u' )  THEN           
+!               DO  m = 1, vmea(l)%ns
+!                  vmea(l)%measured_vars(var,m) = u(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
+!               ENDDO
+!            ENDIF
+!         ENDDO
+
+     ENDDO
+     
+  END SUBROUTINE vm_sampling
+ 
+
+ END MODULE virtual_measurement_mod
