Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3518)
+++ palm/trunk/SOURCE/Makefile	(revision 3522)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Dependencies for virtual measurement module added
+# 
+# 3494 2018-11-06 14:51:27Z suehring
 # Surface output revised
 # 
@@ -1982,8 +1985,11 @@
 virtual_measurement_mod.o: \
 	cpulog_mod.o \
+	chemistry_model_mod.o \
+	chem_modules.o \
 	mod_kinds.o \
 	modules.o \
 	netcdf_data_input_mod.o \
-	netcdf_interface_mod.o
+	netcdf_interface_mod.o \
+	radiation_model_mod.o
 wind_turbine_model_mod.o: \
 	basic_constants_and_equations_mod.o \
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3518)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 3522)
@@ -25,16 +25,15 @@
 ! -----------------
 ! $Id$
+! Sampling of variables
+! 
+! 3494 2018-11-06 14:51:27Z suehring
 ! Bugfixing
 ! 
 ! 3473 2018-10-30 20:50:15Z suehring
 ! Initial revision
-! 
-! 3472 2018-10-30 20:43:50Z suehring
 !
 ! Authors:
 ! --------
-! @author Matthias Suehring and Klaus Ketelsen
-! 
-!
+! @author Matthias Suehring
 !
 ! Description:
@@ -49,13 +48,26 @@
 !> which allows for straight-forward comparison of model results with 
 !> observations. 
+!>
+!> @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 
+!> @note Fluxes are not processed
 !------------------------------------------------------------------------------!
  MODULE virtual_measurement_mod
 
-
     USE arrays_3d,                                                             &
         ONLY:  q, pt, u, v, w, zu, zw
 
+    USE chem_modules,                                                          &
+        ONLY:  nspec
+
+    USE chemistry_model_mod,                                                   &
+        ONLY:  chem_species
+        
     USE control_parameters,                                                    &
-        ONLY:  dz, message_string, virtual_measurement
+        ONLY:  air_chemistry, dz, humidity, neutral, message_string,           &
+               virtual_measurement
 
     USE cpulog,                                                                &
@@ -66,5 +78,5 @@
 
     USE indices,                                                               &
-        ONLY:  nzb, nzt, nxl, nxr, nys, nyn
+        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny
 
     USE kinds
@@ -80,11 +92,9 @@
        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) ::  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), 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
@@ -102,9 +112,5 @@
        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 
        
@@ -124,6 +130,7 @@
     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
+!
+!-- MS: List requires careful revision!
+    CHARACTER(LEN=10), DIMENSION(1:47), 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)
@@ -139,5 +146,5 @@
           'mcbcda    ',  & ! mass concentration of black carbon paritcles (kg/m3)
           'ncaa      ',  & ! number concentation of particles (1/m3)
-          'mfco2     ',  & ! mole fraction of CO (mol/mol)
+          'mfco      ',  & ! mole fraction of CO (mol/mol)
           'mfco2     ',  & ! mole fraction of CO2 (mol/mol)
           'mfch4     ',  & ! mole fraction of methane (mol/mol)
@@ -164,13 +171,7 @@
           '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)
@@ -307,7 +308,4 @@
     USE arrays_3d,                                                             &
         ONLY:  zu, zw
- 
-    USE control_parameters,                                                    &
-        ONLY:  message_string
         
     USE grid_variables,                                                        &
@@ -329,7 +327,5 @@
     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
+    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
@@ -338,17 +334,13 @@
     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) ::  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
-    INTEGER(iwp) ::  k_prev    !< grid index along z for UTM coordinate at previous observation time step 
     INTEGER(iwp) ::  ksurf     !< topography top index
+    INTEGER(iwp) ::  ku        !< upper vertical index of surrounding grid points of an observation coordinate
     INTEGER(iwp) ::  l         !< running index over all stations
     INTEGER(iwp) ::  len_char  !< character length of single measured variables without Null character
@@ -359,4 +351,9 @@
     INTEGER(iwp) ::  t         !< running index over number of trajectories
     
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  meas_flag !< mask array indicating measurement positions
+    
+    LOGICAL ::  chem_include !< flag indicating that chemical species is considered in modelled mechanism
+    LOGICAL ::  on_pe        !< flag indicating that the respective measurement coordinate is on subdomain
+    
     REAL(wp)     ::  fill_eutm !< _FillValue for coordinate array E_UTM
     REAL(wp)     ::  fill_nutm !< _FillValue for coordinate array N_UTM
@@ -370,9 +367,15 @@
     CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm,   &
                                 global_attribute, 'open', '' )
+                                
+!     write(9,*) "num stationi", nvm
+!     flush(9)
 !
 !-- ALLOCATE data structure
     ALLOCATE( vmea(1:nvm) )
-    
-!    print*, "nvm", nvm
+!
+!-- Allocate flag array
+    ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+    meas_flag = 0
+ 
 !
 !-- Read station coordinates and further attributes. 
@@ -466,6 +469,22 @@
           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)
+!
+!--    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,
@@ -502,5 +521,4 @@
 !--       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, 
@@ -536,8 +554,22 @@
           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
+
+!
 !--       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
+          meas_flag = 0
           DO  t = 1, vmea(l)%ntraj
 !
@@ -553,13 +585,5 @@
                      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.
-            
-             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
+
 !
 !--          Compute grid indices relative to origin and check if these are
@@ -567,8 +591,4 @@
 !--          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 )
@@ -579,86 +599,66 @@
                           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
+!--             Check if observation coordinate is on subdomain
                 IF ( on_pe )  THEN
+!
+!--                Determine vertical index which correspond to the observation
+!--                height.
                    ksurf = get_topography_top_index_ji( js, is, 's' )
                    ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
+!
+!--                Set mask array at the observation coordinates. Also, flag the
+!--                surrounding coordinate points, but first check whether the
+!--                surrounding coordinate points are on the subdomain. 
+                   kl = MERGE( ks-1, ks, ks-1 >= nzb  .AND. ks-1 > ksurf )
+                   ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
+                 
+                   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 ) 
                 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
-
+!
+!--       Based on the flag array count the number of observation points.
+          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 ) )
+                ENDDO
+             ENDDO
+          ENDDO
 !
 !--       Store number of observation points on subdomain and allocate index
 !--       arrays.
           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) )
+!
+!--       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 = ns + 1
+                      vmea(l)%i(ns) = is
+                      vmea(l)%j(ns) = js
+                      vmea(l)%k(ns) = ks
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
           
-!           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
-                ENDIF
-!
-!--             Store arrays for next iteration - avoid double counting 
-                i_prev = is
-                j_prev = js
-                k_prev = ks
-             ENDDO
-             
-          ENDDO
+!           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.
@@ -677,5 +677,5 @@
        ENDIF
     ENDDO
-    flush(9)
+!     flush(9)
     
 !
@@ -684,4 +684,8 @@
     CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '',              &
                                 global_attribute, 'close', '' )
+!                               
+!-- Dellocate flag array
+    DEALLOCATE( meas_flag )
+        
   END SUBROUTINE vm_init
   
@@ -694,8 +698,15 @@
   SUBROUTINE vm_sampling
 
-    USE arrays_3d !,                                                             &
-!         ONLY:  pt
-
-    USE surface_mod
+    USE arrays_3d,                                                             &
+        ONLY:  exner, pt, q, u, v, w
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  pi
+    
+    USE radiation_model_mod,                                                   &
+        ONLY:  radiation  
+
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
     
      IMPLICIT NONE
@@ -703,10 +714,12 @@
      CHARACTER(LEN=10) ::  trimvar !< dummy for the measured variable name
       
-     INTEGER(iwp) ::  l !< 
-     INTEGER(iwp) ::  m !<
-     INTEGER(iwp) ::  var !< 
-     
-     INTEGER(iwp) ::  mm, j, i
-
+     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 
@@ -714,44 +727,316 @@
      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
-!                  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
+!--     Loop over all measured variables at this station. Please note, 
+!--     velocity components are interpolated onto scalar grid.  
+        DO  n = 1, vmea(l)%nvar
+        
+           SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
+           
+              CASE ( 'theta' )
+                 IF ( .NOT. neutral )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(n,m) = pt(k,j,i)
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'ta', 't_va' )
+                 IF ( .NOT. neutral )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'hus', 'haa' )
+                 IF ( humidity )  THEN
+                    DO  m = 1, vmea(l)%ns
+                       k = vmea(l)%k(m)
+                       j = vmea(l)%j(m)
+                       i = vmea(l)%i(m)
+                       vmea(l)%measured_vars(n,m) = q(k,j,i)
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'u', 'ua' )
+                 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(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
+                 ENDDO
+                 
+              CASE ( 'v', 'va' )
+                 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(n,m) = 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)
+                    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) )
+                 ENDDO
+                 
+              CASE ( 'wspeed' )
+                 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(n,m) = 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   &
+                                                     )
+                 ENDDO
+                 
+              CASE ( 'wdir' )
+                 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(n,m) = 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                 
+                    DO  nn = 1, nspec                    
+                       IF ( TRIM( vmea(l)%measured_vars_name(m) ) ==           &
+                            TRIM( chem_species(nn)%name ) )  THEN                            
+                          DO  m = 1, vmea(l)%ns             
+                             k = vmea(l)%k(m)
+                             j = vmea(l)%j(m)
+                             i = vmea(l)%i(m)                    
+                             vmea(l)%measured_vars(n,m) =                      &
+                                                   chem_species(nn)%conc(k,j,i)
+                          ENDDO
+                       ENDIF
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'us' )
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not 
+!--                 on ghost points. Hence, limit the indices. 
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                    j = MERGE( j           , nyn, j            > nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                    i = MERGE( i           , nxr, i            > nxr )
+                    
+                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
+                             surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(n,m) = 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)
+                    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)
+                    ENDDO
+                 ENDDO
+                 
+              CASE ( 'ts' )
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not 
+!--                 on ghost points. Hence, limit the indices. 
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                    j = MERGE( j           , nyn, j            > nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                    i = MERGE( i           , nxr, i            > nxr )
+                    
+                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
+                             surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(n,m) = 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)
+                    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)
+                    ENDDO
+                 ENDDO
+                 
+              CASE ( 'hfls' )
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not 
+!--                 on ghost points. Hence, limit the indices. 
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                    j = MERGE( j           , nyn, j            > nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                    i = MERGE( i           , nxr, i            > nxr )
+                    
+                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
+                             surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(n,m) = 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)
+                    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)
+                    ENDDO
+                 ENDDO
+                 
+              CASE ( 'hfss' )
+                 DO  m = 1, vmea(l)%ns
+!
+!--                 Surface data is only available on inner subdomains, not 
+!--                 on ghost points. Hence, limit the indices. 
+                    j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                    j = MERGE( j           , nyn, j            > nyn )
+                    i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                    i = MERGE( i           , nxr, i            > nxr )
+                    
+                    DO  mm = surf_def_h(0)%start_index(j,i),                   &
+                             surf_def_h(0)%end_index(j,i)
+                       vmea(l)%measured_vars(n,m) = 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)
+                    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)
+                    ENDDO
+                 ENDDO
+                 
+              CASE ( 'rnds' )
+                 IF ( radiation )  THEN
+                    DO  m = 1, vmea(l)%ns
+!
+!--                    Surface data is only available on inner subdomains, not 
+!--                    on ghost points. Hence, limit the indices. 
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                       j = MERGE( j           , nyn, j            > nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                       i = MERGE( i           , nxr, i            > nxr )
+                    
+                       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)
+                       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)
+                       ENDDO
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'rsus', 'rsu' )
+                 IF ( radiation )  THEN
+                    DO  m = 1, vmea(l)%ns
+!
+!--                    Surface data is only available on inner subdomains, not 
+!--                    on ghost points. Hence, limit the indices. 
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                       j = MERGE( j           , nyn, j            > nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                       i = MERGE( i           , nxr, i            > nxr )
+                    
+                       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)
+                       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)
+                       ENDDO
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'rsds', 'rsd' )
+                 IF ( radiation )  THEN
+                    DO  m = 1, vmea(l)%ns
+!
+!--                    Surface data is only available on inner subdomains, not 
+!--                    on ghost points. Hence, limit the indices. 
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                       j = MERGE( j           , nyn, j            > nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                       i = MERGE( i           , nxr, i            > nxr )
+                    
+                       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)
+                       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)
+                       ENDDO
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'rlus', 'rlu' )
+                 IF ( radiation )  THEN
+                    DO  m = 1, vmea(l)%ns
+!
+!--                    Surface data is only available on inner subdomains, not 
+!--                    on ghost points. Hence, limit the indices. 
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                       j = MERGE( j           , nyn, j            > nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                       i = MERGE( i           , nxr, i            > nxr )
+                    
+                       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)
+                       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)
+                       ENDDO
+                    ENDDO
+                 ENDIF
+                 
+              CASE ( 'rlds', 'rld' )
+                 IF ( radiation )  THEN
+                    DO  m = 1, vmea(l)%ns
+!
+!--                    Surface data is only available on inner subdomains, not 
+!--                    on ghost points. Hence, limit the indices. 
+                       j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
+                       j = MERGE( j           , nyn, j            > nyn )
+                       i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
+                       i = MERGE( i           , nxr, i            > nxr )
+                    
+                       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)
+                       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)
+                       ENDDO
+                    ENDDO
+                 ENDIF
+!
+!--           More will follow ...
+                 
+           END SELECT
+
         ENDDO
 
