!> @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 .
!
! Copyright 2017 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id: virtual_measurement_mod.f90 3473 2018-10-30 20:50:15Z maronga $
! Initial revision
!
! 3472 2018-10-30 20:43:50Z suehring
!
! 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,
! 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