Home

Context Navigation

source: palm/trunk/SOURCE/virtual_measurement_mod.f90 @ 3704

Last change on this file since 3704 was 3704, checked in by suehring, 5 years ago
Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 73.2 KB

Rev	Line
[3471]	1	!> @virtual_measurement_mod.f90
[3434]	2	!------------------------------------------------------------------------------!
	3	! This file is part of the PALM model system.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 2017 Leibniz Universitaet Hannover
	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
[3704]	22	! - initialization revised
	23	! - binary data output
	24	! - list of allowed variables extended
[3434]	25	!
	26	! Former revisions:
	27	! -----------------
[3472]	28	! $Id: virtual_measurement_mod.f90 3704 2019-01-29 19:51:41Z suehring $
[3522]	29	! Sampling of variables
	30	!
	31	! 3494 2018-11-06 14:51:27Z suehring
[3494]	32	! Bugfixing
	33	!
	34	! 3473 2018-10-30 20:50:15Z suehring
[3473]	35	! Initial revision
[3434]	36	!
	37	! Authors:
	38	! --------
[3522]	39	! @author Matthias Suehring
[3434]	40	!
	41	! Description:
	42	! ------------
[3471]	43	!> The module acts as an interface between 'real-world' observations and
	44	!> model simulations. Virtual measurements will be taken in the model at the
[3704]	45	!> coordinates representative for the 'real-world' observation coordinates.
[3471]	46	!> More precisely, coordinates and measured quanties will be read from a
	47	!> NetCDF file which contains all required information. In the model,
	48	!> the same quantities (as long as all the required components are switched-on)
	49	!> will be sampled at the respective positions and output into an extra file,
	50	!> which allows for straight-forward comparison of model results with
	51	!> observations.
[3522]	52	!>
	53	!> @todo list_of_allowed variables needs careful checking
	54	!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
	55	!> the (UC)2 standard
	56	!> @note Fluxes are not processed
[3434]	57	!------------------------------------------------------------------------------!
[3471]	58	MODULE virtual_measurement_mod
[3434]	59
	60	USE arrays_3d, &
	61	ONLY: q, pt, u, v, w, zu, zw
	62
[3522]	63	USE chem_modules, &
	64	ONLY: nspec
	65
	66	USE chemistry_model_mod, &
	67	ONLY: chem_species
	68
[3434]	69	USE control_parameters, &
[3704]	70	ONLY: air_chemistry, dz, humidity, io_blocks, io_group, neutral, &
	71	message_string, time_since_reference_point, virtual_measurement
[3434]	72
	73	USE cpulog, &
	74	ONLY: cpu_log, log_point
	75
	76	USE grid_variables, &
	77	ONLY: dx, dy
	78
	79	USE indices, &
[3704]	80	ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0
[3434]	81
	82	USE kinds
[3704]	83
	84	USE netcdf_data_input_mod, &
	85	ONLY: init_model
	86
	87	USE pegrid
	88
	89	USE surface_mod, &
	90	ONLY: surf_lsm_h, surf_usm_h
	91
	92	USE land_surface_model_mod, &
	93	ONLY: nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h
	94
	95	USE radiation_model_mod
	96
	97	USE urban_surface_mod, &
	98	ONLY: nzb_wall, nzt_wall, t_wall_h
[3434]	99
	100
	101	IMPLICIT NONE
[3704]	102
	103	TYPE virt_general
	104	INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements
	105	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
	106	END TYPE virt_general
[3434]	107
	108	TYPE virt_mea
	109
[3704]	110	CHARACTER(LEN=100) :: feature_type !< type of the measurement
	111	CHARACTER(LEN=100) :: filename_original !< name of the original file
	112	CHARACTER(LEN=100) :: site !< name of the measurement site
[3434]	113
	114	CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables
	115
[3704]	116	INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements
	117	INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements
	118	INTEGER(iwp) :: ntraj !< number of trajectories of a measurement
	119	INTEGER(iwp) :: nvar !< number of measured variables (atmosphere + soil)
[3434]	120
[3704]	121	INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements
	122	INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
	123
[3434]	124	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues
	125
[3704]	126	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
	127	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
	128	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
	129
	130	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction
	131	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction
	132	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction
[3434]	133
	134	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
	135	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
	136	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
[3704]	137	LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled
[3434]	138
[3704]	139	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
	140	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
	141	REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values
	142	REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building
	143	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
	144	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
	145
	146	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_ag !< measurement height above ground level
	147	REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil
[3522]	148
[3704]	149	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
	150	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables
[3434]	151
	152	END TYPE virt_mea
	153
	154	CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting
	155	CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type
[3704]	156
	157	! This need to generalized
	158	CHARACTER(LEN=8) :: char_filename = "filename" !< attribute name for filename
	159	CHARACTER(LEN=11) :: char_soil = "soil_sample" !< attribute name for soil sampling indication
[3434]	160	CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue
	161	CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names
	162	CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing
	163	CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations
	164	CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x
	165	CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y
	166	CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name
	167	CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable
	168	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
	169	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
	170	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
[3704]	171
	172	CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ & !< list of soil variables
	173	't_soil', &
	174	'm_soil', &
	175	'lwc ', &
	176	'lwcs ', &
	177	'smp ' /)
	178
	179	CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ &
	180	'mcpm1 ', 'kc_PM1 ', &
	181	'mcpm2p5 ', 'kc_PM2.5 ', &
	182	'mcpm25 ', 'kc_PM25 ', &
	183	'mcpm10 ', 'kc_PM10 ', &
	184	'mfno2 ', 'kc_NO2 ', &
	185	'mfno ', 'kc_NO ', &
	186	'tro3 ', 'kc_O3 ', &
	187	'mfco ', 'kc_CO ' &
	188	/), (/ 2, 8 /) )
[3522]	189	!
	190	!-- MS: List requires careful revision!
[3704]	191	CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM
[3471]	192	(/ 'hfls ', & ! surface latent heat flux (W/m2)
	193	'hfss ', & ! surface sensible heat flux (W/m2)
	194	'hur ', & ! relative humidity (-)
	195	'hus ', & ! specific humidity (g/kg)
	196	'haa ', & ! absolute atmospheric humidity (kg/m3)
	197	'mcpm1 ', & ! mass concentration of PM1 (kg/m3)
	198	'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3)
	199	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
	200	'mcco ', & ! mass concentration of CO (kg/m3)
	201	'mcco2 ', & ! mass concentration of CO2 (kg/m3)
	202	'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3)
	203	'ncaa ', & ! number concentation of particles (1/m3)
[3522]	204	'mfco ', & ! mole fraction of CO (mol/mol)
[3471]	205	'mfco2 ', & ! mole fraction of CO2 (mol/mol)
	206	'mfch4 ', & ! mole fraction of methane (mol/mol)
	207	'mfnh3 ', & ! mole fraction of amonia (mol/mol)
	208	'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol)
	209	'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol)
	210	'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol)
	211	'mfh20 ', & ! mole fraction of water (mol/mol)
	212	'plev ', & ! ? air pressure - hydrostaic + perturbation?
	213	'rlds ', & ! surface downward longwave flux (W/m2)
	214	'rlus ', & ! surface upward longwave flux (W/m2)
	215	'rsds ', & ! surface downward shortwave flux (W/m2)
	216	'rsus ', & ! surface upward shortwave flux (W/m2)
	217	'ta ', & ! air temperature (degree C)
	218	't_va ', & ! virtual accoustic temperature (K)
	219	'theta ', & ! potential temperature (K)
	220	'tro3 ', & ! mole fraction of ozone air (mol/mol)
	221	'ts ', & ! scaling parameter of temperature (K)
	222	'wspeed ', & ! ? wind speed - horizontal?
	223	'wdir ', & ! wind direction
	224	'us ', & ! friction velocity
	225	'msoil ', & ! ? soil moisture - which depth?
	226	'tsoil ', & ! ? soil temperature - which depth?
	227	'u ', & ! u-component
[3704]	228	'utheta ', & ! total eastward kinematic heat flux
[3471]	229	'ua ', & ! eastward wind (is there any difference to u?)
	230	'v ', & ! v-component
[3704]	231	'vtheta ', & ! total northward kinematic heat flux
[3471]	232	'va ', & ! northward wind (is there any difference to v?)
	233	'w ', & ! w-component
[3704]	234	'wtheta ', & ! total vertical kinematic heat flux
[3471]	235	'rld ', & ! downward longwave radiative flux (W/m2)
	236	'rlu ', & ! upnward longwave radiative flux (W/m2)
	237	'rsd ', & ! downward shortwave radiative flux (W/m2)
	238	'rsu ', & ! upward shortwave radiative flux (W/m2)
	239	'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2)
[3704]	240	'rnds ', & ! surface net downward radiative flux (W/m2)
	241	't_soil ', &
	242	'm_soil ', &
	243	'lwc ', &
	244	'lwcs ', &
	245	'smp ' &
[3471]	246	/)
[3704]	247
[3434]	248
[3704]	249	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
	250	LOGICAL :: init = .TRUE. !< flag indicating initialization of data output
[3434]	251	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
	252
	253	REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start
	254
[3704]	255	TYPE( virt_general ) :: vmea_general !< data structure which encompass general variables
[3434]	256	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure
	257
	258	INTERFACE vm_check_parameters
	259	MODULE PROCEDURE vm_check_parameters
	260	END INTERFACE vm_check_parameters
	261
[3704]	262	INTERFACE vm_data_output
	263	MODULE PROCEDURE vm_data_output
	264	END INTERFACE vm_data_output
	265
[3434]	266	INTERFACE vm_init
	267	MODULE PROCEDURE vm_init
	268	END INTERFACE vm_init
	269
[3704]	270	INTERFACE vm_last_actions
	271	MODULE PROCEDURE vm_last_actions
	272	END INTERFACE vm_last_actions
	273
[3434]	274	INTERFACE vm_parin
	275	MODULE PROCEDURE vm_parin
	276	END INTERFACE vm_parin
	277
	278	INTERFACE vm_sampling
	279	MODULE PROCEDURE vm_sampling
	280	END INTERFACE vm_sampling
	281
	282	SAVE
	283
	284	PRIVATE
	285
	286	!
	287	!-- Public interfaces
[3704]	288	PUBLIC vm_check_parameters, vm_data_output, vm_init, vm_last_actions, &
	289	vm_parin, vm_sampling
[3434]	290
	291	!
	292	!-- Public variables
[3704]	293	PUBLIC vmea, vmea_general, vm_time_start
[3434]	294
	295	CONTAINS
	296
	297
	298	!------------------------------------------------------------------------------!
	299	! Description:
	300	! ------------
[3471]	301	!> Check parameters for virtual measurement module
[3434]	302	!------------------------------------------------------------------------------!
	303	SUBROUTINE vm_check_parameters
	304
	305	USE control_parameters, &
	306	ONLY: message_string, virtual_measurement
	307
	308	USE netcdf_data_input_mod, &
	309	ONLY: input_pids_static
	310
	311	IMPLICIT NONE
	312
	313	!
	314	!-- In case virtual measurements are taken, a static input file is required.
	315	!-- This is because UTM coordinates for the PALM domain origin are required
	316	!-- for correct mapping of the measurements.
	317	!-- ToDo: Revise this later and remove this requirement.
	318	IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN
[3704]	319	message_string = 'If virtual measurements are taken, a static input ' //&
[3434]	320	'file is mandatory.'
	321	CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
	322	ENDIF
	323
	324	END SUBROUTINE vm_check_parameters
	325
	326	!------------------------------------------------------------------------------!
	327	! Description:
	328	! ------------
[3471]	329	!> Read namelist for the virtual measurement module
[3434]	330	!------------------------------------------------------------------------------!
	331	SUBROUTINE vm_parin
	332
	333	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	334
	335	NAMELIST /virtual_measurement_parameters/ use_virtual_measurement, &
	336	vm_time_start
	337
	338	line = ' '
	339
	340	!
	341	!-- Try to find stg package
	342	REWIND ( 11 )
	343	line = ' '
	344	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
	345	READ ( 11, '(A)', END=20 ) line
	346	ENDDO
	347	BACKSPACE ( 11 )
	348
	349	!
	350	!-- Read namelist
	351	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
	352
	353	!
[3471]	354	!-- Set flag that indicates that the virtual measurement module is switched on
[3434]	355	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
	356
	357	GOTO 20
	358
	359	10 BACKSPACE( 11 )
	360	READ( 11 , '(A)') line
	361	CALL parin_fail_message( 'virtual_measurement_parameters', line )
	362
	363	20 CONTINUE
	364
	365	END SUBROUTINE vm_parin
	366
	367
	368	!------------------------------------------------------------------------------!
	369	! Description:
	370	! ------------
	371	!> Initialize virtual measurements: read coordiante arrays and measured
	372	!> variables, set indicies indicating the measurement points, read further
	373	!> attributes, etc..
	374	!------------------------------------------------------------------------------!
	375	SUBROUTINE vm_init
	376
	377	USE arrays_3d, &
	378	ONLY: zu, zw
	379
	380	USE grid_variables, &
	381	ONLY: ddx, ddy, dx, dy
	382
	383	USE indices, &
	384	ONLY: nxl, nxr, nyn, nys
	385
	386	USE netcdf_data_input_mod, &
[3704]	387	ONLY: init_model, input_file_vm, &
	388	netcdf_data_input_get_dimension_length, &
[3434]	389	netcdf_data_input_att, netcdf_data_input_var
	390
	391	USE surface_mod, &
	392	ONLY: get_topography_top_index_ji
	393
	394	IMPLICIT NONE
	395
	396	CHARACTER(LEN=5) :: dum !< dummy string indicate station id
[3704]	397	CHARACTER(LEN=5) :: dummy_read !< dummy string indicate station id
[3434]	398	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
[3522]	399	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
[3434]	400
[3704]	401	INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate
	402	INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate
[3434]	403	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
[3704]	404	INTEGER(iwp) :: dim_zag !< dimension size of height coordinate
	405	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
	406	INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement
[3434]	407	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
[3704]	408	INTEGER(iwp) :: j !< grid index of observation point in x-direction
[3434]	409	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
[3704]	410	INTEGER(iwp) :: k !< grid index of observation point in x-direction
[3522]	411	INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate
[3434]	412	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
	413	INTEGER(iwp) :: ksurf !< topography top index
[3522]	414	INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate
[3434]	415	INTEGER(iwp) :: l !< running index over all stations
	416	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
	417	INTEGER(iwp) :: ll !< running index over all measured variables in file
	418	INTEGER(iwp) :: lll !< running index over all allowed variables
	419	INTEGER(iwp) :: n !< running index over trajectory coordinates
	420	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
	421	INTEGER(iwp) :: t !< running index over number of trajectories
[3704]	422	INTEGER(iwp) :: m
[3434]	423
[3704]	424	INTEGER(KIND=1):: soil_dum
	425
	426	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates
	427
[3522]	428	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions
	429
	430	LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism
	431	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
	432
[3434]	433	REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM
	434	REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM
	435	REAL(wp) :: fill_zag !< _FillValue for height coordinate
	436
[3437]	437	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable
	438	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable,
	439	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable
[3434]	440	!
[3704]	441	!-- Obtain number of sites. Also, pass the 'open' string, in order to initially
	442	!-- open the measurement driver.
	443	CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, &
	444	vmea_general%id_vm, input_file_vm, &
[3434]	445	global_attribute, 'open', '' )
[3522]	446
[3434]	447	!
[3704]	448	!-- Allocate data structure which encompass all required information, such as
	449	!-- grid points indicies, absolute UTM coordinates, the measured quantities,
	450	!-- etc. .
	451	ALLOCATE( vmea(1:vmea_general%nvm) )
[3434]	452	!
[3704]	453	!-- Allocate flag array. This dummy array is used to identify grid points
	454	!-- where virtual measurements should be taken. Please note, at least one
	455	!-- ghost point is required, in order to include also the surrounding
	456	!-- grid points of the original coordinate.
[3522]	457	ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	458	meas_flag = 0
	459	!
[3704]	460	!-- Loop over all sites.
	461	DO l = 1, vmea_general%nvm
[3434]	462	!
[3704]	463	!-- Determine suffix which contains the ID, ordered according to the number
	464	!-- of measurements.
[3434]	465	IF( l < 10 ) THEN
	466	WRITE( dum, '(I1)') l
	467	ELSEIF( l < 100 ) THEN
	468	WRITE( dum, '(I2)') l
	469	ELSEIF( l < 1000 ) THEN
	470	WRITE( dum, '(I3)') l
	471	ELSEIF( l < 10000 ) THEN
	472	WRITE( dum, '(I4)') l
	473	ELSEIF( l < 100000 ) THEN
	474	WRITE( dum, '(I5)') l
	475	ENDIF
[3704]	476	!
	477	!-- Read site coordinates (UTM).
	478	CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx // &
	479	TRIM( dum ), vmea_general%id_vm, '', &
	480	global_attribute, '', '' )
	481	CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy // &
	482	TRIM( dum ), vmea_general%id_vm, '', &
	483	global_attribute, '', '' )
	484	!
	485	!-- Read site name
	486	CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ), &
	487	vmea_general%id_vm, '', global_attribute, &
[3434]	488	'', '' )
[3704]	489	!
	490	!-- Read type of the measurement (trajectory, profile, timeseries).
	491	CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature // &
	492	TRIM( dum ), vmea_general%id_vm, '', &
	493	global_attribute, '', '' )
	494	!
	495	!-- Read the name of the original file where observational data is stored.
	496	CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
	497	TRIM( dum ), vmea_general%id_vm, '', &
	498	global_attribute, '', '' )
	499	!
	500	!-- Read a flag which indicates that also soil quantities are take at the
	501	!-- respective site (is part of the virtual measurement driver).
	502	CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ), &
	503	vmea_general%id_vm, '', global_attribute, &
[3434]	504	'', '' )
	505	!
[3704]	506	!-- Set flag for soil-sampling.
	507	IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE.
	508	!
[3434]	509	!--- Set logicals depending on the type of the measurement
	510	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
	511	vmea(l)%timseries_profile = .TRUE.
	512	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
	513	vmea(l)%timseries = .TRUE.
	514	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
	515	vmea(l)%trajectory = .TRUE.
[3704]	516	!
	517	!-- Give error message in case the type matches non of the pre-defined types.
[3434]	518	ELSE
	519	message_string = 'Attribue featureType = ' // &
	520	TRIM( vmea(l)%feature_type ) // &
	521	' is not allowed.'
	522	CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
	523	ENDIF
	524	!
[3704]	525	!-- Read string with all measured variables at this site
[3434]	526	measured_variables_file = ''
	527	CALL netcdf_data_input_var( measured_variables_file, &
[3704]	528	char_mv // TRIM( dum ), vmea_general%id_vm )
[3434]	529	!
[3704]	530	!-- Count the number of measured variables. Only count variables that match
	531	!-- with the allowed variables.
	532	!-- Please note, for some NetCDF interal reasons characters end with a NULL,
	533	!-- i.e. also empty characters contain a NULL. Therefore, check the strings
	534	!-- for a NULL to get the correct character length in order to compare
	535	!-- them with the list of allowed variables.
	536	vmea(l)%nvar = 0
[3434]	537	DO ll = 1, SIZE( measured_variables_file )
	538	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
	539	measured_variables_file(ll)(1:1) /= ' ') THEN
	540	!
	541	!-- Obtain character length of the character
	542	len_char = 1
	543	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
	544	.AND. measured_variables_file(ll)(len_char:len_char) /= ' ' )
	545	len_char = len_char + 1
	546	ENDDO
	547	len_char = len_char - 1
	548	!
	549	!-- Now, compare the measured variable with the list of allowed
	550	!-- variables.
	551	DO lll= 1, SIZE( list_allowed_variables )
	552	IF ( measured_variables_file(ll)(1:len_char) == &
	553	TRIM( list_allowed_variables(lll) ) ) THEN
	554	vmea(l)%nvar = vmea(l)%nvar + 1
	555	measured_variables(vmea(l)%nvar) = &
	556	measured_variables_file(ll)(1:len_char)
	557	ENDIF
	558	ENDDO
	559	ENDIF
	560	ENDDO
	561	!
[3704]	562	!-- Allocate array for the measured variables names for the respective site.
[3434]	563	ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
	564
	565	DO ll = 1, vmea(l)%nvar
	566	vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
	567	ENDDO
	568	!
[3522]	569	!-- In case of chemistry, check if species is considered in the modelled
	570	!-- chemistry mechanism.
[3704]	571	! IF ( air_chemistry ) THEN
	572	! DO ll = 1, vmea(l)%nvar
	573	! chem_include = .FALSE.
	574	! DO n = 1, nspec
	575	! IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == &
	576	! TRIM( chem_species(n)%name ) ) chem_include = .TRUE.
	577	! ENDDO
	578	! !
	579	! !-- Revise this. It should only check for chemistry variables and not for all!
	580	! IF ( .NOT. chem_include ) THEN
	581	! message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // &
	582	! ' is not considered in the modelled ' // &
	583	! 'chemistry mechanism'
	584	! CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
	585	! ENDIF
	586	! ENDDO
	587	! ENDIF
[3522]	588	!
[3704]	589	!-- Read the UTM coordinates for the actual site. Based on the coordinates,
	590	!-- define the grid-index space on each subdomain where virtual measurements
	591	!-- should be taken. Note, the entire coordinate arrays will not be stored
	592	!-- as this would exceed memory requirements, particularly for trajectory
	593	!-- measurements.
[3434]	594	IF ( vmea(l)%nvar > 0 ) THEN
	595	!
	596	!-- For stationary measurements UTM coordinates are just one value and
	597	!-- its dimension is "station", while for mobile measurements UTM
[3704]	598	!-- coordinates are arrays depending on the number of trajectories and
	599	!-- time, according to (UC)2 standard. First, inquire dimension length
	600	!-- of the UTM coordinates.
[3434]	601	IF ( vmea(l)%trajectory ) THEN
	602	!
	603	!-- For non-stationary measurements read the number of trajectories
[3704]	604	!-- and the number of time coordinates.
	605	CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
[3434]	606	vmea(l)%ntraj, &
	607	"traj" // &
	608	TRIM( dum ) )
[3704]	609	CALL netcdf_data_input_get_dimension_length( vmea_general%id_vm, &
	610	dim_ntime, &
[3434]	611	"ntime" // &
	612	TRIM( dum ) )
	613	!
[3704]	614	!-- For stationary measurements the dimension for UTM and time
	615	!-- coordinates is 1.
[3434]	616	ELSE
	617	vmea(l)%ntraj = 1
	618	dim_ntime = 1
	619	ENDIF
	620	!
	621	!- Allocate array which defines individual time frame for each
[3704]	622	!-- trajectory or station.
[3434]	623	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
	624	!
	625	!-- Allocate temporary arrays for UTM and height coordinates. Note,
	626	!-- on file UTM coordinates might be 1D or 2D variables
[3437]	627	ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	628	ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	629	ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) )
[3434]	630	!
[3704]	631	!-- Read _FillValue attributes of the coordinate dimensions.
[3434]	632	CALL netcdf_data_input_att( fill_eutm, char_fillvalue, &
[3704]	633	vmea_general%id_vm, '', &
	634	.NOT. global_attribute, '', &
[3434]	635	char_eutm // TRIM( dum ) )
	636	CALL netcdf_data_input_att( fill_nutm, char_fillvalue, &
[3704]	637	vmea_general%id_vm, '', &
	638	.NOT. global_attribute, '', &
[3434]	639	char_nutm // TRIM( dum ) )
	640	CALL netcdf_data_input_att( fill_zag, char_fillvalue, &
[3704]	641	vmea_general%id_vm, '', &
	642	.NOT. global_attribute, '', &
[3434]	643	char_zag // TRIM( dum ) )
	644	!
	645	!-- Read UTM and height coordinates coordinates for all trajectories and
	646	!-- times.
[3437]	647	IF ( vmea(l)%trajectory ) THEN
[3704]	648	CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), &
	649	vmea_general%id_vm, &
[3437]	650	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
[3704]	651	CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), &
	652	vmea_general%id_vm, &
[3437]	653	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
[3704]	654	CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), &
	655	vmea_general%id_vm, &
[3437]	656	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	657	ELSE
[3704]	658	CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
	659	vmea_general%id_vm )
	660	CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
	661	vmea_general%id_vm )
	662	CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), &
	663	vmea_general%id_vm )
	664	ENDIF
[3434]	665	!
	666	!-- Based on UTM coordinates, check if the measurement station or parts
	667	!-- of the trajectory is on subdomain. This case, setup grid index space
	668	!-- sample these quantities.
[3522]	669	meas_flag = 0
[3434]	670	DO t = 1, vmea(l)%ntraj
[3704]	671	!
	672	!-- First, compute relative x- and y-coordinates with respect to the
	673	!-- lower-left origin of the model domain, which is the difference
	674	!-- betwen UTM coordinates. Note, if the origin is not correct, the
	675	!-- virtual sites will be misplaced.
	676	e_utm(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
	677	n_utm(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
[3434]	678	!
	679	!-- Determine the individual time coordinate length for each station and
	680	!-- trajectory. This is required as several stations and trajectories
	681	!-- are merged into one file but they do not have the same number of
	682	!-- points in time, hence, missing values may occur and cannot be
[3704]	683	!-- processed further. This is actually a work-around for the specific
	684	!-- (UC)2 dataset, but it won't harm in anyway.
[3434]	685	vmea(l)%dim_t(t) = 0
	686	DO n = 1, dim_ntime
[3437]	687	IF ( e_utm(t,n) /= fill_eutm .AND. &
	688	n_utm(t,n) /= fill_nutm .AND. &
	689	z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n
[3434]	690	ENDDO
	691	!
	692	!-- Compute grid indices relative to origin and check if these are
	693	!-- on the subdomain. Note, virtual measurements will be taken also
	694	!-- at grid points surrounding the station, hence, check also for
	695	!-- these grid points.
[3437]	696	DO n = 1, vmea(l)%dim_t(t)
	697	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
	698	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
[3434]	699	!
	700	!-- Is the observation point on subdomain?
	701	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
	702	js >= nys .AND. js <= nyn )
	703	!
[3522]	704	!-- Check if observation coordinate is on subdomain
[3434]	705	IF ( on_pe ) THEN
[3522]	706	!
	707	!-- Determine vertical index which correspond to the observation
	708	!-- height.
[3434]	709	ksurf = get_topography_top_index_ji( js, is, 's' )
[3437]	710	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
[3434]	711	!
[3522]	712	!-- Set mask array at the observation coordinates. Also, flag the
	713	!-- surrounding coordinate points, but first check whether the
	714	!-- surrounding coordinate points are on the subdomain.
[3704]	715	kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 >= ksurf )
	716	ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
[3522]	717
[3704]	718	DO i = is-1, is+1
	719	DO j = js-1, js+1
	720	DO k = kl, ku
	721	meas_flag(k,j,i) = MERGE( &
	722	IBSET( meas_flag(k,j,i), 0 ), &
	723	0, &
	724	BTEST( wall_flags_0(k,j,i), 0 ) &
	725	)
	726	ENDDO
	727	ENDDO
	728	ENDDO
[3434]	729	ENDIF
	730	ENDDO
	731
	732	ENDDO
	733	!
[3704]	734	!-- Based on the flag array count the number of of sampling coordinates.
	735	!-- Please note, sampling coordinates in atmosphere and soil may be
	736	!-- different, as within the soil all levels will be measured.
	737	!-- Hence, count individually. Start with atmoshere.
[3522]	738	ns = 0
[3704]	739	DO i = nxl-1, nxr+1
	740	DO j = nys-1, nyn+1
	741	DO k = nzb, nzt+1
	742	ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
[3522]	743	ENDDO
	744	ENDDO
	745	ENDDO
[3704]	746
[3522]	747	!
[3434]	748	!-- Store number of observation points on subdomain and allocate index
[3704]	749	!-- arrays as well as array containing height information.
[3434]	750	vmea(l)%ns = ns
	751
	752	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
	753	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
	754	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
[3704]	755	ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )
[3434]	756	!
[3522]	757	!-- Based on the flag array store the grid indices which correspond to
	758	!-- the observation coordinates.
[3704]	759	ns = 0
	760	DO i = nxl-1, nxr+1
	761	DO j = nys-1, nyn+1
	762	DO k = nzb, nzt+1
	763	IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN
[3522]	764	ns = ns + 1
[3704]	765	vmea(l)%i(ns) = i
	766	vmea(l)%j(ns) = j
	767	vmea(l)%k(ns) = k
	768	vmea(l)%z_ag(ns) = zu(k) - &
	769	zw(get_topography_top_index_ji( j, i, 's' ))
[3522]	770	ENDIF
	771	ENDDO
[3434]	772	ENDDO
	773	ENDDO
	774	!
[3704]	775	!-- Same for the soil. Based on the flag array, count the number of
	776	!-- sampling coordinates in soil. Sample at all soil levels in this case.
	777	IF ( vmea(l)%soil_sampling ) THEN
	778	DO i = nxl, nxr
	779	DO j = nys, nyn
	780	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
	781	IF ( surf_lsm_h%start_index(j,i) <= &
	782	surf_lsm_h%end_index(j,i) ) THEN
	783	vmea(l)%ns_soil = vmea(l)%ns_soil + &
	784	nzt_soil - nzb_soil + 1
	785	ENDIF
	786	IF ( surf_usm_h%start_index(j,i) <= &
	787	surf_usm_h%end_index(j,i) ) THEN
	788	vmea(l)%ns_soil = vmea(l)%ns_soil + &
	789	nzt_wall - nzb_wall + 1
	790	ENDIF
	791	ENDIF
	792	ENDDO
	793	ENDDO
	794	ENDIF
	795	!
	796	!-- Allocate index arrays as well as array containing height information
	797	!-- for soil.
	798	IF ( vmea(l)%soil_sampling ) THEN
	799	ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
	800	ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
	801	ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
	802	ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
	803	ENDIF
	804	!
	805	!-- For soil, store the grid indices.
	806	ns = 0
	807	IF ( vmea(l)%soil_sampling ) THEN
	808	DO i = nxl, nxr
	809	DO j = nys, nyn
	810	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
	811	IF ( surf_lsm_h%start_index(j,i) <= &
	812	surf_lsm_h%end_index(j,i) ) THEN
	813	m = surf_lsm_h%start_index(j,i)
	814	DO k = nzb_soil, nzt_soil
	815	ns = ns + 1
	816	vmea(l)%i_soil(ns) = i
	817	vmea(l)%j_soil(ns) = j
	818	vmea(l)%k_soil(ns) = k
	819	vmea(l)%depth(ns) = zs(k)
	820	ENDDO
	821	ENDIF
	822
	823	IF ( surf_usm_h%start_index(j,i) <= &
	824	surf_usm_h%end_index(j,i) ) THEN
	825	m = surf_usm_h%start_index(j,i)
	826	DO k = nzb_wall, nzt_wall
	827	ns = ns + 1
	828	vmea(l)%i_soil(ns) = i
	829	vmea(l)%j_soil(ns) = j
	830	vmea(l)%k_soil(ns) = k
	831	vmea(l)%depth(ns) = surf_usm_h%zw(k,m)
	832	ENDDO
	833	ENDIF
	834	ENDIF
	835	ENDDO
	836	ENDDO
	837	ENDIF
	838	!
[3434]	839	!-- Allocate array to save the sampled values.
[3704]	840	ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) )
	841
	842	IF ( vmea(l)%soil_sampling ) &
	843	ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, &
	844	1:vmea(l)%nvar) )
[3434]	845	!
[3704]	846	!-- Initialize with _FillValues
	847	vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nvar) = vmea(l)%fillout
	848	IF ( vmea(l)%soil_sampling ) &
	849	vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nvar) = &
	850	vmea(l)%fillout
[3434]	851	!
	852	!-- Deallocate temporary coordinate arrays
	853	IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm )
	854	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
	855	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag )
[3704]	856	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( vmea(l)%dim_t )
[3434]	857	ENDIF
	858	ENDDO
	859	!
	860	!-- Close input file for virtual measurements. Therefore, just call
	861	!-- the read attribute routine with the "close" option.
[3704]	862	CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, &
	863	vmea_general%id_vm, '', &
[3434]	864	global_attribute, 'close', '' )
[3704]	865	!
	866	!-- Sum-up the number of observation coordiates, for atmosphere first.
	867	!-- This is actually only required for data output.
	868	ALLOCATE( ns_all(1:vmea_general%nvm) )
	869	ns_all = 0
	870	#if defined( __parallel )
	871	CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, &
	872	MPI_SUM, comm2d, ierr )
	873	#else
	874	ns_all(:) = vmea(:)%ns
	875	#endif
	876	vmea(:)%ns_tot = ns_all(:)
	877	!
	878	!-- Now for soil
	879	ns_all = 0
	880	#if defined( __parallel )
	881	CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, &
	882	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	883	#else
	884	ns_all(:) = vmea(:)%ns_soil
	885	#endif
	886	vmea(:)%ns_soil_tot = ns_all(:)
	887
	888	DEALLOCATE( ns_all )
[3522]	889	!
	890	!-- Dellocate flag array
	891	DEALLOCATE( meas_flag )
[3704]	892	!
	893	!-- Initialize binary data output of virtual measurements.
	894	!-- Open binary output file.
	895	CALL check_open( 27 )
	896	!
	897	!-- Output header information.
	898	CALL vm_data_output
[3522]	899
[3434]	900	END SUBROUTINE vm_init
	901
	902
	903	!------------------------------------------------------------------------------!
	904	! Description:
	905	! ------------
[3704]	906	!> Binary data output.
	907	!------------------------------------------------------------------------------!
	908	SUBROUTINE vm_data_output
	909
	910	USE pegrid
	911
	912	IMPLICIT NONE
	913
	914	INTEGER(iwp) :: i !< running index over IO blocks
	915	INTEGER(iwp) :: l !< running index over all stations
	916	INTEGER(iwp) :: n !< running index over all measured variables at a station
	917	!
	918	!-- Header output on each PE
	919	IF ( init ) THEN
	920
	921	DO i = 0, io_blocks-1
	922	IF ( i == io_group ) THEN
	923	WRITE ( 27 ) 'number of measurements '
	924	WRITE ( 27 ) vmea_general%nvm
	925
	926	DO l = 1, vmea_general%nvm
	927	WRITE ( 27 ) 'site '
	928	WRITE ( 27 ) vmea(l)%site
	929	WRITE ( 27 ) 'file '
	930	WRITE ( 27 ) vmea(l)%filename_original
	931	WRITE ( 27 ) 'feature_type '
	932	WRITE ( 27 ) vmea(l)%feature_type
	933	WRITE ( 27 ) 'origin_x_obs '
	934	WRITE ( 27 ) vmea(l)%origin_x_obs
	935	WRITE ( 27 ) 'origin_y_obs '
	936	WRITE ( 27 ) vmea(l)%origin_y_obs
	937	WRITE ( 27 ) 'total number of observation points'
	938	WRITE ( 27 ) vmea(l)%ns_tot
	939	WRITE ( 27 ) 'number of measured variables '
	940	WRITE ( 27 ) vmea(l)%nvar
	941	WRITE ( 27 ) 'variables '
	942	WRITE ( 27 ) vmea(l)%measured_vars_name(:)
	943	WRITE ( 27 ) 'number of observation points '
	944	WRITE ( 27 ) vmea(l)%ns
	945	WRITE ( 27 ) 'E_UTM '
	946	WRITE ( 27 ) init_model%origin_x + &
	947	REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx
	948	WRITE ( 27 ) 'N_UTM '
	949	WRITE ( 27 ) init_model%origin_y + &
	950	REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy
	951	WRITE ( 27 ) 'Z_AG '
	952	WRITE ( 27 ) vmea(l)%z_ag(1:vmea(l)%ns)
	953	WRITE ( 27 ) 'soil sampling '
	954	WRITE ( 27 ) MERGE( 'yes ', &
	955	'no ', &
	956	vmea(l)%soil_sampling )
	957
	958	IF ( vmea(l)%soil_sampling ) THEN
	959	WRITE ( 27 ) 'total number of soil points '
	960	WRITE ( 27 ) vmea(l)%ns_soil_tot
	961	print*, "vmea(l)%ns_soil_tot", vmea(l)%ns_soil_tot
	962	WRITE ( 27 ) 'number of soil points '
	963	WRITE ( 27 ) vmea(l)%ns_soil
	964	WRITE ( 27 ) 'E_UTM soil '
	965	WRITE ( 27 ) init_model%origin_x + &
	966	REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
	967	KIND = wp ) * dx
	968	WRITE ( 27 ) 'N_UTM soil '
	969	WRITE ( 27 ) init_model%origin_y + &
	970	REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
	971	KIND = wp ) * dy
	972	WRITE ( 27 ) 'DEPTH '
	973	WRITE ( 27 ) vmea(l)%depth(1:vmea(l)%ns_soil)
	974	ENDIF
	975	ENDDO
	976
	977	ENDIF
	978	ENDDO
	979
	980	#if defined( __parallel )
	981	CALL MPI_BARRIER( comm2d, ierr )
	982	#endif
	983	!
	984	!-- After header information is written, set control flag to .FALSE.
	985	init = .FALSE.
	986	!
	987	!-- Data output at each measurement timestep on each PE
	988	ELSE
	989	DO i = 0, io_blocks-1
	990
	991	IF ( i == io_group ) THEN
	992	WRITE( 27 ) 'output time '
	993	WRITE( 27 ) time_since_reference_point
	994	DO l = 1, vmea_general%nvm
	995	!
	996	!-- Skip binary writing if no observation points are defined on PE
	997	IF ( vmea(l)%ns < 1 .AND. vmea(l)%ns_soil < 1) CYCLE
	998	DO n = 1, vmea(l)%nvar
	999	WRITE( 27 ) vmea(l)%measured_vars_name(n)
	1000	IF ( vmea(l)%soil_sampling .AND. &
	1001	ANY( TRIM( vmea(l)%measured_vars_name(n)) == &
	1002	soil_vars ) ) THEN
	1003	WRITE( 27 ) vmea(l)%measured_vars_soil(:,n)
	1004	ELSE
	1005	WRITE( 27 ) vmea(l)%measured_vars(:,n)
	1006	ENDIF
	1007	ENDDO
	1008
	1009	ENDDO
	1010	ENDIF
	1011	ENDDO
	1012	#if defined( __parallel )
	1013	CALL MPI_BARRIER( comm2d, ierr )
	1014	#endif
	1015	ENDIF
	1016
	1017	END SUBROUTINE vm_data_output
	1018
	1019
	1020	!------------------------------------------------------------------------------!
	1021	! Description:
	1022	! ------------
	1023	!> Write end-of-file statement as last action.
	1024	!------------------------------------------------------------------------------!
	1025	SUBROUTINE vm_last_actions
	1026
	1027	USE pegrid
	1028
	1029	IMPLICIT NONE
	1030
	1031	INTEGER(iwp) :: i !< running index over IO blocks
	1032	INTEGER(iwp) :: l !< running index over all stations
	1033	INTEGER(iwp) :: n !< running index over all measured variables at a station
	1034
	1035	DO i = 0, io_blocks-1
	1036	IF ( i == io_group ) THEN
	1037	WRITE( 27 ) 'EOF '
	1038	ENDIF
	1039	ENDDO
	1040	#if defined( __parallel )
	1041	CALL MPI_BARRIER( comm2d, ierr )
	1042	#endif
	1043	!
	1044	!-- Close binary file
	1045	CALL close_file( 27 )
	1046
	1047	END SUBROUTINE vm_last_actions
	1048
	1049	!------------------------------------------------------------------------------!
	1050	! Description:
	1051	! ------------
[3434]	1052	!> Sampling of the actual quantities along the observation coordinates
	1053	!------------------------------------------------------------------------------!
[3471]	1054	SUBROUTINE vm_sampling
[3434]	1055
[3522]	1056	USE arrays_3d, &
	1057	ONLY: exner, pt, q, u, v, w
[3471]	1058
[3522]	1059	USE basic_constants_and_equations_mod, &
	1060	ONLY: pi
[3434]	1061
[3522]	1062	USE radiation_model_mod, &
	1063	ONLY: radiation
	1064
	1065	USE surface_mod, &
	1066	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	1067
[3434]	1068	IMPLICIT NONE
	1069
[3704]	1070	CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name
[3434]	1071
[3704]	1072	INTEGER(iwp) :: i !< grid index in x-direction
	1073	INTEGER(iwp) :: j !< grid index in y-direction
	1074	INTEGER(iwp) :: k !< grid index in z-direction
	1075	INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
	1076	INTEGER(iwp) :: l !< running index over the number of stations
	1077	INTEGER(iwp) :: m !< running index over all virtual observation coordinates
	1078	INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate
	1079	INTEGER(iwp) :: n !< running index over all measured variables at a station
	1080	INTEGER(iwp) :: nn !< running index over the number of chemcal species
	1081
	1082	LOGICAL :: match_lsm !< flag indicating natural-type surface
	1083	LOGICAL :: match_usm !< flag indicating urban-type surface
[3434]	1084	!
[3704]	1085	!-- Loop over all sites.
	1086	DO l = 1, vmea_general%nvm
[3434]	1087	!
[3704]	1088	!-- At the beginning, set _FillValues
	1089	IF ( ALLOCATED( vmea(l)%measured_vars ) ) &
	1090	vmea(l)%measured_vars = vmea(l)%fillout
	1091	IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) &
	1092	vmea(l)%measured_vars_soil = vmea(l)%fillout
	1093	!
	1094	!-- Loop over all variables measured at this site.
[3522]	1095	DO n = 1, vmea(l)%nvar
	1096
	1097	SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
	1098
	1099	CASE ( 'theta' )
	1100	IF ( .NOT. neutral ) THEN
	1101	DO m = 1, vmea(l)%ns
	1102	k = vmea(l)%k(m)
	1103	j = vmea(l)%j(m)
	1104	i = vmea(l)%i(m)
[3704]	1105	vmea(l)%measured_vars(m,n) = pt(k,j,i)
[3522]	1106	ENDDO
	1107	ENDIF
	1108
[3704]	1109	CASE ( 'ta' )
[3522]	1110	IF ( .NOT. neutral ) THEN
	1111	DO m = 1, vmea(l)%ns
	1112	k = vmea(l)%k(m)
	1113	j = vmea(l)%j(m)
	1114	i = vmea(l)%i(m)
[3704]	1115	vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
[3522]	1116	ENDDO
	1117	ENDIF
[3704]	1118
	1119	CASE ( 't_va' )
[3522]	1120
	1121	CASE ( 'hus', 'haa' )
	1122	IF ( humidity ) THEN
	1123	DO m = 1, vmea(l)%ns
	1124	k = vmea(l)%k(m)
	1125	j = vmea(l)%j(m)
	1126	i = vmea(l)%i(m)
[3704]	1127	vmea(l)%measured_vars(m,n) = q(k,j,i)
[3522]	1128	ENDDO
	1129	ENDIF
	1130
	1131	CASE ( 'u', 'ua' )
	1132	DO m = 1, vmea(l)%ns
	1133	k = vmea(l)%k(m)
	1134	j = vmea(l)%j(m)
	1135	i = vmea(l)%i(m)
[3704]	1136	vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
[3522]	1137	ENDDO
	1138
	1139	CASE ( 'v', 'va' )
	1140	DO m = 1, vmea(l)%ns
	1141	k = vmea(l)%k(m)
	1142	j = vmea(l)%j(m)
	1143	i = vmea(l)%i(m)
[3704]	1144	vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
[3522]	1145	ENDDO
	1146
	1147	CASE ( 'w' )
	1148	DO m = 1, vmea(l)%ns
[3704]	1149	k = MAX ( 1, vmea(l)%k(m) )
[3522]	1150	j = vmea(l)%j(m)
	1151	i = vmea(l)%i(m)
[3704]	1152	vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[3522]	1153	ENDDO
	1154
	1155	CASE ( 'wspeed' )
	1156	DO m = 1, vmea(l)%ns
	1157	k = vmea(l)%k(m)
	1158	j = vmea(l)%j(m)
	1159	i = vmea(l)%i(m)
[3704]	1160	vmea(l)%measured_vars(m,n) = SQRT( &
[3522]	1161	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
	1162	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
	1163	)
	1164	ENDDO
	1165
	1166	CASE ( 'wdir' )
	1167	DO m = 1, vmea(l)%ns
	1168	k = vmea(l)%k(m)
	1169	j = vmea(l)%j(m)
	1170	i = vmea(l)%i(m)
	1171
[3704]	1172	vmea(l)%measured_vars(m,n) = ATAN2( &
[3522]	1173	- 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), &
	1174	- 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) &
	1175	) * 180.0_wp / pi
	1176	ENDDO
[3704]	1177
	1178	CASE ( 'utheta' )
	1179	DO m = 1, vmea(l)%ns
	1180	k = vmea(l)%k(m)
	1181	j = vmea(l)%j(m)
	1182	i = vmea(l)%i(m)
	1183	vmea(l)%measured_vars(m,n) = 0.5_wp * &
	1184	( u(k,j,i) + u(k,j,i+1) ) * &
	1185	pt(k,j,i)
	1186	ENDDO
	1187
	1188	CASE ( 'vtheta' )
	1189	DO m = 1, vmea(l)%ns
	1190	k = vmea(l)%k(m)
	1191	j = vmea(l)%j(m)
	1192	i = vmea(l)%i(m)
	1193	vmea(l)%measured_vars(m,n) = 0.5_wp * &
	1194	( v(k,j,i) + v(k,j+1,i) ) * &
	1195	pt(k,j,i)
	1196	ENDDO
	1197
	1198	CASE ( 'wtheta' )
	1199	DO m = 1, vmea(l)%ns
	1200	k = MAX ( 1, vmea(l)%k(m) )
	1201	j = vmea(l)%j(m)
	1202	i = vmea(l)%i(m)
	1203	vmea(l)%measured_vars(m,n) = 0.5_wp * &
	1204	( w(k-1,j,i) + w(k,j,i) ) * &
	1205	pt(k,j,i)
	1206	ENDDO
	1207
	1208	CASE ( 'uw' )
	1209	DO m = 1, vmea(l)%ns
	1210	k = MAX ( 1, vmea(l)%k(m) )
	1211	j = vmea(l)%j(m)
	1212	i = vmea(l)%i(m)
	1213	vmea(l)%measured_vars(m,n) = 0.25_wp * &
	1214	( w(k-1,j,i) + w(k,j,i) ) * &
	1215	( u(k,j,i) + u(k,j,i+1) )
	1216	ENDDO
	1217
	1218	CASE ( 'vw' )
	1219	DO m = 1, vmea(l)%ns
	1220	k = MAX ( 1, vmea(l)%k(m) )
	1221	j = vmea(l)%j(m)
	1222	i = vmea(l)%i(m)
	1223	vmea(l)%measured_vars(m,n) = 0.25_wp * &
	1224	( w(k-1,j,i) + w(k,j,i) ) * &
	1225	( v(k,j,i) + v(k,j+1,i) )
	1226	ENDDO
	1227
	1228	CASE ( 'uv' )
	1229	DO m = 1, vmea(l)%ns
	1230	k = MAX ( 1, vmea(l)%k(m) )
	1231	j = vmea(l)%j(m)
	1232	i = vmea(l)%i(m)
	1233	vmea(l)%measured_vars(m,n) = 0.25_wp * &
	1234	( u(k,j,i) + u(k,j,i+1) ) * &
	1235	( v(k,j,i) + v(k,j+1,i) )
	1236	ENDDO
[3522]	1237	!
[3704]	1238	!-- List of variables may need extension.
	1239	CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2', &
	1240	'tro3' )
	1241	IF ( air_chemistry ) THEN
	1242	!
	1243	!-- First, search for the measured variable in the chem_vars
	1244	!-- list, in order to get the internal name of the variable.
	1245	DO nn = 1, UBOUND( chem_vars, 2 )
	1246	IF ( TRIM( vmea(l)%measured_vars_name(m) ) == &
	1247	TRIM( chem_vars(0,nn) ) ) ind_chem = nn
	1248	ENDDO
	1249	!
	1250	!-- Run loop over all chemical species, if the measured
	1251	!-- variable matches the interal name, sample the variable.
[3522]	1252	DO nn = 1, nspec
[3704]	1253	IF ( TRIM( chem_vars(1,ind_chem) ) == &
[3522]	1254	TRIM( chem_species(nn)%name ) ) THEN
	1255	DO m = 1, vmea(l)%ns
	1256	k = vmea(l)%k(m)
	1257	j = vmea(l)%j(m)
	1258	i = vmea(l)%i(m)
[3704]	1259	vmea(l)%measured_vars(m,n) = &
[3522]	1260	chem_species(nn)%conc(k,j,i)
	1261	ENDDO
	1262	ENDIF
	1263	ENDDO
	1264	ENDIF
	1265
	1266	CASE ( 'us' )
	1267	DO m = 1, vmea(l)%ns
	1268	!
	1269	!-- Surface data is only available on inner subdomains, not
	1270	!-- on ghost points. Hence, limit the indices.
	1271	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1272	j = MERGE( j , nyn, j > nyn )
	1273	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1274	i = MERGE( i , nxr, i > nxr )
	1275
	1276	DO mm = surf_def_h(0)%start_index(j,i), &
	1277	surf_def_h(0)%end_index(j,i)
[3704]	1278	vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
[3522]	1279	ENDDO
	1280	DO mm = surf_lsm_h%start_index(j,i), &
	1281	surf_lsm_h%end_index(j,i)
[3704]	1282	vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
[3522]	1283	ENDDO
	1284	DO mm = surf_usm_h%start_index(j,i), &
	1285	surf_usm_h%end_index(j,i)
[3704]	1286	vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
[3522]	1287	ENDDO
	1288	ENDDO
	1289
	1290	CASE ( 'ts' )
	1291	DO m = 1, vmea(l)%ns
	1292	!
	1293	!-- Surface data is only available on inner subdomains, not
	1294	!-- on ghost points. Hence, limit the indices.
	1295	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1296	j = MERGE( j , nyn, j > nyn )
	1297	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1298	i = MERGE( i , nxr, i > nxr )
	1299
	1300	DO mm = surf_def_h(0)%start_index(j,i), &
	1301	surf_def_h(0)%end_index(j,i)
[3704]	1302	vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
[3522]	1303	ENDDO
	1304	DO mm = surf_lsm_h%start_index(j,i), &
	1305	surf_lsm_h%end_index(j,i)
[3704]	1306	vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
[3522]	1307	ENDDO
	1308	DO mm = surf_usm_h%start_index(j,i), &
	1309	surf_usm_h%end_index(j,i)
[3704]	1310	vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
[3522]	1311	ENDDO
	1312	ENDDO
	1313
	1314	CASE ( 'hfls' )
	1315	DO m = 1, vmea(l)%ns
	1316	!
	1317	!-- Surface data is only available on inner subdomains, not
	1318	!-- on ghost points. Hence, limit the indices.
	1319	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1320	j = MERGE( j , nyn, j > nyn )
	1321	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1322	i = MERGE( i , nxr, i > nxr )
	1323
	1324	DO mm = surf_def_h(0)%start_index(j,i), &
	1325	surf_def_h(0)%end_index(j,i)
[3704]	1326	vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
[3522]	1327	ENDDO
	1328	DO mm = surf_lsm_h%start_index(j,i), &
	1329	surf_lsm_h%end_index(j,i)
[3704]	1330	vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
[3522]	1331	ENDDO
	1332	DO mm = surf_usm_h%start_index(j,i), &
	1333	surf_usm_h%end_index(j,i)
[3704]	1334	vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
[3522]	1335	ENDDO
	1336	ENDDO
	1337
	1338	CASE ( 'hfss' )
	1339	DO m = 1, vmea(l)%ns
	1340	!
	1341	!-- Surface data is only available on inner subdomains, not
	1342	!-- on ghost points. Hence, limit the indices.
	1343	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1344	j = MERGE( j , nyn, j > nyn )
	1345	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1346	i = MERGE( i , nxr, i > nxr )
	1347
	1348	DO mm = surf_def_h(0)%start_index(j,i), &
	1349	surf_def_h(0)%end_index(j,i)
[3704]	1350	vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
[3522]	1351	ENDDO
	1352	DO mm = surf_lsm_h%start_index(j,i), &
	1353	surf_lsm_h%end_index(j,i)
[3704]	1354	vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
[3522]	1355	ENDDO
	1356	DO mm = surf_usm_h%start_index(j,i), &
	1357	surf_usm_h%end_index(j,i)
[3704]	1358	vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
[3522]	1359	ENDDO
	1360	ENDDO
	1361
	1362	CASE ( 'rnds' )
	1363	IF ( radiation ) THEN
	1364	DO m = 1, vmea(l)%ns
	1365	!
	1366	!-- Surface data is only available on inner subdomains, not
	1367	!-- on ghost points. Hence, limit the indices.
	1368	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1369	j = MERGE( j , nyn, j > nyn )
	1370	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1371	i = MERGE( i , nxr, i > nxr )
	1372
	1373	DO mm = surf_lsm_h%start_index(j,i), &
	1374	surf_lsm_h%end_index(j,i)
[3704]	1375	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
[3522]	1376	ENDDO
	1377	DO mm = surf_usm_h%start_index(j,i), &
	1378	surf_usm_h%end_index(j,i)
[3704]	1379	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
[3522]	1380	ENDDO
	1381	ENDDO
	1382	ENDIF
	1383
[3704]	1384	CASE ( 'rsus' )
[3522]	1385	IF ( radiation ) THEN
	1386	DO m = 1, vmea(l)%ns
	1387	!
	1388	!-- Surface data is only available on inner subdomains, not
	1389	!-- on ghost points. Hence, limit the indices.
	1390	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1391	j = MERGE( j , nyn, j > nyn )
	1392	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1393	i = MERGE( i , nxr, i > nxr )
	1394
	1395	DO mm = surf_lsm_h%start_index(j,i), &
	1396	surf_lsm_h%end_index(j,i)
[3704]	1397	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
[3522]	1398	ENDDO
	1399	DO mm = surf_usm_h%start_index(j,i), &
	1400	surf_usm_h%end_index(j,i)
[3704]	1401	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
[3522]	1402	ENDDO
	1403	ENDDO
	1404	ENDIF
	1405
[3704]	1406	CASE ( 'rsds' )
[3522]	1407	IF ( radiation ) THEN
	1408	DO m = 1, vmea(l)%ns
	1409	!
	1410	!-- Surface data is only available on inner subdomains, not
	1411	!-- on ghost points. Hence, limit the indices.
	1412	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1413	j = MERGE( j , nyn, j > nyn )
	1414	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1415	i = MERGE( i , nxr, i > nxr )
	1416
	1417	DO mm = surf_lsm_h%start_index(j,i), &
	1418	surf_lsm_h%end_index(j,i)
[3704]	1419	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
[3522]	1420	ENDDO
	1421	DO mm = surf_usm_h%start_index(j,i), &
	1422	surf_usm_h%end_index(j,i)
[3704]	1423	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
[3522]	1424	ENDDO
	1425	ENDDO
	1426	ENDIF
	1427
[3704]	1428	CASE ( 'rlus' )
[3522]	1429	IF ( radiation ) THEN
	1430	DO m = 1, vmea(l)%ns
	1431	!
	1432	!-- Surface data is only available on inner subdomains, not
	1433	!-- on ghost points. Hence, limit the indices.
	1434	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1435	j = MERGE( j , nyn, j > nyn )
	1436	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1437	i = MERGE( i , nxr, i > nxr )
	1438
	1439	DO mm = surf_lsm_h%start_index(j,i), &
	1440	surf_lsm_h%end_index(j,i)
[3704]	1441	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
[3522]	1442	ENDDO
	1443	DO mm = surf_usm_h%start_index(j,i), &
	1444	surf_usm_h%end_index(j,i)
[3704]	1445	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
[3522]	1446	ENDDO
	1447	ENDDO
	1448	ENDIF
	1449
[3704]	1450	CASE ( 'rlds' )
[3522]	1451	IF ( radiation ) THEN
	1452	DO m = 1, vmea(l)%ns
	1453	!
	1454	!-- Surface data is only available on inner subdomains, not
	1455	!-- on ghost points. Hence, limit the indices.
	1456	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1457	j = MERGE( j , nyn, j > nyn )
	1458	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1459	i = MERGE( i , nxr, i > nxr )
	1460
	1461	DO mm = surf_lsm_h%start_index(j,i), &
	1462	surf_lsm_h%end_index(j,i)
[3704]	1463	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
[3522]	1464	ENDDO
	1465	DO mm = surf_usm_h%start_index(j,i), &
	1466	surf_usm_h%end_index(j,i)
[3704]	1467	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
[3522]	1468	ENDDO
	1469	ENDDO
	1470	ENDIF
[3704]	1471
	1472	CASE ( 'rsd' )
	1473	IF ( radiation ) THEN
	1474	DO m = 1, vmea(l)%ns
	1475	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
	1476	j = vmea(l)%j(m)
	1477	i = vmea(l)%i(m)
	1478
	1479	vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
	1480	ENDDO
	1481	ENDIF
	1482
	1483	CASE ( 'rsu' )
	1484	IF ( radiation ) THEN
	1485	DO m = 1, vmea(l)%ns
	1486	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
	1487	j = vmea(l)%j(m)
	1488	i = vmea(l)%i(m)
	1489
	1490	vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
	1491	ENDDO
	1492	ENDIF
	1493
	1494	CASE ( 'rlu' )
	1495	IF ( radiation ) THEN
	1496	DO m = 1, vmea(l)%ns
	1497	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
	1498	j = vmea(l)%j(m)
	1499	i = vmea(l)%i(m)
	1500
	1501	vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
	1502	ENDDO
	1503	ENDIF
	1504
	1505	CASE ( 'rld' )
	1506	IF ( radiation ) THEN
	1507	DO m = 1, vmea(l)%ns
	1508	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
	1509	j = vmea(l)%j(m)
	1510	i = vmea(l)%i(m)
	1511
	1512	vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
	1513	ENDDO
	1514	ENDIF
	1515
	1516	CASE ( 'rsddif' )
	1517	IF ( radiation ) THEN
	1518	DO m = 1, vmea(l)%ns
	1519	j = vmea(l)%j(m)
	1520	i = vmea(l)%i(m)
	1521
	1522	vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
	1523	ENDDO
	1524	ENDIF
	1525
	1526	CASE ( 't_soil' )
	1527	DO m = 1, vmea(l)%ns_soil
	1528	i = vmea(l)%i_soil(m)
	1529	j = vmea(l)%j_soil(m)
	1530	k = vmea(l)%k_soil(m)
	1531
	1532	match_lsm = surf_lsm_h%start_index(j,i) <= &
	1533	surf_lsm_h%end_index(j,i)
	1534	match_usm = surf_usm_h%start_index(j,i) <= &
	1535	surf_usm_h%end_index(j,i)
	1536
	1537	IF ( match_lsm ) THEN
	1538	mm = surf_lsm_h%start_index(j,i)
	1539	vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
	1540	ENDIF
	1541
	1542	IF ( match_usm ) THEN
	1543	mm = surf_usm_h%start_index(j,i)
	1544	vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
	1545	ENDIF
	1546	ENDDO
	1547
	1548	CASE ( 'm_soil' )
	1549	DO m = 1, vmea(l)%ns_soil
	1550	i = vmea(l)%i_soil(m)
	1551	j = vmea(l)%j_soil(m)
	1552	k = vmea(l)%k_soil(m)
	1553
	1554	match_lsm = surf_lsm_h%start_index(j,i) <= &
	1555	surf_lsm_h%end_index(j,i)
	1556
	1557	IF ( match_lsm ) THEN
	1558	mm = surf_lsm_h%start_index(j,i)
	1559	vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
	1560	ENDIF
	1561
	1562	ENDDO
[3522]	1563	!
	1564	!-- More will follow ...
[3704]	1565
	1566	!
	1567	!-- No match found - just set a fill value
	1568	CASE DEFAULT
	1569	vmea(l)%measured_vars(:,n) = vmea(l)%fillout
[3522]	1570	END SELECT
	1571
[3494]	1572	ENDDO
[3434]	1573
	1574	ENDDO
[3704]	1575
[3471]	1576	END SUBROUTINE vm_sampling
[3434]	1577
	1578
[3471]	1579	END MODULE virtual_measurement_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |