Home

Context Navigation

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

Last change on this file since 4744 was 4707, checked in by suehring, 4 years ago
Bugfix for previous commit + minor formatting adjustments
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 157.9 KB

Rev	Line
[3471]	1	!> @virtual_measurement_mod.f90
[4498]	2	!--------------------------------------------------------------------------------------------------!
[3434]	3	! This file is part of the PALM model system.
	4	!
[4498]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[3434]	8	!
[4498]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[3434]	12	!
[4498]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[3434]	15	!
[4498]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
	17	!--------------------------------------------------------------------------------------------------!
[3434]	18	!
[4498]	19	!
[3434]	20	! Current revisions:
	21	! -----------------
[4641]	22	!
	23	!
[3705]	24	! Former revisions:
	25	! -----------------
	26	! $Id: virtual_measurement_mod.f90 4707 2020-09-28 14:25:28Z suehring $
[4707]	27	! Bugfix for previous commit + minor formatting adjustments
	28	!
	29	! 4706 2020-09-28 13:15:23Z suehring
[4706]	30	! Revise setting of measurement height above ground for trajectory measurements
	31	!
	32	! 4695 2020-09-24 11:30:03Z suehring
[4695]	33	! Introduce additional namelist parameters to further customize sampling in the horizontal and
	34	! vertical surroundings of the original observation coordinates
	35	!
	36	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	37	! Implementation of downward facing USM and LSM surfaces
	38	!
	39	! 4645 2020-08-24 13:55:58Z suehring
[4645]	40	! Bugfix in output of E_UTM_soil coordinate
	41	!
	42	! 4642 2020-08-13 15:47:33Z suehring
[4642]	43	! Do not set attribute bounds for time variable, as it refers to time_bounds which is not defined
	44	! for non-aggregated quantities (according to data standard)
	45	!
	46	! 4641 2020-08-13 09:57:07Z suehring
[4641]	47	! - To be in agreement with (UC)2 data standard do not list the measured variables in attribute
	48	! data_content but simply set 'airmeteo'
	49	! - Bugfix in setting long_name attribute for variable t_va and for global attribute creation_time
	50	!
	51	! 4536 2020-05-17 17:24:13Z raasch
[4536]	52	! bugfix: preprocessor directive adjusted
	53	!
	54	! 4504 2020-04-20 12:11:24Z raasch
[4498]	55	! file re-formatted to follow the PALM coding standard
	56	!
	57	! 4481 2020-03-31 18:55:54Z maronga
[4444]	58	! bugfix: cpp-directives for serial mode added
[4498]	59	!
[4444]	60	! 4438 2020-03-03 20:49:28Z suehring
[4438]	61	! Add cpu-log points
[4498]	62	!
[4438]	63	! 4422 2020-02-24 22:45:13Z suehring
[4422]	64	! Missing trim()
[4498]	65	!
[4422]	66	! 4408 2020-02-14 10:04:39Z gronemeier
[4421]	67	! - Output of character string station_name after DOM has been enabled to
	68	! output character variables
[4422]	69	! - Bugfix, missing coupling_char statement when opening the input file
[4498]	70	!
[4421]	71	! 4408 2020-02-14 10:04:39Z gronemeier
[4408]	72	! write fill_value attribute
	73	!
	74	! 4406 2020-02-13 20:06:29Z knoop
[4406]	75	! Bugix: removed oro_rel wrong loop bounds and removed unnecessary restart method
[4408]	76	!
[4406]	77	! 4400 2020-02-10 20:32:41Z suehring
[4400]	78	! Revision of the module:
	79	! - revised input from NetCDF setup file
	80	! - parallel NetCDF output via data-output module ( Tobias Gronemeier )
	81	! - variable attributes added
	82	! - further variables defined
	83	!
	84	! 4346 2019-12-18 11:55:56Z motisi
[4346]	85	! Introduction of wall_flags_total_0, which currently sets bits based on static
	86	! topography information used in wall_flags_static_0
[4400]	87	!
[4346]	88	! 4329 2019-12-10 15:46:36Z motisi
[4329]	89	! Renamed wall_flags_0 to wall_flags_static_0
[4400]	90	!
[4329]	91	! 4226 2019-09-10 17:03:24Z suehring
[4226]	92	! Netcdf input routine for dimension length renamed
[4400]	93	!
[4226]	94	! 4182 2019-08-22 15:20:23Z scharf
[4182]	95	! Corrected "Former revisions" section
[4400]	96	!
[4182]	97	! 4168 2019-08-16 13:50:17Z suehring
[4168]	98	! Replace function get_topography_top_index by topo_top_ind
[4400]	99	!
[4168]	100	! 3988 2019-05-22 11:32:37Z kanani
[3988]	101	! Add variables to enable steering of output interval for virtual measurements
[4400]	102	!
[3988]	103	! 3913 2019-04-17 15:12:28Z gronemeier
[3913]	104	! Bugfix: rotate positions of measurements before writing them into file
[4400]	105	!
[3913]	106	! 3910 2019-04-17 11:46:56Z suehring
[3910]	107	! Bugfix in rotation of UTM coordinates
[4400]	108	!
[3910]	109	! 3904 2019-04-16 18:22:51Z gronemeier
[3904]	110	! Rotate coordinates of stations by given rotation_angle
[4400]	111	!
[3904]	112	! 3876 2019-04-08 18:41:49Z knoop
[3855]	113	! Remove print statement
[4400]	114	!
[3855]	115	! 3854 2019-04-02 16:59:33Z suehring
[3854]	116	! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
[4400]	117	! nspec by nvar
	118	!
[3833]	119	! 3766 2019-02-26 16:23:41Z raasch
[3766]	120	! unused variables removed
[4400]	121	!
[3766]	122	! 3718 2019-02-06 11:08:28Z suehring
[3718]	123	! Adjust variable name connections between UC2 and chemistry variables
[4400]	124	!
[3718]	125	! 3717 2019-02-05 17:21:16Z suehring
[3717]	126	! Additional check + error numbers adjusted
[4400]	127	!
[3717]	128	! 3706 2019-01-29 20:02:26Z suehring
[3706]	129	! unused variables removed
[4400]	130	!
[3706]	131	! 3705 2019-01-29 19:56:39Z suehring
[3704]	132	! - initialization revised
	133	! - binary data output
	134	! - list of allowed variables extended
[4400]	135	!
[3705]	136	! 3704 2019-01-29 19:51:41Z suehring
[3522]	137	! Sampling of variables
[4400]	138	!
[4182]	139	! 3473 2018-10-30 20:50:15Z suehring
	140	! Initial revision
[3434]	141	!
[4182]	142	! Authors:
	143	! --------
	144	! @author Matthias Suehring
[4400]	145	! @author Tobias Gronemeier
[4182]	146	!
[3434]	147	! Description:
	148	! ------------
[4498]	149	!> The module acts as an interface between 'real-world' observations and model simulations.
	150	!> Virtual measurements will be taken in the model at the coordinates representative for the
	151	!> 'real-world' observation coordinates. More precisely, coordinates and measured quanties will be
	152	!> read from a NetCDF file which contains all required information. In the model, the same
	153	!> quantities (as long as all the required components are switched-on) will be sampled at the
	154	!> respective positions and output into an extra file, which allows for straight-forward comparison
	155	!> of model results with observations.
	156	!--------------------------------------------------------------------------------------------------!
[3471]	157	MODULE virtual_measurement_mod
[3434]	158
[4498]	159	USE arrays_3d, &
	160	ONLY: dzw, &
	161	exner, &
	162	hyp, &
	163	q, &
	164	ql, &
	165	pt, &
	166	rho_air, &
	167	u, &
	168	v, &
	169	w, &
	170	zu, &
[4400]	171	zw
[3434]	172
[4498]	173	USE basic_constants_and_equations_mod, &
	174	ONLY: convert_utm_to_geographic, &
	175	degc_to_k, &
	176	magnus, &
	177	pi, &
[4400]	178	rd_d_rv
[3904]	179
[4498]	180	USE chem_gasphase_mod, &
[3833]	181	ONLY: nvar
[3522]	182
[4498]	183	USE chem_modules, &
[3522]	184	ONLY: chem_species
[4400]	185
[4498]	186	USE control_parameters, &
	187	ONLY: air_chemistry, &
	188	coupling_char, &
	189	dz, &
	190	end_time, &
	191	humidity, &
	192	message_string, &
	193	neutral, &
	194	origin_date_time, &
	195	rho_surface, &
	196	surface_pressure, &
	197	time_since_reference_point, &
[4406]	198	virtual_measurement
[3434]	199
[4498]	200	USE cpulog, &
	201	ONLY: cpu_log, &
[4438]	202	log_point_s
[4400]	203
	204	USE data_output_module
	205
[4498]	206	USE grid_variables, &
	207	ONLY: ddx, &
	208	ddy, &
	209	dx, &
[4400]	210	dy
[3434]	211
[4498]	212	USE indices, &
	213	ONLY: nbgp, &
	214	nzb, &
	215	nzt, &
	216	nxl, &
	217	nxlg, &
	218	nxr, &
	219	nxrg, &
	220	nys, &
	221	nysg, &
	222	nyn, &
	223	nyng, &
	224	topo_top_ind, &
[4346]	225	wall_flags_total_0
[3434]	226
	227	USE kinds
[4400]	228
[4498]	229	USE netcdf_data_input_mod, &
	230	ONLY: close_input_file, &
	231	coord_ref_sys, &
	232	crs_list, &
	233	get_attribute, &
	234	get_dimension_length, &
	235	get_variable, &
	236	init_model, &
	237	input_file_atts, &
	238	input_file_vm, &
	239	input_pids_static, &
	240	input_pids_vm, &
	241	inquire_fill_value, &
	242	open_read_file, &
[4400]	243	pids_id
	244
[3704]	245	USE pegrid
[4400]	246
[4498]	247	USE surface_mod, &
	248	ONLY: surf_lsm_h, &
[4400]	249	surf_usm_h
	250
[4498]	251	USE land_surface_model_mod, &
	252	ONLY: m_soil_h, &
	253	nzb_soil, &
	254	nzt_soil, &
	255	t_soil_h, &
[4400]	256	zs
	257
[4498]	258	USE radiation_model_mod, &
	259	ONLY: rad_lw_in, &
	260	rad_lw_out, &
	261	rad_sw_in, &
	262	rad_sw_in_diff, &
	263	rad_sw_out, &
[4400]	264	radiation_scheme
	265
[4498]	266	USE urban_surface_mod, &
	267	ONLY: nzb_wall, &
	268	nzt_wall, &
[4400]	269	t_wall_h
[3434]	270
	271
	272	IMPLICIT NONE
[4400]	273
[3704]	274	TYPE virt_general
[4498]	275	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
[3704]	276	END TYPE virt_general
[3434]	277
[4400]	278	TYPE virt_var_atts
[4498]	279	CHARACTER(LEN=100) :: coordinates !< defined longname of the variable
	280	CHARACTER(LEN=100) :: grid_mapping !< defined longname of the variable
	281	CHARACTER(LEN=100) :: long_name !< defined longname of the variable
	282	CHARACTER(LEN=100) :: name !< variable name
	283	CHARACTER(LEN=100) :: standard_name !< defined standard name of the variable
	284	CHARACTER(LEN=100) :: units !< unit of the output variable
[4400]	285
[4498]	286	REAL(wp) :: fill_value = -9999.0 !< _FillValue attribute
[4400]	287	END TYPE virt_var_atts
	288
[3434]	289	TYPE virt_mea
[4498]	290	CHARACTER(LEN=100) :: feature_type !< type of the real-world measurement
	291	CHARACTER(LEN=100) :: feature_type_out = 'timeSeries' !< type of the virtual measurement
	292	!< (all will be timeSeries, even trajectories)
	293	CHARACTER(LEN=100) :: nc_filename !< name of the NetCDF output file for the station
	294	CHARACTER(LEN=100) :: site !< name of the measurement site
[4400]	295
[4498]	296	CHARACTER(LEN=1000) :: data_content = REPEAT(' ', 1000) !< string of measured variables (data output only)
[4400]	297
[4498]	298	INTEGER(iwp) :: end_coord_a = 0 !< end coordinate in NetCDF file for local atmosphere observations
	299	INTEGER(iwp) :: end_coord_s = 0 !< end coordinate in NetCDF file for local soil observations
	300	INTEGER(iwp) :: file_time_index = 0 !< time index in NetCDF output file
	301	INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements
	302	INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements
	303	INTEGER(iwp) :: n_tr_st !< number of trajectories / station of a measurement
	304	INTEGER(iwp) :: nmeas !< number of measured variables (atmosphere + soil)
	305	INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements
	306	INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
	307	INTEGER(iwp) :: start_coord_a = 0 !< start coordinate in NetCDF file for local atmosphere observations
	308	INTEGER(iwp) :: start_coord_s = 0 !< start coordinate in NetCDF file for local soil observations
[4400]	309
[4498]	310	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory
	311	!< or station that are no _FillValues
[4400]	312
[3704]	313	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
	314	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
	315	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
[4400]	316
[3704]	317	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction
	318	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction
	319	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction
[4400]	320
[4498]	321	LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled
	322	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
	323	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
	324	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
[4400]	325
[4498]	326	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
	327	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
	328	REAL(wp) :: fill_zar !< fill value for heigth coordinates in case of missing values
	329	REAL(wp) :: fillout = -9999.0 !< fill value for output in case an observation is taken e.g. from inside a building
	330	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
	331	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
[4400]	332
[4706]	333	REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil
	334	REAL(wp), DIMENSION(:), ALLOCATABLE :: zar !< measurement height above ground level
[4400]	335
	336	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
	337	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables
	338
	339	TYPE( virt_var_atts ), DIMENSION(:), ALLOCATABLE :: var_atts !< variable attributes
[3434]	340	END TYPE virt_mea
	341
[4707]	342	CHARACTER(LEN=5) :: char_eutm = 'E_UTM' !< dimension name for UTM coordinate easting
	343	CHARACTER(LEN=11) :: char_feature = 'featureType' !< attribute name for feature type
[4400]	344
	345	! This need to be generalized
[4408]	346	CHARACTER(LEN=10) :: char_fill = '_FillValue' !< attribute name for fill value
[4707]	347	CHARACTER(LEN=6) :: char_height = 'height' !< attribute name indicating height above surface (for trajectories only)
[4400]	348	CHARACTER(LEN=9) :: char_long = 'long_name' !< attribute name for long_name
[4707]	349	CHARACTER(LEN=18) :: char_mv = 'measured_variables' !< variable name for the array with the measured variable names
	350	CHARACTER(LEN=5) :: char_nutm = 'N_UTM' !< dimension name for UTM coordinate northing
	351	CHARACTER(LEN=18) :: char_numstations = 'number_of_stations' !< attribute name for number of stations
	352	CHARACTER(LEN=8) :: char_origx = 'origin_x' !< attribute name for station coordinate in x
	353	CHARACTER(LEN=8) :: char_origy = 'origin_y' !< attribute name for station coordinate in y
	354	CHARACTER(LEN=4) :: char_site = 'site' !< attribute name for site name
	355	CHARACTER(LEN=11) :: char_soil = 'soil_sample' !< attribute name for soil sampling indication
[4498]	356	CHARACTER(LEN=13) :: char_standard = 'standard_name' !< attribute name for standard_name
[4707]	357	CHARACTER(LEN=9) :: char_station_h = 'station_h' !< variable name indicating height of the site
[4498]	358	CHARACTER(LEN=5) :: char_unit = 'units' !< attribute name for standard_name
[4707]	359	CHARACTER(LEN=1) :: char_zar = 'z' !< attribute name indicating height above reference level
[3434]	360	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
	361	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
	362	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
[4400]	363
[4498]	364	CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ 't_soil', & !< list of soil variables
	365	'm_soil', &
	366	'lwc ', &
	367	'lwcs ', &
	368	'smp ' /)
[4400]	369
[4498]	370	CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ 'mcpm1 ', 'PM1 ', &
	371	'mcpm2p5 ', 'PM2.5 ', &
	372	'mcpm10 ', 'PM10 ', &
	373	'mfno2 ', 'NO2 ', &
	374	'mfno ', 'NO ', &
	375	'mcno2 ', 'NO2 ', &
	376	'mcno ', 'NO ', &
	377	'tro3 ', 'O3 ' &
	378	/), (/ 2, 8 /) )
[3434]	379
[4498]	380	INTEGER(iwp) :: maximum_name_length = 32 !< maximum name length of station names
	381	INTEGER(iwp) :: ntimesteps !< number of timesteps defined in NetCDF output file
[4695]	382	INTEGER(iwp) :: off_pr = 1 !< number of neighboring grid points (in horizontal direction) where virtual profile
[4498]	383	!< measurements shall be taken, in addition to the given coordinates in the driver
[4695]	384	INTEGER(iwp) :: off_pr_z = 0 !< number of additional grid points (in each upwardd and downward direction) where
	385	!< virtual profile measurements shall be taken, in addition to the given z coordinates in the driver
	386	INTEGER(iwp) :: off_ts = 1 !< number of neighboring grid points (in horizontal direction) where virtual profile
[4498]	387	!< measurements shall be taken, in addition to the given coordinates in the driver
[4695]	388	INTEGER(iwp) :: off_ts_z = 50 !< number of additional grid points (in each upwardd and downward direction) where
	389	!< virtual profile measurements shall be taken, in addition to the given z coordinates in the driver
	390	INTEGER(iwp) :: off_tr = 1 !< number of neighboring grid points (in horizontal direction) where virtual profile
[4498]	391	!< measurements shall be taken, in addition to the given coordinates in the driver
[4695]	392	INTEGER(iwp) :: off_tr_z = 1 !< number of additional grid points (in each upwardd and downward direction) where
	393	!< virtual profile measurements shall be taken, in addition to the given z coordinates in the driver
[4498]	394	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
	395	LOGICAL :: initial_write_coordinates = .FALSE. !< flag indicating a global attribute
	396	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
[3988]	397
[4498]	398	REAL(wp) :: dt_virtual_measurement = 0.0_wp !< sampling interval
[4400]	399	REAL(wp) :: time_virtual_measurement = 0.0_wp !< time since last sampling
[4498]	400	REAL(wp) :: vm_time_start = 0.0 !< time after which sampling shall start
[4400]	401
[4498]	402	TYPE( virt_general ) :: vmea_general !< data structure which encompasses global variables
	403	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< data structure containing station-specific variables
[4400]	404
[3434]	405	INTERFACE vm_check_parameters
	406	MODULE PROCEDURE vm_check_parameters
	407	END INTERFACE vm_check_parameters
[4400]	408
[3704]	409	INTERFACE vm_data_output
	410	MODULE PROCEDURE vm_data_output
	411	END INTERFACE vm_data_output
[4400]	412
[3434]	413	INTERFACE vm_init
	414	MODULE PROCEDURE vm_init
	415	END INTERFACE vm_init
[4400]	416
	417	INTERFACE vm_init_output
	418	MODULE PROCEDURE vm_init_output
	419	END INTERFACE vm_init_output
	420
[3434]	421	INTERFACE vm_parin
	422	MODULE PROCEDURE vm_parin
	423	END INTERFACE vm_parin
[4400]	424
[3434]	425	INTERFACE vm_sampling
	426	MODULE PROCEDURE vm_sampling
	427	END INTERFACE vm_sampling
	428
	429	SAVE
	430
	431	PRIVATE
	432
	433	!
	434	!-- Public interfaces
[4498]	435	PUBLIC vm_check_parameters, &
	436	vm_data_output, &
	437	vm_init, &
	438	vm_init_output, &
	439	vm_parin, &
[4400]	440	vm_sampling
[3434]	441
	442	!
	443	!-- Public variables
[4498]	444	PUBLIC dt_virtual_measurement, &
	445	time_virtual_measurement, &
	446	vmea, &
	447	vmea_general, &
[4400]	448	vm_time_start
[3434]	449
	450	CONTAINS
	451
	452
[4498]	453	!--------------------------------------------------------------------------------------------------!
[3434]	454	! Description:
	455	! ------------
[3471]	456	!> Check parameters for virtual measurement module
[4498]	457	!--------------------------------------------------------------------------------------------------!
[3434]	458	SUBROUTINE vm_check_parameters
	459
[4400]	460	IF ( .NOT. virtual_measurement ) RETURN
[3434]	461	!
[4400]	462	!-- Virtual measurements require a setup file.
	463	IF ( .NOT. input_pids_vm ) THEN
[4498]	464	message_string = 'If virtual measurements are taken, a setup input ' // &
[3717]	465	'file for the site locations is mandatory.'
	466	CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
[4400]	467	ENDIF
[3717]	468	!
[3434]	469	!-- In case virtual measurements are taken, a static input file is required.
[4498]	470	!-- This is because UTM coordinates for the PALM domain origin are required for correct mapping of
	471	!-- the measurements.
[3434]	472	!-- ToDo: Revise this later and remove this requirement.
[4400]	473	IF ( .NOT. input_pids_static ) THEN
[4498]	474	message_string = 'If virtual measurements are taken, a static input file is mandatory.'
[3717]	475	CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
[3434]	476	ENDIF
[4400]	477
	478	#if !defined( __netcdf4_parallel )
	479	!
	480	!-- In case of non-parallel NetCDF the virtual measurement output is not
	481	!-- working. This is only designed for parallel NetCDF.
[4498]	482	message_string = 'If virtual measurements are taken, parallel NetCDF is required.'
[4400]	483	CALL message( 'vm_check_parameters', 'PA0708', 1, 2, 0, 6, 0 )
	484	#endif
	485	!
[4498]	486	!-- Check if the given number of neighboring grid points do not exceed the number
[4400]	487	!-- of ghost points.
[4498]	488	IF ( off_pr > nbgp - 1 .OR. off_ts > nbgp - 1 .OR. off_tr > nbgp - 1 ) THEN
	489	WRITE(message_string,*) &
	490	'If virtual measurements are taken, the number ' // &
	491	'of surrounding grid points must not be larger ' // &
[4504]	492	'than the number of ghost points - 1, which is: ', nbgp - 1
[4400]	493	CALL message( 'vm_check_parameters', 'PA0705', 1, 2, 0, 6, 0 )
	494	ENDIF
[4406]	495
	496	IF ( dt_virtual_measurement <= 0.0 ) THEN
	497	message_string = 'dt_virtual_measurement must be > 0.0'
[4400]	498	CALL message( 'check_parameters', 'PA0706', 1, 2, 0, 6, 0 )
	499	ENDIF
	500
[3434]	501	END SUBROUTINE vm_check_parameters
[4408]	502
[4498]	503	!--------------------------------------------------------------------------------------------------!
[3434]	504	! Description:
	505	! ------------
[4498]	506	!> Subroutine defines variable attributes according to UC2 standard. Note, later this list can be
	507	!> moved to the data-output module where it can be re-used also for other output.
	508	!--------------------------------------------------------------------------------------------------!
	509	SUBROUTINE vm_set_attributes( output_variable )
[4400]	510
[4498]	511	TYPE( virt_var_atts ), INTENT(INOUT) :: output_variable !< data structure with attributes that need to be set
[4400]	512
[4498]	513	output_variable%long_name = 'none'
	514	output_variable%standard_name = 'none'
	515	output_variable%units = 'none'
	516	output_variable%coordinates = 'lon lat E_UTM N_UTM x y z time station_name'
	517	output_variable%grid_mapping = 'crs'
[4400]	518
[4498]	519	SELECT CASE ( TRIM( output_variable%name ) )
[4400]	520
[4498]	521	CASE ( 'u' )
	522	output_variable%long_name = 'u wind component'
	523	output_variable%units = 'm s-1'
[4400]	524
[4498]	525	CASE ( 'ua' )
	526	output_variable%long_name = 'eastward wind'
	527	output_variable%standard_name = 'eastward_wind'
	528	output_variable%units = 'm s-1'
[4400]	529
[4498]	530	CASE ( 'v' )
	531	output_variable%long_name = 'v wind component'
	532	output_variable%units = 'm s-1'
[4400]	533
[4498]	534	CASE ( 'va' )
	535	output_variable%long_name = 'northward wind'
	536	output_variable%standard_name = 'northward_wind'
	537	output_variable%units = 'm s-1'
[4400]	538
[4498]	539	CASE ( 'w' )
	540	output_variable%long_name = 'w wind component'
	541	output_variable%standard_name = 'upward_air_velocity'
	542	output_variable%units = 'm s-1'
[4400]	543
[4498]	544	CASE ( 'wspeed' )
	545	output_variable%long_name = 'wind speed'
	546	output_variable%standard_name = 'wind_speed'
	547	output_variable%units = 'm s-1'
[4400]	548
[4498]	549	CASE ( 'wdir' )
	550	output_variable%long_name = 'wind from direction'
	551	output_variable%standard_name = 'wind_from_direction'
	552	output_variable%units = 'degrees'
[4400]	553
[4498]	554	CASE ( 'theta' )
	555	output_variable%long_name = 'air potential temperature'
	556	output_variable%standard_name = 'air_potential_temperature'
	557	output_variable%units = 'K'
[4400]	558
[4498]	559	CASE ( 'utheta' )
	560	output_variable%long_name = 'eastward kinematic sensible heat flux in air'
	561	output_variable%units = 'K m s-1'
[4400]	562
[4498]	563	CASE ( 'vtheta' )
	564	output_variable%long_name = 'northward kinematic sensible heat flux in air'
	565	output_variable%units = 'K m s-1'
[4400]	566
[4498]	567	CASE ( 'wtheta' )
	568	output_variable%long_name = 'upward kinematic sensible heat flux in air'
	569	output_variable%units = 'K m s-1'
[4400]	570
[4498]	571	CASE ( 'ta' )
	572	output_variable%long_name = 'air temperature'
	573	output_variable%standard_name = 'air_temperature'
	574	output_variable%units = 'degree_C'
[4400]	575
[4641]	576	CASE ( 't_va' )
[4498]	577	output_variable%long_name = 'virtual acoustic temperature'
	578	output_variable%units = 'K'
[4400]	579
[4498]	580	CASE ( 'haa' )
	581	output_variable%long_name = 'absolute atmospheric humidity'
	582	output_variable%units = 'kg m-3'
[4400]	583
[4498]	584	CASE ( 'hus' )
	585	output_variable%long_name = 'specific humidity'
	586	output_variable%standard_name = 'specific_humidity'
	587	output_variable%units = 'kg kg-1'
[4400]	588
[4498]	589	CASE ( 'hur' )
	590	output_variable%long_name = 'relative humidity'
	591	output_variable%standard_name = 'relative_humidity'
	592	output_variable%units = '1'
[4400]	593
[4498]	594	CASE ( 'rlu' )
	595	output_variable%long_name = 'upwelling longwave flux in air'
	596	output_variable%standard_name = 'upwelling_longwave_flux_in_air'
	597	output_variable%units = 'W m-2'
[4400]	598
[4498]	599	CASE ( 'rlus' )
	600	output_variable%long_name = 'surface upwelling longwave flux in air'
	601	output_variable%standard_name = 'surface_upwelling_longwave_flux_in_air'
	602	output_variable%units = 'W m-2'
[4400]	603
[4498]	604	CASE ( 'rld' )
	605	output_variable%long_name = 'downwelling longwave flux in air'
	606	output_variable%standard_name = 'downwelling_longwave_flux_in_air'
	607	output_variable%units = 'W m-2'
[4400]	608
[4498]	609	CASE ( 'rsddif' )
	610	output_variable%long_name = 'diffuse downwelling shortwave flux in air'
	611	output_variable%standard_name = 'diffuse_downwelling_shortwave_flux_in_air'
	612	output_variable%units = 'W m-2'
[4400]	613
[4498]	614	CASE ( 'rsd' )
	615	output_variable%long_name = 'downwelling shortwave flux in air'
	616	output_variable%standard_name = 'downwelling_shortwave_flux_in_air'
	617	output_variable%units = 'W m-2'
[4400]	618
[4498]	619	CASE ( 'rnds' )
	620	output_variable%long_name = 'surface net downward radiative flux'
	621	output_variable%standard_name = 'surface_net_downward_radiative_flux'
	622	output_variable%units = 'W m-2'
[4400]	623
[4498]	624	CASE ( 'rsu' )
	625	output_variable%long_name = 'upwelling shortwave flux in air'
	626	output_variable%standard_name = 'upwelling_shortwave_flux_in_air'
	627	output_variable%units = 'W m-2'
[4400]	628
[4498]	629	CASE ( 'rsus' )
	630	output_variable%long_name = 'surface upwelling shortwave flux in air'
	631	output_variable%standard_name = 'surface_upwelling_shortwave_flux_in_air'
	632	output_variable%units = 'W m-2'
[4400]	633
[4498]	634	CASE ( 'rsds' )
	635	output_variable%long_name = 'surface downwelling shortwave flux in air'
	636	output_variable%standard_name = 'surface_downwelling_shortwave_flux_in_air'
	637	output_variable%units = 'W m-2'
[4400]	638
[4498]	639	CASE ( 'hfss' )
	640	output_variable%long_name = 'surface upward sensible heat flux'
	641	output_variable%standard_name = 'surface_upward_sensible_heat_flux'
	642	output_variable%units = 'W m-2'
[4400]	643
[4498]	644	CASE ( 'hfls' )
	645	output_variable%long_name = 'surface upward latent heat flux'
	646	output_variable%standard_name = 'surface_upward_latent_heat_flux'
	647	output_variable%units = 'W m-2'
[4400]	648
[4498]	649	CASE ( 'ts' )
	650	output_variable%long_name = 'surface temperature'
	651	output_variable%standard_name = 'surface_temperature'
	652	output_variable%units = 'K'
[4400]	653
[4498]	654	CASE ( 'thetas' )
	655	output_variable%long_name = 'surface layer temperature scale'
	656	output_variable%units = 'K'
[4400]	657
[4498]	658	CASE ( 'us' )
	659	output_variable%long_name = 'friction velocity'
	660	output_variable%units = 'm s-1'
[4400]	661
[4498]	662	CASE ( 'uw' )
	663	output_variable%long_name = 'upward eastward kinematic momentum flux in air'
	664	output_variable%units = 'm2 s-2'
[4400]	665
[4498]	666	CASE ( 'vw' )
	667	output_variable%long_name = 'upward northward kinematic momentum flux in air'
	668	output_variable%units = 'm2 s-2'
[4400]	669
[4498]	670	CASE ( 'uv' )
	671	output_variable%long_name = 'eastward northward kinematic momentum flux in air'
	672	output_variable%units = 'm2 s-2'
[4400]	673
[4498]	674	CASE ( 'plev' )
	675	output_variable%long_name = 'air pressure'
	676	output_variable%standard_name = 'air_pressure'
	677	output_variable%units = 'Pa'
[4400]	678
[4498]	679	CASE ( 'm_soil' )
	680	output_variable%long_name = 'soil moisture volumetric'
	681	output_variable%units = 'm3 m-3'
[4400]	682
[4498]	683	CASE ( 't_soil' )
	684	output_variable%long_name = 'soil temperature'
	685	output_variable%standard_name = 'soil_temperature'
	686	output_variable%units = 'degree_C'
[4400]	687
[4498]	688	CASE ( 'hfdg' )
	689	output_variable%long_name = 'downward heat flux at ground level in soil'
	690	output_variable%standard_name = 'downward_heat_flux_at_ground_level_in_soil'
	691	output_variable%units = 'W m-2'
[4400]	692
[4498]	693	CASE ( 'hfds' )
	694	output_variable%long_name = 'downward heat flux in soil'
	695	output_variable%standard_name = 'downward_heat_flux_in_soil'
	696	output_variable%units = 'W m-2'
[4400]	697
[4498]	698	CASE ( 'hfla' )
	699	output_variable%long_name = 'upward latent heat flux in air'
	700	output_variable%standard_name = 'upward_latent_heat_flux_in_air'
	701	output_variable%units = 'W m-2'
[4400]	702
[4498]	703	CASE ( 'hfsa' )
	704	output_variable%long_name = 'upward latent heat flux in air'
	705	output_variable%standard_name = 'upward_sensible_heat_flux_in_air'
	706	output_variable%units = 'W m-2'
[4400]	707
[4498]	708	CASE ( 'jno2' )
	709	output_variable%long_name = 'photolysis rate of nitrogen dioxide'
	710	output_variable%standard_name = 'photolysis_rate_of_nitrogen_dioxide'
	711	output_variable%units = 's-1'
[4400]	712
[4498]	713	CASE ( 'lwcs' )
	714	output_variable%long_name = 'liquid water content of soil layer'
	715	output_variable%standard_name = 'liquid_water_content_of_soil_layer'
	716	output_variable%units = 'kg m-2'
[4400]	717
[4498]	718	CASE ( 'lwp' )
	719	output_variable%long_name = 'liquid water path'
	720	output_variable%standard_name = 'atmosphere_mass_content_of_cloud_liquid_water'
	721	output_variable%units = 'kg m-2'
[4400]	722
[4498]	723	CASE ( 'ps' )
	724	output_variable%long_name = 'surface air pressure'
	725	output_variable%standard_name = 'surface_air_pressure'
	726	output_variable%units = 'hPa'
[4400]	727
[4498]	728	CASE ( 'pswrtg' )
	729	output_variable%long_name = 'platform speed wrt ground'
	730	output_variable%standard_name = 'platform_speed_wrt_ground'
	731	output_variable%units = 'm s-1'
[4400]	732
[4498]	733	CASE ( 'pswrta' )
	734	output_variable%long_name = 'platform speed wrt air'
	735	output_variable%standard_name = 'platform_speed_wrt_air'
	736	output_variable%units = 'm s-1'
[4400]	737
[4498]	738	CASE ( 'pwv' )
	739	output_variable%long_name = 'water vapor partial pressure in air'
	740	output_variable%standard_name = 'water_vapor_partial_pressure_in_air'
	741	output_variable%units = 'hPa'
[4400]	742
[4498]	743	CASE ( 'ssdu' )
	744	output_variable%long_name = 'duration of sunshine'
	745	output_variable%standard_name = 'duration_of_sunshine'
	746	output_variable%units = 's'
[4400]	747
[4498]	748	CASE ( 't_lw' )
	749	output_variable%long_name = 'land water temperature'
	750	output_variable%units = 'degree_C'
[4400]	751
[4498]	752	CASE ( 'tb' )
	753	output_variable%long_name = 'brightness temperature'
	754	output_variable%standard_name = 'brightness_temperature'
	755	output_variable%units = 'K'
[4400]	756
[4498]	757	CASE ( 'uqv' )
	758	output_variable%long_name = 'eastward kinematic latent heat flux in air'
	759	output_variable%units = 'g kg-1 m s-1'
[4400]	760
[4498]	761	CASE ( 'vqv' )
	762	output_variable%long_name = 'northward kinematic latent heat flux in air'
	763	output_variable%units = 'g kg-1 m s-1'
[4400]	764
[4498]	765	CASE ( 'wqv' )
	766	output_variable%long_name = 'upward kinematic latent heat flux in air'
	767	output_variable%units = 'g kg-1 m s-1'
[4400]	768
[4498]	769	CASE ( 'zcb' )
	770	output_variable%long_name = 'cloud base altitude'
	771	output_variable%standard_name = 'cloud_base_altitude'
	772	output_variable%units = 'm'
[4400]	773
[4498]	774	CASE ( 'zmla' )
	775	output_variable%long_name = 'atmosphere boundary layer thickness'
	776	output_variable%standard_name = 'atmosphere_boundary_layer_thickness'
	777	output_variable%units = 'm'
[4400]	778
[4498]	779	CASE ( 'mcpm1' )
	780	output_variable%long_name = 'mass concentration of pm1 ambient aerosol particles in air'
	781	output_variable%standard_name = 'mass_concentration_of_pm1_ambient_aerosol_particles_in_air'
	782	output_variable%units = 'kg m-3'
[4400]	783
[4498]	784	CASE ( 'mcpm10' )
	785	output_variable%long_name = 'mass concentration of pm10 ambient aerosol particles in air'
	786	output_variable%standard_name = 'mass_concentration_of_pm10_ambient_aerosol_particles_in_air'
	787	output_variable%units = 'kg m-3'
[4400]	788
[4498]	789	CASE ( 'mcpm2p5' )
	790	output_variable%long_name = 'mass concentration of pm2p5 ambient aerosol particles in air'
	791	output_variable%standard_name = 'mass_concentration_of_pm2p5_ambient_aerosol_particles_in_air'
	792	output_variable%units = 'kg m-3'
[4400]	793
[4498]	794	CASE ( 'mfno', 'mcno' )
	795	output_variable%long_name = 'mole fraction of nitrogen monoxide in air'
	796	output_variable%standard_name = 'mole_fraction_of_nitrogen_monoxide_in_air'
	797	output_variable%units = 'ppm' !'mol mol-1'
[4400]	798
[4498]	799	CASE ( 'mfno2', 'mcno2' )
	800	output_variable%long_name = 'mole fraction of nitrogen dioxide in air'
	801	output_variable%standard_name = 'mole_fraction_of_nitrogen_dioxide_in_air'
	802	output_variable%units = 'ppm' !'mol mol-1'
[4400]	803
[4641]	804	CASE ( 'ncaa' )
	805	output_variable%long_name = 'number concentration of ambient aerosol particles in air'
	806	output_variable%standard_name = 'number_concentration_of_ambient_aerosol_particles_in_air'
	807	output_variable%units = 'm-3' !'mol mol-1'
	808
[4498]	809	CASE ( 'tro3' )
	810	output_variable%long_name = 'mole fraction of ozone in air'
	811	output_variable%standard_name = 'mole_fraction_of_ozone_in_air'
	812	output_variable%units = 'ppm' !'mol mol-1'
[4400]	813
[4498]	814	CASE DEFAULT
[4400]	815
[4498]	816	END SELECT
[4400]	817
[4498]	818	END SUBROUTINE vm_set_attributes
[4400]	819
	820
[4498]	821	!--------------------------------------------------------------------------------------------------!
[4400]	822	! Description:
	823	! ------------
[3471]	824	!> Read namelist for the virtual measurement module
[4498]	825	!--------------------------------------------------------------------------------------------------!
[3434]	826	SUBROUTINE vm_parin
[4400]	827
[4498]	828	CHARACTER(LEN=80) :: line !< dummy string that contains the current line of the parameter file
[4400]	829
[4498]	830	NAMELIST /virtual_measurement_parameters/ dt_virtual_measurement, &
	831	off_pr, &
[4706]	832	off_pr_z, &
[4498]	833	off_tr, &
[4706]	834	off_tr_z, &
	835	off_ts, &
	836	off_ts_z, &
[4498]	837	use_virtual_measurement, &
[3434]	838	vm_time_start
	839
	840	line = ' '
	841	!
	842	!-- Try to find stg package
	843	REWIND ( 11 )
	844	line = ' '
[4498]	845	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
[3434]	846	READ ( 11, '(A)', END=20 ) line
	847	ENDDO
	848	BACKSPACE ( 11 )
	849
	850	!
	851	!-- Read namelist
	852	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
	853
	854	!
[3471]	855	!-- Set flag that indicates that the virtual measurement module is switched on
[3434]	856	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
	857	GOTO 20
	858
	859	10 BACKSPACE( 11 )
	860	READ( 11 , '(A)') line
	861	CALL parin_fail_message( 'virtual_measurement_parameters', line )
	862
	863	20 CONTINUE
[4400]	864
[3434]	865	END SUBROUTINE vm_parin
	866
	867
[4498]	868	!--------------------------------------------------------------------------------------------------!
[3434]	869	! Description:
	870	! ------------
[4498]	871	!> Initialize virtual measurements: read coordiante arrays and measured variables, set indicies
	872	!> indicating the measurement points, read further attributes, etc..
	873	!--------------------------------------------------------------------------------------------------!
[3434]	874	SUBROUTINE vm_init
	875
[4498]	876	CHARACTER(LEN=5) :: dum !< dummy string indicating station id
	877	CHARACTER(LEN=100), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
	878	CHARACTER(LEN=100), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
[4400]	879
[4498]	880	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
	881	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
	882	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
	883	INTEGER(iwp) :: j !< grid index of observation point in x-direction
	884	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
	885	INTEGER(iwp) :: k !< grid index of observation point in x-direction
	886	INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate
	887	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
	888	INTEGER(iwp) :: ksurf !< topography top index
	889	INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate
	890	INTEGER(iwp) :: l !< running index over all stations
	891	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
	892	INTEGER(iwp) :: ll !< running index over all measured variables in file
	893	INTEGER(iwp) :: m !< running index for surface elements
	894	INTEGER(iwp) :: n !< running index over trajectory coordinates
	895	INTEGER(iwp) :: nofill !< dummy for nofill return value (not used)
	896	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
[4695]	897	INTEGER(iwp) :: off !< number of horizontally surrounding grid points to be sampled
	898	INTEGER(iwp) :: off_z !< number of vertically surrounding grid points to be sampled
[4498]	899	INTEGER(iwp) :: t !< running index over number of trajectories
[4400]	900
[4498]	901	INTEGER(KIND=1) :: soil_dum !< dummy variable to input a soil flag
[4400]	902
[4498]	903	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates
[4400]	904
[4536]	905	#if defined( __netcdf4_parallel )
[4498]	906	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ns_atmos !< number of observation points for each station on each mpi rank
	907	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: ns_soil !< number of observation points for each station on each mpi rank
[4444]	908	#endif
[4400]	909
[4498]	910	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions
[4400]	911
[4498]	912	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
[4400]	913
[4706]	914	REAL(wp) :: fill_height !< _FillValue for height coordinate (for trajectories)
	915	REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM
	916	REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM
	917	REAL(wp) :: fill_zar !< _FillValue for zar coordinate
[4400]	918
[4498]	919	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable
	920	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm_tmp !< EUTM coordinate before rotation
[4706]	921	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: height !< observation height above ground (for trajectories)
[4498]	922	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable
	923	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm_tmp !< NUTM coordinate before rotation
	924	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: station_h !< station height above reference
	925	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zar !< observation height above reference
[4400]	926	#if defined( __netcdf )
[3434]	927	!
[4400]	928	!-- Open the input file.
[4422]	929	CALL open_read_file( TRIM( input_file_vm ) // TRIM( coupling_char ), pids_id )
[3434]	930	!
[4400]	931	!-- Obtain number of sites.
[4498]	932	CALL get_attribute( pids_id, char_numstations, vmea_general%nvm, global_attribute )
[4400]	933	!
[4498]	934	!-- Allocate data structure which encompasses all required information, such as grid points indicies,
	935	!-- absolute UTM coordinates, the measured quantities, etc. .
[3704]	936	ALLOCATE( vmea(1:vmea_general%nvm) )
[3434]	937	!
[4498]	938	!-- Allocate flag array. This dummy array is used to identify grid points where virtual measurements
	939	!-- should be taken. Please note, in order to include also the surrounding grid points of the
	940	!-- original coordinate, ghost points are required.
[4400]	941	ALLOCATE( meas_flag(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[3522]	942	meas_flag = 0
	943	!
[4400]	944	!-- Loop over all sites in the setup file.
[3704]	945	DO l = 1, vmea_general%nvm
[3434]	946	!
[4498]	947	!-- Determine suffix which contains the ID, ordered according to the number of measurements.
[3434]	948	IF( l < 10 ) THEN
	949	WRITE( dum, '(I1)') l
	950	ELSEIF( l < 100 ) THEN
	951	WRITE( dum, '(I2)') l
	952	ELSEIF( l < 1000 ) THEN
	953	WRITE( dum, '(I3)') l
	954	ELSEIF( l < 10000 ) THEN
	955	WRITE( dum, '(I4)') l
	956	ELSEIF( l < 100000 ) THEN
	957	WRITE( dum, '(I5)') l
	958	ENDIF
[3704]	959	!
[4400]	960	!-- Read the origin site coordinates (UTM).
[4498]	961	CALL get_attribute( pids_id, char_origx // TRIM( dum ), vmea(l)%origin_x_obs, global_attribute )
	962	CALL get_attribute( pids_id, char_origy // TRIM( dum ), vmea(l)%origin_y_obs, global_attribute )
[3704]	963	!
[4400]	964	!-- Read site name.
[4498]	965	CALL get_attribute( pids_id, char_site // TRIM( dum ), vmea(l)%site, global_attribute )
[3704]	966	!
[4498]	967	!-- Read a flag which indicates that also soil quantities are take at the respective site
	968	!-- (is part of the virtual measurement driver).
	969	CALL get_attribute( pids_id, char_soil // TRIM( dum ), soil_dum, global_attribute )
[3704]	970	!
[4400]	971	!-- Set flag indicating soil-sampling.
	972	IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE.
[3704]	973	!
[4400]	974	!-- Read type of the measurement (trajectory, profile, timeseries).
[4498]	975	CALL get_attribute( pids_id, char_feature // TRIM( dum ), vmea(l)%feature_type, global_attribute )
[3434]	976	!
	977	!--- Set logicals depending on the type of the measurement
	978	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
	979	vmea(l)%timseries_profile = .TRUE.
	980	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
	981	vmea(l)%timseries = .TRUE.
	982	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
	983	vmea(l)%trajectory = .TRUE.
[3704]	984	!
[4400]	985	!-- Give error message in case the type matches non of the pre-defined types.
[3434]	986	ELSE
[4498]	987	message_string = 'Attribue featureType = ' // TRIM( vmea(l)%feature_type ) // ' is not allowed.'
[3717]	988	CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
[3434]	989	ENDIF
	990	!
[4400]	991	!-- Read string with all measured variables at this site.
[3434]	992	measured_variables_file = ''
[4498]	993	CALL get_variable( pids_id, char_mv // TRIM( dum ), measured_variables_file )
[3434]	994	!
[4400]	995	!-- Count the number of measured variables.
[4498]	996	!-- Please note, for some NetCDF interal reasons, characters end with a NULL, i.e. also empty
	997	!-- characters contain a NULL. Therefore, check the strings for a NULL to get the correct
	998	!-- character length in order to compare them with the list of allowed variables.
[4400]	999	vmea(l)%nmeas = 1
[4498]	1000	DO ll = 1, SIZE( measured_variables_file )
	1001	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
[3434]	1002	measured_variables_file(ll)(1:1) /= ' ') THEN
	1003	!
	1004	!-- Obtain character length of the character
	1005	len_char = 1
[4498]	1006	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0) .AND. &
	1007	measured_variables_file(ll)(len_char:len_char) /= ' ' )
[3434]	1008	len_char = len_char + 1
	1009	ENDDO
	1010	len_char = len_char - 1
[4400]	1011
[4498]	1012	measured_variables(vmea(l)%nmeas) = measured_variables_file(ll)(1:len_char)
[4400]	1013	vmea(l)%nmeas = vmea(l)%nmeas + 1
	1014
[3434]	1015	ENDIF
	1016	ENDDO
[4400]	1017	vmea(l)%nmeas = vmea(l)%nmeas - 1
[3434]	1018	!
[4498]	1019	!-- Allocate data-type array for the measured variables names and attributes at the respective
	1020	!-- site.
[4400]	1021	ALLOCATE( vmea(l)%var_atts(1:vmea(l)%nmeas) )
	1022	!
[4498]	1023	!-- Store the variable names in a data structure, which assigns further attributes to this name.
	1024	!-- Further, for data output reasons, create a string of output variables, which will be written
	1025	!-- into the attribute data_content.
[4400]	1026	DO ll = 1, vmea(l)%nmeas
	1027	vmea(l)%var_atts(ll)%name = TRIM( measured_variables(ll) )
[3434]	1028
[4641]	1029	! vmea(l)%data_content = TRIM( vmea(l)%data_content ) // " " // &
	1030	! TRIM( vmea(l)%var_atts(ll)%name )
[3434]	1031	ENDDO
	1032	!
[4498]	1033	!-- Read all the UTM coordinates for the site. Based on the coordinates, define the grid-index
	1034	!-- space on each subdomain where virtual measurements should be taken. Note, the entire
	1035	!-- coordinate array (on the entire model domain) won't be stored as this would exceed memory
	1036	!-- requirements, particularly for trajectories.
[3833]	1037	IF ( vmea(l)%nmeas > 0 ) THEN
[3434]	1038	!
[4498]	1039	!-- For stationary measurements UTM coordinates are just one value and its dimension is
	1040	!-- "station", while for mobile measurements UTM coordinates are arrays depending on the
	1041	!-- number of trajectories and time, according to (UC)2 standard. First, inquire dimension
	1042	!-- length of the UTM coordinates.
[3434]	1043	IF ( vmea(l)%trajectory ) THEN
	1044	!
[4498]	1045	!-- For non-stationary measurements read the number of trajectories and the number of time
	1046	!-- coordinates.
	1047	CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "traj" // TRIM( dum ) )
	1048	CALL get_dimension_length( pids_id, dim_ntime, "ntime" // TRIM( dum ) )
[3434]	1049	!
[4498]	1050	!-- For stationary measurements the dimension for UTM is station and for the time-coordinate
	1051	!-- it is one.
[3434]	1052	ELSE
[4498]	1053	CALL get_dimension_length( pids_id, vmea(l)%n_tr_st, "station" // TRIM( dum ) )
[3434]	1054	dim_ntime = 1
	1055	ENDIF
	1056	!
[4498]	1057	!- Allocate array which defines individual time/space frame for each trajectory or station.
[4400]	1058	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%n_tr_st) )
[3434]	1059	!
[4498]	1060	!-- Allocate temporary arrays for UTM and height coordinates. Note, on file UTM coordinates
	1061	!-- might be 1D or 2D variables
[4400]	1062	ALLOCATE( e_utm(1:vmea(l)%n_tr_st,1:dim_ntime) )
	1063	ALLOCATE( n_utm(1:vmea(l)%n_tr_st,1:dim_ntime) )
	1064	ALLOCATE( station_h(1:vmea(l)%n_tr_st,1:dim_ntime) )
	1065	ALLOCATE( zar(1:vmea(l)%n_tr_st,1:dim_ntime) )
[4706]	1066	IF ( vmea(l)%trajectory ) ALLOCATE( height(1:vmea(l)%n_tr_st,1:dim_ntime) )
[4400]	1067	e_utm = 0.0_wp
	1068	n_utm = 0.0_wp
	1069	station_h = 0.0_wp
	1070	zar = 0.0_wp
[4706]	1071	IF ( vmea(l)%trajectory ) height = 0.0_wp
[4400]	1072
	1073	ALLOCATE( e_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
	1074	ALLOCATE( n_utm_tmp(1:vmea(l)%n_tr_st,1:dim_ntime) )
[3434]	1075	!
[4498]	1076	!-- Read UTM and height coordinates for all trajectories and times. Note, in case
	1077	!-- these obtain any missing values, replace them with default _FillValues.
[4706]	1078	CALL inquire_fill_value( pids_id, char_eutm // TRIM( dum ), nofill, fill_eutm )
	1079	CALL inquire_fill_value( pids_id, char_nutm // TRIM( dum ), nofill, fill_nutm )
	1080	CALL inquire_fill_value( pids_id, char_zar // TRIM( dum ), nofill, fill_zar )
	1081	IF ( vmea(l)%trajectory ) &
	1082	CALL inquire_fill_value( pids_id, char_height // TRIM( dum ), nofill, fill_height )
[3434]	1083	!
[4498]	1084	!-- Further line is just to avoid compiler warnings. nofill might be used in future.
[4400]	1085	IF ( nofill == 0 .OR. nofill /= 0 ) CONTINUE
	1086	!
[4498]	1087	!-- Read observation coordinates. Please note, for trajectories the observation height is
	1088	!-- stored directly in z, while for timeSeries it is stored in z - station_h, according to
	1089	!-- UC2-standard.
[3437]	1090	IF ( vmea(l)%trajectory ) THEN
[4706]	1091	CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm, 0, dim_ntime-1, 0, &
[4498]	1092	vmea(l)%n_tr_st-1 )
[4706]	1093	CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm, 0, dim_ntime-1, 0, &
[4498]	1094	vmea(l)%n_tr_st-1 )
[4706]	1095	CALL get_variable( pids_id, char_zar // TRIM( dum ), zar, 0, dim_ntime-1, 0, &
[4498]	1096	vmea(l)%n_tr_st-1 )
[4706]	1097	CALL get_variable( pids_id, char_height // TRIM( dum ), height, 0, dim_ntime-1, 0, &
	1098	vmea(l)%n_tr_st-1 )
[3437]	1099	ELSE
[4706]	1100	CALL get_variable( pids_id, char_eutm // TRIM( dum ), e_utm(:,1) )
	1101	CALL get_variable( pids_id, char_nutm // TRIM( dum ), n_utm(:,1) )
[4498]	1102	CALL get_variable( pids_id, char_station_h // TRIM( dum ), station_h(:,1) )
[4706]	1103	CALL get_variable( pids_id, char_zar // TRIM( dum ), zar(:,1) )
[4400]	1104	ENDIF
	1105
	1106	e_utm = MERGE( e_utm, vmea(l)%fillout, e_utm /= fill_eutm )
	1107	n_utm = MERGE( n_utm, vmea(l)%fillout, n_utm /= fill_nutm )
	1108	zar = MERGE( zar, vmea(l)%fillout, zar /= fill_zar )
[4706]	1109	IF ( vmea(l)%trajectory ) &
	1110	height = MERGE( height, vmea(l)%fillout, height /= fill_height )
[3434]	1111	!
[4706]	1112	!-- Compute observation height above ground. Note, for trajectory measurements the height
	1113	!-- above the surface is actually stored in 'height'
	1114	IF ( vmea(l)%trajectory ) THEN
	1115	zar = height
	1116	fill_zar = fill_height
	1117	ELSE
	1118	zar = zar - station_h
	1119	ENDIF
[4400]	1120	!
[4498]	1121	!-- Based on UTM coordinates, check if the measurement station or parts of the trajectory are
	1122	!-- on subdomain. This case, setup grid index space sample these quantities.
[3522]	1123	meas_flag = 0
[4400]	1124	DO t = 1, vmea(l)%n_tr_st
	1125	!
[4498]	1126	!-- First, compute relative x- and y-coordinates with respect to the lower-left origin of
	1127	!-- the model domain, which is the difference between UTM coordinates. Note, if the origin
	1128	!-- is not correct, the virtual sites will be misplaced. Further, in case of an rotated
	1129	!-- model domain, the UTM coordinates must also be rotated.
[3910]	1130	e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
	1131	n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
[4498]	1132	e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp ) &
	1133	* e_utm_tmp(t,1:dim_ntime) &
	1134	- SIN( init_model%rotation_angle * pi / 180.0_wp ) &
[3910]	1135	* n_utm_tmp(t,1:dim_ntime)
[4498]	1136	n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp ) &
	1137	* e_utm_tmp(t,1:dim_ntime) &
	1138	+ COS( init_model%rotation_angle * pi / 180.0_wp ) &
[3910]	1139	* n_utm_tmp(t,1:dim_ntime)
[3434]	1140	!
[4498]	1141	!-- Determine the individual time coordinate length for each station and trajectory. This
	1142	!-- is required as several stations and trajectories are merged into one file but they do
	1143	!-- not have the same number of points in time, hence, missing values may occur and cannot
	1144	!-- be processed further. This is actually a work-around for the specific (UC)2 dataset,
	1145	!-- but it won't harm anyway.
[3434]	1146	vmea(l)%dim_t(t) = 0
	1147	DO n = 1, dim_ntime
[4498]	1148	IF ( e_utm(t,n) /= fill_eutm .AND. n_utm(t,n) /= fill_nutm .AND. &
[4400]	1149	zar(t,n) /= fill_zar ) vmea(l)%dim_t(t) = n
[3434]	1150	ENDDO
	1151	!
[4498]	1152	!-- Compute grid indices relative to origin and check if these are on the subdomain. Note,
	1153	!-- virtual measurements will be taken also at grid points surrounding the station, hence,
	1154	!-- check also for these grid points. The number of surrounding grid points is set
	1155	!-- according to the featureType.
[4400]	1156	IF ( vmea(l)%timseries_profile ) THEN
[4695]	1157	off = off_pr
	1158	off_z = off_pr_z
[4400]	1159	ELSEIF ( vmea(l)%timseries ) THEN
[4695]	1160	off = off_ts
	1161	off_z = off_ts_z
[4400]	1162	ELSEIF ( vmea(l)%trajectory ) THEN
[4695]	1163	off = off_tr
	1164	off_z = off_tr_z
[4400]	1165	ENDIF
	1166
[3437]	1167	DO n = 1, vmea(l)%dim_t(t)
[4498]	1168	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
	1169	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
[3434]	1170	!
	1171	!-- Is the observation point on subdomain?
[4498]	1172	on_pe = ( is >= nxl .AND. is <= nxr .AND. js >= nys .AND. js <= nyn )
[3434]	1173	!
[4498]	1174	!-- Check if observation coordinate is on subdomain.
[3434]	1175	IF ( on_pe ) THEN
[3522]	1176	!
[4498]	1177	!-- Determine vertical index which corresponds to the observation height.
[4168]	1178	ksurf = topo_top_ind(js,is,0)
[4400]	1179	ks = MINLOC( ABS( zu - zw(ksurf) - zar(t,n) ), DIM = 1 ) - 1
[3434]	1180	!
[4498]	1181	!-- Set mask array at the observation coordinates. Also, flag the surrounding
	1182	!-- coordinate points, but first check whether the surrounding coordinate points are
	1183	!-- on the subdomain.
[4695]	1184	kl = MERGE( ks-off_z, ksurf, ks-off_z >= nzb .AND. ks-off_z >= ksurf )
	1185	ku = MERGE( ks+off_z, nzt, ks+off_z < nzt+1 )
[4400]	1186
	1187	DO i = is-off, is+off
	1188	DO j = js-off, js+off
[3704]	1189	DO k = kl, ku
[4498]	1190	meas_flag(k,j,i) = MERGE( IBSET( meas_flag(k,j,i), 0 ), 0, &
[4645]	1191	BTEST( wall_flags_total_0(k,j,i), 0 ) )
[3704]	1192	ENDDO
	1193	ENDDO
	1194	ENDDO
[3434]	1195	ENDIF
	1196	ENDDO
[4400]	1197
[3434]	1198	ENDDO
	1199	!
[4498]	1200	!-- Based on the flag array, count the number of sampling coordinates. Please note, sampling
	1201	!-- coordinates in atmosphere and soil may be different, as within the soil all levels will be
	1202	!-- measured. Hence, count individually. Start with atmoshere.
[3522]	1203	ns = 0
[4400]	1204	DO i = nxl-off, nxr+off
	1205	DO j = nys-off, nyn+off
[3704]	1206	DO k = nzb, nzt+1
	1207	ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
[3522]	1208	ENDDO
	1209	ENDDO
	1210	ENDDO
[4400]	1211
[3522]	1212	!
[4498]	1213	!-- Store number of observation points on subdomain and allocate index arrays as well as array
	1214	!-- containing height information.
[3434]	1215	vmea(l)%ns = ns
[4400]	1216
[3434]	1217	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
	1218	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
	1219	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
[4400]	1220	ALLOCATE( vmea(l)%zar(1:vmea(l)%ns) )
[3434]	1221	!
[4498]	1222	!-- Based on the flag array store the grid indices which correspond to the observation
	1223	!-- coordinates.
[3704]	1224	ns = 0
[4400]	1225	DO i = nxl-off, nxr+off
	1226	DO j = nys-off, nyn+off
[3704]	1227	DO k = nzb, nzt+1
	1228	IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN
[3522]	1229	ns = ns + 1
[3704]	1230	vmea(l)%i(ns) = i
	1231	vmea(l)%j(ns) = j
	1232	vmea(l)%k(ns) = k
[4400]	1233	vmea(l)%zar(ns) = zu(k) - zw(topo_top_ind(j,i,0))
[3522]	1234	ENDIF
	1235	ENDDO
[3434]	1236	ENDDO
	1237	ENDDO
	1238	!
[4498]	1239	!-- Same for the soil. Based on the flag array, count the number of sampling coordinates in
	1240	!-- soil. Sample at all soil levels in this case. Please note, soil variables can only be
	1241	!-- sampled on subdomains, not on ghost layers.
[3704]	1242	IF ( vmea(l)%soil_sampling ) THEN
	1243	DO i = nxl, nxr
	1244	DO j = nys, nyn
	1245	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
[4671]	1246	IF ( surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i) ) THEN
[4498]	1247	vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_soil - nzb_soil + 1
[3704]	1248	ENDIF
[4671]	1249	IF ( surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i) ) THEN
[4498]	1250	vmea(l)%ns_soil = vmea(l)%ns_soil + nzt_wall - nzb_wall + 1
[3704]	1251	ENDIF
	1252	ENDIF
	1253	ENDDO
	1254	ENDDO
[4400]	1255	ENDIF
[3704]	1256	!
[4498]	1257	!-- Allocate index arrays as well as array containing height information for soil.
[3704]	1258	IF ( vmea(l)%soil_sampling ) THEN
	1259	ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
	1260	ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
	1261	ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
[4400]	1262	ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
	1263	ENDIF
[3704]	1264	!
	1265	!-- For soil, store the grid indices.
	1266	ns = 0
	1267	IF ( vmea(l)%soil_sampling ) THEN
	1268	DO i = nxl, nxr
	1269	DO j = nys, nyn
	1270	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
[4671]	1271	IF ( surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i) ) THEN
	1272	m = surf_lsm_h(0)%start_index(j,i)
[3704]	1273	DO k = nzb_soil, nzt_soil
	1274	ns = ns + 1
	1275	vmea(l)%i_soil(ns) = i
	1276	vmea(l)%j_soil(ns) = j
	1277	vmea(l)%k_soil(ns) = k
[4400]	1278	vmea(l)%depth(ns) = - zs(k)
[3704]	1279	ENDDO
	1280	ENDIF
[4400]	1281
[4671]	1282	IF ( surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i) ) THEN
	1283	m = surf_usm_h(0)%start_index(j,i)
[3704]	1284	DO k = nzb_wall, nzt_wall
	1285	ns = ns + 1
	1286	vmea(l)%i_soil(ns) = i
	1287	vmea(l)%j_soil(ns) = j
	1288	vmea(l)%k_soil(ns) = k
[4671]	1289	vmea(l)%depth(ns) = - surf_usm_h(0)%zw(k,m)
[3704]	1290	ENDDO
	1291	ENDIF
	1292	ENDIF
	1293	ENDDO
	1294	ENDDO
	1295	ENDIF
	1296	!
[3434]	1297	!-- Allocate array to save the sampled values.
[3833]	1298	ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
[4400]	1299
[4498]	1300	IF ( vmea(l)%soil_sampling ) &
	1301	ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, 1:vmea(l)%nmeas) )
[3434]	1302	!
[3704]	1303	!-- Initialize with _FillValues
[3833]	1304	vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout
[4498]	1305	IF ( vmea(l)%soil_sampling ) &
	1306	vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = vmea(l)%fillout
[3434]	1307	!
	1308	!-- Deallocate temporary coordinate arrays
[3910]	1309	IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm )
	1310	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
	1311	IF ( ALLOCATED( e_utm_tmp ) ) DEALLOCATE( e_utm_tmp )
	1312	IF ( ALLOCATED( n_utm_tmp ) ) DEALLOCATE( n_utm_tmp )
	1313	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
[4400]	1314	IF ( ALLOCATED( zar ) ) DEALLOCATE( vmea(l)%dim_t )
	1315	IF ( ALLOCATED( zar ) ) DEALLOCATE( zar )
[4706]	1316	IF ( ALLOCATED( height ) ) DEALLOCATE( height )
[4400]	1317	IF ( ALLOCATED( station_h ) ) DEALLOCATE( station_h )
	1318
[3434]	1319	ENDIF
	1320	ENDDO
	1321	!
[4400]	1322	!-- Dellocate flag array
	1323	DEALLOCATE( meas_flag )
[3704]	1324	!
[4408]	1325	!-- Close input file for virtual measurements.
[4400]	1326	CALL close_input_file( pids_id )
	1327	!
	1328	!-- Sum-up the number of observation coordiates, for atmosphere first.
[3704]	1329	!-- This is actually only required for data output.
	1330	ALLOCATE( ns_all(1:vmea_general%nvm) )
[4400]	1331	ns_all = 0
[3704]	1332	#if defined( __parallel )
[4498]	1333	CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, &
	1334	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[3704]	1335	#else
	1336	ns_all(:) = vmea(:)%ns
	1337	#endif
	1338	vmea(:)%ns_tot = ns_all(:)
	1339	!
	1340	!-- Now for soil
[4400]	1341	ns_all = 0
[3704]	1342	#if defined( __parallel )
[4498]	1343	CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, &
[3704]	1344	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	1345	#else
	1346	ns_all(:) = vmea(:)%ns_soil
	1347	#endif
	1348	vmea(:)%ns_soil_tot = ns_all(:)
[4400]	1349
[3704]	1350	DEALLOCATE( ns_all )
	1351	!
[4498]	1352	!-- In case of parallel NetCDF the start coordinate for each mpi rank needs to be defined, so that
	1353	!-- each processor knows where to write the data.
[4400]	1354	#if defined( __netcdf4_parallel )
	1355	ALLOCATE( ns_atmos(0:numprocs-1,1:vmea_general%nvm) )
	1356	ALLOCATE( ns_soil(0:numprocs-1,1:vmea_general%nvm) )
	1357	ns_atmos = 0
	1358	ns_soil = 0
	1359
	1360	DO l = 1, vmea_general%nvm
	1361	ns_atmos(myid,l) = vmea(l)%ns
	1362	ns_soil(myid,l) = vmea(l)%ns_soil
	1363	ENDDO
	1364
	1365	#if defined( __parallel )
[4498]	1366	CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_atmos, numprocs * vmea_general%nvm, &
[4400]	1367	MPI_INTEGER, MPI_SUM, comm2d, ierr )
[4498]	1368	CALL MPI_ALLREDUCE( MPI_IN_PLACE, ns_soil, numprocs * vmea_general%nvm, &
[4400]	1369	MPI_INTEGER, MPI_SUM, comm2d, ierr )
	1370	#else
	1371	ns_atmos(0,:) = vmea(:)%ns
	1372	ns_soil(0,:) = vmea(:)%ns_soil
	1373	#endif
	1374
[3704]	1375	!
[4498]	1376	!-- Determine the start coordinate in NetCDF file for the local arrays. Note, start coordinates are
	1377	!-- initialized with zero for sake of simplicity in summation. However, in NetCDF the start
	1378	!-- coordinates must be >= 1, so that a one needs to be added at the end.
[4400]	1379	DO l = 1, vmea_general%nvm
	1380	DO n = 0, myid - 1
	1381	vmea(l)%start_coord_a = vmea(l)%start_coord_a + ns_atmos(n,l)
	1382	vmea(l)%start_coord_s = vmea(l)%start_coord_s + ns_soil(n,l)
	1383	ENDDO
	1384	!
	1385	!-- Start coordinate in NetCDF starts always at one not at 0.
	1386	vmea(l)%start_coord_a = vmea(l)%start_coord_a + 1
	1387	vmea(l)%start_coord_s = vmea(l)%start_coord_s + 1
	1388	!
	1389	!-- Determine the local end coordinate
	1390	vmea(l)%end_coord_a = vmea(l)%start_coord_a + vmea(l)%ns - 1
	1391	vmea(l)%end_coord_s = vmea(l)%start_coord_s + vmea(l)%ns_soil - 1
	1392	ENDDO
	1393
	1394	DEALLOCATE( ns_atmos )
	1395	DEALLOCATE( ns_soil )
	1396
	1397	#endif
	1398
	1399	#endif
	1400
[4498]	1401	END SUBROUTINE vm_init
[4400]	1402
	1403
[4498]	1404	!--------------------------------------------------------------------------------------------------!
[3434]	1405	! Description:
	1406	! ------------
[4400]	1407	!> Initialize output using data-output module
[4498]	1408	!--------------------------------------------------------------------------------------------------!
[4400]	1409	SUBROUTINE vm_init_output
	1410
	1411	CHARACTER(LEN=100) :: variable_name !< name of output variable
	1412
[4498]	1413	INTEGER(iwp) :: l !< loop index
	1414	INTEGER(iwp) :: n !< loop index
	1415	INTEGER :: return_value !< returned status value of called function
[4400]	1416
	1417	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ndim !< dummy to write dimension
	1418
[4498]	1419	REAL(wp) :: dum_lat !< transformed geographical coordinate (latitude)
	1420	REAL(wp) :: dum_lon !< transformed geographical coordinate (longitude)
[4400]	1421
[3704]	1422	!
[4438]	1423	!-- Determine the number of output timesteps.
[4504]	1424	ntimesteps = CEILING( ( end_time - MAX( vm_time_start, time_since_reference_point ) &
	1425	) / dt_virtual_measurement )
[4400]	1426	!
	1427	!-- Create directory where output files will be stored.
	1428	CALL local_system( 'mkdir -p VM_OUTPUT' // TRIM( coupling_char ) )
	1429	!
	1430	!-- Loop over all sites.
	1431	DO l = 1, vmea_general%nvm
	1432	!
	1433	!-- Skip if no observations will be taken for this site.
	1434	IF ( vmea(l)%ns_tot == 0 .AND. vmea(l)%ns_soil_tot == 0 ) CYCLE
	1435	!
	1436	!-- Define output file.
[4498]	1437	WRITE( vmea(l)%nc_filename, '(A,I4.4)' ) 'VM_OUTPUT' // TRIM( coupling_char ) // '/' // &
	1438	'site', l
[3704]	1439
[4400]	1440	return_value = dom_def_file( vmea(l)%nc_filename, 'netcdf4-parallel' )
	1441	!
	1442	!-- Define global attributes.
	1443	!-- Before, transform UTM into geographical coordinates.
[4498]	1444	CALL convert_utm_to_geographic( crs_list, vmea(l)%origin_x_obs, vmea(l)%origin_y_obs, &
	1445	dum_lon, dum_lat )
[4400]	1446
[4504]	1447	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'site', &
[4400]	1448	value = TRIM( vmea(l)%site ) )
[4504]	1449	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'title', &
[4400]	1450	value = 'Virtual measurement output')
[4504]	1451	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'source', &
[4400]	1452	value = 'PALM-4U')
[4504]	1453	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'institution', &
[4400]	1454	value = input_file_atts%institution )
[4504]	1455	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'acronym', &
[4400]	1456	value = input_file_atts%acronym )
[4504]	1457	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'author', &
[4400]	1458	value = input_file_atts%author )
[4504]	1459	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'contact_person', &
[4641]	1460	value = input_file_atts%author )
[4504]	1461	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'iop', &
[4400]	1462	value = input_file_atts%campaign )
[4504]	1463	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'campaign', &
[4400]	1464	value = 'PALM-4U' )
[4504]	1465	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_time ', &
[4400]	1466	value = origin_date_time)
[4504]	1467	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'location', &
[4400]	1468	value = input_file_atts%location )
[4504]	1469	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_x', &
[4400]	1470	value = vmea(l)%origin_x_obs )
[4504]	1471	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_y', &
[4400]	1472	value = vmea(l)%origin_y_obs )
[4504]	1473	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_lon', &
[4400]	1474	value = dum_lon )
[4504]	1475	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_lat', &
[4400]	1476	value = dum_lat )
[4504]	1477	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'origin_z', value = 0.0 )
	1478	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'rotation_angle', &
[4400]	1479	value = input_file_atts%rotation_angle )
[4504]	1480	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'featureType', &
[4400]	1481	value = TRIM( vmea(l)%feature_type_out ) )
[4504]	1482	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'data_content', &
[4400]	1483	value = TRIM( vmea(l)%data_content ) )
[4504]	1484	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'creation_time', &
[4400]	1485	value = input_file_atts%creation_time )
[4504]	1486	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'version', value = 1 ) !input_file_atts%version
	1487	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'Conventions', &
[4400]	1488	value = input_file_atts%conventions )
[4504]	1489	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'dependencies', &
[4400]	1490	value = input_file_atts%dependencies )
[4504]	1491	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'history', &
[4400]	1492	value = input_file_atts%history )
[4504]	1493	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'references', &
[4400]	1494	value = input_file_atts%references )
[4504]	1495	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'comment', &
[4400]	1496	value = input_file_atts%comment )
[4504]	1497	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'keywords', &
[4400]	1498	value = input_file_atts%keywords )
[4504]	1499	return_value = dom_def_att( vmea(l)%nc_filename, attribute_name = 'licence', &
[4498]	1500	value = '[UC]2 Open Licence; see [UC]2 ' // &
	1501	'data policy available at ' // &
[4400]	1502	'www.uc2-program.org/uc2_data_policy.pdf' )
	1503	!
	1504	!-- Define dimensions.
	1505	!-- station
	1506	ALLOCATE( ndim(1:vmea(l)%ns_tot) )
	1507	DO n = 1, vmea(l)%ns_tot
	1508	ndim(n) = n
	1509	ENDDO
[4498]	1510	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'station', &
	1511	output_type = 'int32', bounds = (/1_iwp, vmea(l)%ns_tot/), &
[4400]	1512	values_int32 = ndim )
	1513	DEALLOCATE( ndim )
	1514	!
	1515	!-- ntime
	1516	ALLOCATE( ndim(1:ntimesteps) )
	1517	DO n = 1, ntimesteps
	1518	ndim(n) = n
	1519	ENDDO
	1520
[4498]	1521	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'ntime', &
	1522	output_type = 'int32', bounds = (/1_iwp, ntimesteps/), &
[4400]	1523	values_int32 = ndim )
	1524	DEALLOCATE( ndim )
	1525	!
	1526	!-- nv
	1527	ALLOCATE( ndim(1:2) )
	1528	DO n = 1, 2
	1529	ndim(n) = n
	1530	ENDDO
	1531
[4498]	1532	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'nv', &
	1533	output_type = 'int32', bounds = (/1_iwp, 2_iwp/), &
[4400]	1534	values_int32 = ndim )
	1535	DEALLOCATE( ndim )
	1536	!
	1537	!-- maximum name length
	1538	ALLOCATE( ndim(1:maximum_name_length) )
	1539	DO n = 1, maximum_name_length
	1540	ndim(n) = n
	1541	ENDDO
	1542
[4498]	1543	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'max_name_len', &
	1544	output_type = 'int32', &
	1545	bounds = (/1_iwp, maximum_name_length /), values_int32 = ndim )
[4400]	1546	DEALLOCATE( ndim )
	1547	!
	1548	!-- Define coordinate variables.
	1549	!-- time
	1550	variable_name = 'time'
[4498]	1551	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1552	dimension_names = (/ 'station ', 'ntime '/), &
[4400]	1553	output_type = 'real32' )
	1554	!
[4421]	1555	!-- station_name
	1556	variable_name = 'station_name'
[4498]	1557	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1558	dimension_names = (/ 'max_name_len', 'station ' /), &
[4421]	1559	output_type = 'char' )
[4400]	1560	!
	1561	!-- vrs (vertical reference system)
	1562	variable_name = 'vrs'
[4498]	1563	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1564	dimension_names = (/ 'station' /), output_type = 'int8' )
[4400]	1565	!
	1566	!-- crs (coordinate reference system)
	1567	variable_name = 'crs'
[4498]	1568	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1569	dimension_names = (/ 'station' /), output_type = 'int8' )
[4400]	1570	!
	1571	!-- z
	1572	variable_name = 'z'
[4498]	1573	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1574	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1575	!
	1576	!-- station_h
	1577	variable_name = 'station_h'
[4498]	1578	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1579	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1580	!
	1581	!-- x
	1582	variable_name = 'x'
[4498]	1583	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1584	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1585	!
	1586	!-- y
	1587	variable_name = 'y'
[4498]	1588	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1589	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1590	!
	1591	!-- E-UTM
	1592	variable_name = 'E_UTM'
[4498]	1593	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1594	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1595	!
	1596	!-- N-UTM
	1597	variable_name = 'N_UTM'
[4498]	1598	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1599	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1600	!
	1601	!-- latitude
	1602	variable_name = 'lat'
[4498]	1603	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1604	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1605	!
	1606	!-- longitude
	1607	variable_name = 'lon'
[4498]	1608	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1609	dimension_names = (/'station'/), output_type = 'real32' )
[4400]	1610	!
[4498]	1611	!-- Set attributes for the coordinate variables. Note, not all coordinates have the same number
	1612	!-- of attributes.
[4400]	1613	!-- Units
[4504]	1614	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1615	attribute_name = char_unit, value = 'seconds since ' // &
	1616	origin_date_time )
	1617	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z', &
	1618	attribute_name = char_unit, value = 'm' )
	1619	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h', &
	1620	attribute_name = char_unit, value = 'm' )
	1621	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x', &
	1622	attribute_name = char_unit, value = 'm' )
	1623	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y', &
	1624	attribute_name = char_unit, value = 'm' )
	1625	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM', &
	1626	attribute_name = char_unit, value = 'm' )
	1627	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM', &
	1628	attribute_name = char_unit, value = 'm' )
	1629	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat', &
	1630	attribute_name = char_unit, value = 'degrees_north' )
	1631	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon', &
	1632	attribute_name = char_unit, value = 'degrees_east' )
[4400]	1633	!
	1634	!-- Long name
[4504]	1635	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name', &
	1636	attribute_name = char_long, value = 'station name')
	1637	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1638	attribute_name = char_long, value = 'time')
	1639	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z', &
	1640	attribute_name = char_long, value = 'height above origin' )
	1641	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h', &
	1642	attribute_name = char_long, value = 'surface altitude' )
	1643	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x', &
[4498]	1644	attribute_name = char_long, &
	1645	value = 'distance to origin in x-direction')
[4504]	1646	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y', &
[4498]	1647	attribute_name = char_long, &
	1648	value = 'distance to origin in y-direction')
[4504]	1649	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM', &
	1650	attribute_name = char_long, value = 'easting' )
	1651	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM', &
	1652	attribute_name = char_long, value = 'northing' )
	1653	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat', &
	1654	attribute_name = char_long, value = 'latitude' )
	1655	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon', &
	1656	attribute_name = char_long, value = 'longitude' )
[4400]	1657	!
	1658	!-- Standard name
[4504]	1659	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name', &
	1660	attribute_name = char_standard, value = 'platform_name')
	1661	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1662	attribute_name = char_standard, value = 'time')
	1663	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z', &
[4498]	1664	attribute_name = char_standard, &
[4400]	1665	value = 'height_above_mean_sea_level' )
[4504]	1666	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h', &
	1667	attribute_name = char_standard, value = 'surface_altitude' )
	1668	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM', &
[4498]	1669	attribute_name = char_standard, &
[4400]	1670	value = 'projection_x_coordinate' )
[4504]	1671	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM', &
[4498]	1672	attribute_name = char_standard, &
[4400]	1673	value = 'projection_y_coordinate' )
[4504]	1674	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat', &
	1675	attribute_name = char_standard, value = 'latitude' )
	1676	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon', &
	1677	attribute_name = char_standard, value = 'longitude' )
[4400]	1678	!
	1679	!-- Axis
[4504]	1680	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1681	attribute_name = 'axis', value = 'T')
	1682	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z', &
	1683	attribute_name = 'axis', value = 'Z' )
	1684	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x', &
	1685	attribute_name = 'axis', value = 'X' )
	1686	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y', &
	1687	attribute_name = 'axis', value = 'Y' )
[4400]	1688	!
	1689	!-- Set further individual attributes for the coordinate variables.
	1690	!-- For station name
[4504]	1691	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name', &
	1692	attribute_name = 'cf_role', value = 'timeseries_id' )
[4400]	1693	!
	1694	!-- For time
[4504]	1695	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1696	attribute_name = 'calendar', value = 'proleptic_gregorian' )
[4642]	1697	! return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time', &
	1698	! attribute_name = 'bounds', value = 'time_bounds' )
[4400]	1699	!
	1700	!-- For vertical reference system
[4504]	1701	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'vrs', &
	1702	attribute_name = char_long, value = 'vertical reference system' )
	1703	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'vrs', &
	1704	attribute_name = 'system_name', value = 'DHHN2016' )
[4400]	1705	!
	1706	!-- For z
[4504]	1707	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z', &
	1708	attribute_name = 'positive', value = 'up' )
[4400]	1709	!
	1710	!-- For coordinate reference system
[4504]	1711	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
	1712	attribute_name = 'epsg_code', value = coord_ref_sys%epsg_code )
	1713	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1714	attribute_name = 'false_easting', &
[4400]	1715	value = coord_ref_sys%false_easting )
[4504]	1716	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1717	attribute_name = 'false_northing', &
[4400]	1718	value = coord_ref_sys%false_northing )
[4504]	1719	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1720	attribute_name = 'grid_mapping_name', &
[4400]	1721	value = coord_ref_sys%grid_mapping_name )
[4504]	1722	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1723	attribute_name = 'inverse_flattening', &
[4400]	1724	value = coord_ref_sys%inverse_flattening )
[4504]	1725	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4400]	1726	attribute_name = 'latitude_of_projection_origin',&
	1727	value = coord_ref_sys%latitude_of_projection_origin )
[4504]	1728	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
	1729	attribute_name = char_long, value = coord_ref_sys%long_name )
	1730	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1731	attribute_name = 'longitude_of_central_meridian', &
[4400]	1732	value = coord_ref_sys%longitude_of_central_meridian )
[4504]	1733	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1734	attribute_name = 'longitude_of_prime_meridian', &
[4400]	1735	value = coord_ref_sys%longitude_of_prime_meridian )
[4504]	1736	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1737	attribute_name = 'scale_factor_at_central_meridian', &
[4400]	1738	value = coord_ref_sys%scale_factor_at_central_meridian )
[4504]	1739	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
[4498]	1740	attribute_name = 'semi_major_axis', &
[4400]	1741	value = coord_ref_sys%semi_major_axis )
[4504]	1742	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'crs', &
	1743	attribute_name = char_unit, value = coord_ref_sys%units )
[4400]	1744	!
[4498]	1745	!-- In case of sampled soil quantities, define further dimensions and coordinates.
[4400]	1746	IF ( vmea(l)%soil_sampling ) THEN
	1747	!
	1748	!-- station for soil
	1749	ALLOCATE( ndim(1:vmea(l)%ns_soil_tot) )
	1750	DO n = 1, vmea(l)%ns_soil_tot
	1751	ndim(n) = n
	1752	ENDDO
	1753
[4504]	1754	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'station_soil', &
[4498]	1755	output_type = 'int32', &
[4504]	1756	bounds = (/1_iwp,vmea(l)%ns_soil_tot/), values_int32 = ndim )
[4400]	1757	DEALLOCATE( ndim )
	1758	!
	1759	!-- ntime for soil
	1760	ALLOCATE( ndim(1:ntimesteps) )
	1761	DO n = 1, ntimesteps
	1762	ndim(n) = n
	1763	ENDDO
	1764
[4504]	1765	return_value = dom_def_dim( vmea(l)%nc_filename, dimension_name = 'ntime_soil', &
	1766	output_type = 'int32', bounds = (/1_iwp,ntimesteps/), &
[4400]	1767	values_int32 = ndim )
	1768	DEALLOCATE( ndim )
	1769	!
	1770	!-- time for soil
	1771	variable_name = 'time_soil'
[4504]	1772	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1773	dimension_names = (/'station_soil', 'ntime_soil '/), &
[4400]	1774	output_type = 'real32' )
	1775	!
	1776	!-- station_name for soil
	1777	variable_name = 'station_name_soil'
[4504]	1778	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1779	dimension_names = (/ 'max_name_len', 'station_soil' /), &
[4400]	1780	output_type = 'char' )
	1781	!
	1782	!-- z
	1783	variable_name = 'z_soil'
[4504]	1784	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1785	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1786	!
	1787	!-- station_h for soil
	1788	variable_name = 'station_h_soil'
[4504]	1789	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1790	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1791	!
	1792	!-- x soil
	1793	variable_name = 'x_soil'
[4504]	1794	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1795	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1796	!
	1797	!- y soil
	1798	variable_name = 'y_soil'
[4504]	1799	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1800	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1801	!
	1802	!-- E-UTM soil
	1803	variable_name = 'E_UTM_soil'
[4504]	1804	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1805	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1806	!
	1807	!-- N-UTM soil
	1808	variable_name = 'N_UTM_soil'
[4504]	1809	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1810	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1811	!
	1812	!-- latitude soil
	1813	variable_name = 'lat_soil'
[4504]	1814	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1815	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1816	!
	1817	!-- longitude soil
	1818	variable_name = 'lon_soil'
[4504]	1819	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1820	dimension_names = (/'station_soil'/), output_type = 'real32' )
[4400]	1821	!
[4498]	1822	!-- Set attributes for the coordinate variables. Note, not all coordinates have the same
	1823	!-- number of attributes.
[4400]	1824	!-- Units
[4504]	1825	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1826	attribute_name = char_unit, value = 'seconds since ' // &
	1827	origin_date_time )
	1828	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil', &
	1829	attribute_name = char_unit, value = 'm' )
	1830	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil', &
	1831	attribute_name = char_unit, value = 'm' )
	1832	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil', &
	1833	attribute_name = char_unit, value = 'm' )
	1834	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil', &
	1835	attribute_name = char_unit, value = 'm' )
	1836	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil', &
	1837	attribute_name = char_unit, value = 'm' )
	1838	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil', &
	1839	attribute_name = char_unit, value = 'm' )
	1840	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil', &
	1841	attribute_name = char_unit, value = 'degrees_north' )
	1842	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil', &
	1843	attribute_name = char_unit, value = 'degrees_east' )
[4400]	1844	!
	1845	!-- Long name
[4504]	1846	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil', &
	1847	attribute_name = char_long, value = 'station name')
	1848	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1849	attribute_name = char_long, value = 'time')
	1850	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil', &
	1851	attribute_name = char_long, value = 'height above origin' )
	1852	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil', &
	1853	attribute_name = char_long, value = 'surface altitude' )
	1854	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil', &
[4498]	1855	attribute_name = char_long, &
[4400]	1856	value = 'distance to origin in x-direction' )
[4504]	1857	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil', &
[4498]	1858	attribute_name = char_long, &
[4400]	1859	value = 'distance to origin in y-direction' )
[4504]	1860	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil', &
	1861	attribute_name = char_long, value = 'easting' )
	1862	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil', &
	1863	attribute_name = char_long, value = 'northing' )
	1864	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil', &
	1865	attribute_name = char_long, value = 'latitude' )
	1866	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil', &
	1867	attribute_name = char_long, value = 'longitude' )
[4400]	1868	!
	1869	!-- Standard name
[4504]	1870	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil', &
	1871	attribute_name = char_standard, value = 'platform_name')
	1872	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1873	attribute_name = char_standard, value = 'time')
	1874	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil', &
[4498]	1875	attribute_name = char_standard, &
[4400]	1876	value = 'height_above_mean_sea_level' )
[4504]	1877	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_h_soil', &
	1878	attribute_name = char_standard, value = 'surface_altitude' )
	1879	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'E_UTM_soil', &
[4498]	1880	attribute_name = char_standard, &
[4400]	1881	value = 'projection_x_coordinate' )
[4504]	1882	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'N_UTM_soil', &
[4498]	1883	attribute_name = char_standard, &
[4400]	1884	value = 'projection_y_coordinate' )
[4504]	1885	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lat_soil', &
	1886	attribute_name = char_standard, value = 'latitude' )
	1887	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'lon_soil', &
	1888	attribute_name = char_standard, value = 'longitude' )
[4400]	1889	!
	1890	!-- Axis
[4504]	1891	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1892	attribute_name = 'axis', value = 'T')
	1893	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil', &
	1894	attribute_name = 'axis', value = 'Z' )
	1895	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'x_soil', &
	1896	attribute_name = 'axis', value = 'X' )
	1897	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'y_soil', &
	1898	attribute_name = 'axis', value = 'Y' )
[4400]	1899	!
	1900	!-- Set further individual attributes for the coordinate variables.
	1901	!-- For station name soil
[4504]	1902	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'station_name_soil', &
	1903	attribute_name = 'cf_role', value = 'timeseries_id' )
[4400]	1904	!
	1905	!-- For time soil
[4504]	1906	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1907	attribute_name = 'calendar', value = 'proleptic_gregorian' )
[4642]	1908	! return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'time_soil', &
	1909	! attribute_name = 'bounds', value = 'time_bounds' )
[4400]	1910	!
	1911	!-- For z soil
[4504]	1912	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = 'z_soil', &
	1913	attribute_name = 'positive', value = 'up' )
[4400]	1914	ENDIF
	1915	!
	1916	!-- Define variables that shall be sampled.
	1917	DO n = 1, vmea(l)%nmeas
	1918	variable_name = TRIM( vmea(l)%var_atts(n)%name )
	1919	!
[4498]	1920	!-- In order to link the correct dimension names, atmosphere and soil variables need to be
	1921	!-- distinguished.
	1922	IF ( vmea(l)%soil_sampling .AND. &
[4400]	1923	ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) ) THEN
	1924
[4504]	1925	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1926	dimension_names = (/'station_soil', 'ntime_soil '/), &
[4400]	1927	output_type = 'real32' )
	1928	ELSE
	1929
[4504]	1930	return_value = dom_def_var( vmea(l)%nc_filename, variable_name = variable_name, &
	1931	dimension_names = (/'station', 'ntime '/), &
[4400]	1932	output_type = 'real32' )
	1933	ENDIF
	1934	!
[4498]	1935	!-- Set variable attributes. Please note, for some variables not all attributes are defined,
	1936	!-- e.g. standard_name for the horizontal wind components.
[4400]	1937	CALL vm_set_attributes( vmea(l)%var_atts(n) )
	1938
	1939	IF ( vmea(l)%var_atts(n)%long_name /= 'none' ) THEN
[4504]	1940	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1941	attribute_name = char_long, &
[4400]	1942	value = TRIM( vmea(l)%var_atts(n)%long_name ) )
	1943	ENDIF
	1944	IF ( vmea(l)%var_atts(n)%standard_name /= 'none' ) THEN
[4504]	1945	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1946	attribute_name = char_standard, &
[4400]	1947	value = TRIM( vmea(l)%var_atts(n)%standard_name ) )
	1948	ENDIF
	1949	IF ( vmea(l)%var_atts(n)%units /= 'none' ) THEN
[4504]	1950	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1951	attribute_name = char_unit, &
[4400]	1952	value = TRIM( vmea(l)%var_atts(n)%units ) )
	1953	ENDIF
	1954
[4504]	1955	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1956	attribute_name = 'grid_mapping', &
[4400]	1957	value = TRIM( vmea(l)%var_atts(n)%grid_mapping ) )
	1958
[4504]	1959	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1960	attribute_name = 'coordinates', &
[4400]	1961	value = TRIM( vmea(l)%var_atts(n)%coordinates ) )
	1962
[4504]	1963	return_value = dom_def_att( vmea(l)%nc_filename, variable_name = variable_name, &
[4498]	1964	attribute_name = char_fill, &
[4408]	1965	value = REAL( vmea(l)%var_atts(n)%fill_value, KIND=4 ) )
[4400]	1966
	1967	ENDDO ! loop over variables per site
	1968
	1969	ENDDO ! loop over sites
	1970
	1971
	1972	END SUBROUTINE vm_init_output
	1973
[4498]	1974	!--------------------------------------------------------------------------------------------------!
[4400]	1975	! Description:
	1976	! ------------
	1977	!> Parallel NetCDF output via data-output module.
[4498]	1978	!--------------------------------------------------------------------------------------------------!
[4400]	1979	SUBROUTINE vm_data_output
	1980
[4498]	1981	CHARACTER(LEN=100) :: variable_name !< name of output variable
	1982	CHARACTER(LEN=maximum_name_length), DIMENSION(:), ALLOCATABLE :: station_name !< string for station name, consecutively ordered
[4400]	1983
[4421]	1984	CHARACTER(LEN=1), DIMENSION(:,:), ALLOCATABLE, TARGET :: output_values_2d_char_target !< target for output name arrays
	1985	CHARACTER(LEN=1), DIMENSION(:,:), POINTER :: output_values_2d_char_pointer !< pointer for output name arrays
	1986
	1987	INTEGER(iwp) :: l !< loop index for the number of sites
	1988	INTEGER(iwp) :: n !< loop index for observation points
	1989	INTEGER(iwp) :: nn !< loop index for number of characters in a name
[4400]	1990	INTEGER :: return_value !< returned status value of called function
	1991	INTEGER(iwp) :: t_ind !< time index
	1992
[4641]	1993	REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_lat !< transformed geographical coordinate (latitude)
	1994	REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_lon !< transformed geographical coordinate (longitude)
[4498]	1995	REAL(wp), DIMENSION(:), ALLOCATABLE :: oro_rel !< relative altitude of model surface
	1996	REAL(wp), DIMENSION(:), POINTER :: output_values_1d_pointer !< pointer for 1d output array
	1997	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: output_values_1d_target !< target for 1d output array
	1998	REAL(wp), DIMENSION(:,:), POINTER :: output_values_2d_pointer !< pointer for 2d output array
	1999	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: output_values_2d_target !< target for 2d output array
[4400]	2000
[4438]	2001	CALL cpu_log( log_point_s(26), 'VM output', 'start' )
[4400]	2002	!
	2003	!-- At the first call of this routine write the spatial coordinates.
	2004	IF ( .NOT. initial_write_coordinates ) THEN
	2005	!
	2006	!-- Write spatial coordinates.
	2007	DO l = 1, vmea_general%nvm
	2008	!
	2009	!-- Skip if no observations were taken.
	2010	IF ( vmea(l)%ns_tot == 0 .AND. vmea(l)%ns_soil_tot == 0 ) CYCLE
	2011
	2012	ALLOCATE( output_values_1d_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
	2013	!
	2014	!-- Output of Easting coordinate. Before output, recalculate EUTM.
[4498]	2015	output_values_1d_target = init_model%origin_x &
	2016	+ REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
	2017	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
	2018	+ REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
[3913]	2019	* SIN( init_model%rotation_angle * pi / 180.0_wp )
[4400]	2020
	2021	output_values_1d_pointer => output_values_1d_target
	2022
[4504]	2023	return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM', &
	2024	values_realwp_1d = output_values_1d_pointer, &
	2025	bounds_start = (/vmea(l)%start_coord_a/), &
	2026	bounds_end = (/vmea(l)%end_coord_a /) )
[4400]	2027	!
	2028	!-- Output of Northing coordinate. Before output, recalculate NUTM.
[4498]	2029	output_values_1d_target = init_model%origin_y &
	2030	- REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
	2031	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
	2032	+ REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
[3913]	2033	* COS( init_model%rotation_angle * pi / 180.0_wp )
[3704]	2034
[4400]	2035	output_values_1d_pointer => output_values_1d_target
[4504]	2036	return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM', &
	2037	values_realwp_1d = output_values_1d_pointer, &
	2038	bounds_start = (/vmea(l)%start_coord_a/), &
	2039	bounds_end = (/vmea(l)%end_coord_a /) )
[3704]	2040	!
[4641]	2041	!-- Output of longitude and latitude coordinate. Before output, convert it.
	2042	ALLOCATE( dum_lat(1:vmea(l)%ns) )
	2043	ALLOCATE( dum_lon(1:vmea(l)%ns) )
	2044
	2045	DO n = 1, vmea(l)%ns
	2046	CALL convert_utm_to_geographic( crs_list, &
	2047	init_model%origin_x &
	2048	+ REAL( vmea(l)%i(n) + 0.5_wp, KIND = wp ) * dx &
	2049	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
	2050	+ REAL( vmea(l)%j(n) + 0.5_wp, KIND = wp ) * dy &
	2051	* SIN( init_model%rotation_angle * pi / 180.0_wp ), &
	2052	init_model%origin_y &
	2053	- REAL( vmea(l)%i(n) + 0.5_wp, KIND = wp ) * dx &
	2054	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
	2055	+ REAL( vmea(l)%j(n) + 0.5_wp, KIND = wp ) * dy &
	2056	* COS( init_model%rotation_angle * pi / 180.0_wp ), &
	2057	dum_lon(n), dum_lat(n) )
	2058	ENDDO
	2059
	2060	output_values_1d_target = dum_lat
	2061	output_values_1d_pointer => output_values_1d_target
	2062	return_value = dom_write_var( vmea(l)%nc_filename, 'lat', &
	2063	values_realwp_1d = output_values_1d_pointer, &
	2064	bounds_start = (/vmea(l)%start_coord_a/), &
	2065	bounds_end = (/vmea(l)%end_coord_a /) )
	2066
	2067	output_values_1d_target = dum_lon
	2068	output_values_1d_pointer => output_values_1d_target
	2069	return_value = dom_write_var( vmea(l)%nc_filename, 'lon', &
	2070	values_realwp_1d = output_values_1d_pointer, &
	2071	bounds_start = (/vmea(l)%start_coord_a/), &
	2072	bounds_end = (/vmea(l)%end_coord_a /) )
	2073	DEALLOCATE( dum_lat )
	2074	DEALLOCATE( dum_lon )
	2075	!
[4400]	2076	!-- Output of relative height coordinate.
[4498]	2077	!-- Before this is output, first define the relative orographie height and add this to z.
[4400]	2078	ALLOCATE( oro_rel(1:vmea(l)%ns) )
	2079	DO n = 1, vmea(l)%ns
	2080	oro_rel(n) = zw(topo_top_ind(vmea(l)%j(n),vmea(l)%i(n),3))
	2081	ENDDO
	2082
	2083	output_values_1d_target = vmea(l)%zar(1:vmea(l)%ns) + oro_rel(:)
	2084	output_values_1d_pointer => output_values_1d_target
[4504]	2085	return_value = dom_write_var( vmea(l)%nc_filename, 'z', &
	2086	values_realwp_1d = output_values_1d_pointer, &
	2087	bounds_start = (/vmea(l)%start_coord_a/), &
	2088	bounds_end = (/vmea(l)%end_coord_a /) )
[3704]	2089	!
[4498]	2090	!-- Write surface altitude for the station. Note, since z is already a relative observation
	2091	!-- height, station_h must be zero, in order to obtain the observation level.
[4400]	2092	output_values_1d_target = oro_rel(:)
	2093	output_values_1d_pointer => output_values_1d_target
[4504]	2094	return_value = dom_write_var( vmea(l)%nc_filename, 'station_h', &
	2095	values_realwp_1d = output_values_1d_pointer, &
	2096	bounds_start = (/vmea(l)%start_coord_a/), &
	2097	bounds_end = (/vmea(l)%end_coord_a /) )
[3704]	2098
[4400]	2099	DEALLOCATE( oro_rel )
	2100	DEALLOCATE( output_values_1d_target )
[3704]	2101	!
[4421]	2102	!-- Write station name
	2103	ALLOCATE ( station_name(vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
	2104	ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_a:vmea(l)%end_coord_a, &
	2105	1:maximum_name_length) )
	2106
	2107	DO n = vmea(l)%start_coord_a, vmea(l)%end_coord_a
	2108	station_name(n) = REPEAT( ' ', maximum_name_length )
	2109	WRITE( station_name(n), '(A,I10.10)') "station", n
	2110	DO nn = 1, maximum_name_length
	2111	output_values_2d_char_target(n,nn) = station_name(n)(nn:nn)
	2112	ENDDO
	2113	ENDDO
	2114
	2115	output_values_2d_char_pointer => output_values_2d_char_target
	2116
[4504]	2117	return_value = dom_write_var( vmea(l)%nc_filename, 'station_name', &
	2118	values_char_2d = output_values_2d_char_pointer, &
	2119	bounds_start = (/ 1, vmea(l)%start_coord_a /), &
	2120	bounds_end = (/ maximum_name_length, &
	2121	vmea(l)%end_coord_a /) )
[4421]	2122
	2123	DEALLOCATE( station_name )
	2124	DEALLOCATE( output_values_2d_char_target )
	2125	!
[4498]	2126	!-- In case of sampled soil quantities, output also the respective coordinate arrays.
[4400]	2127	IF ( vmea(l)%soil_sampling ) THEN
	2128	ALLOCATE( output_values_1d_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
	2129	!
	2130	!-- Output of Easting coordinate. Before output, recalculate EUTM.
[4498]	2131	output_values_1d_target = init_model%origin_x &
[4645]	2132	+ REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
[4498]	2133	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
[4645]	2134	+ REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
[4400]	2135	* SIN( init_model%rotation_angle * pi / 180.0_wp )
	2136	output_values_1d_pointer => output_values_1d_target
[4504]	2137	return_value = dom_write_var( vmea(l)%nc_filename, 'E_UTM_soil', &
	2138	values_realwp_1d = output_values_1d_pointer, &
	2139	bounds_start = (/vmea(l)%start_coord_s/), &
	2140	bounds_end = (/vmea(l)%end_coord_s /) )
[4400]	2141	!
	2142	!-- Output of Northing coordinate. Before output, recalculate NUTM.
[4498]	2143	output_values_1d_target = init_model%origin_y &
[4641]	2144	- REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dx &
[4498]	2145	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
[4641]	2146	+ REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, KIND = wp ) * dy &
[4400]	2147	* COS( init_model%rotation_angle * pi / 180.0_wp )
	2148
	2149	output_values_1d_pointer => output_values_1d_target
[4504]	2150	return_value = dom_write_var( vmea(l)%nc_filename, 'N_UTM_soil', &
	2151	values_realwp_1d = output_values_1d_pointer, &
	2152	bounds_start = (/vmea(l)%start_coord_s/), &
	2153	bounds_end = (/vmea(l)%end_coord_s /) )
[4400]	2154	!
[4641]	2155	!-- Output of longitude and latitude coordinate. Before output, convert it.
	2156	ALLOCATE( dum_lat(1:vmea(l)%ns_soil) )
	2157	ALLOCATE( dum_lon(1:vmea(l)%ns_soil) )
	2158
	2159	DO n = 1, vmea(l)%ns_soil
	2160	CALL convert_utm_to_geographic( crs_list, &
	2161	init_model%origin_x &
	2162	+ REAL( vmea(l)%i_soil(n) + 0.5_wp, KIND = wp ) * dx &
	2163	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
	2164	+ REAL( vmea(l)%j_soil(n) + 0.5_wp, KIND = wp ) * dy &
	2165	* SIN( init_model%rotation_angle * pi / 180.0_wp ), &
	2166	init_model%origin_y &
	2167	- REAL( vmea(l)%i_soil(n) + 0.5_wp, KIND = wp ) * dx &
	2168	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
	2169	+ REAL( vmea(l)%j_soil(n) + 0.5_wp, KIND = wp ) * dy &
	2170	* COS( init_model%rotation_angle * pi / 180.0_wp ), &
	2171	dum_lon(n), dum_lat(n) )
	2172	ENDDO
	2173
	2174	output_values_1d_target = dum_lat
	2175	output_values_1d_pointer => output_values_1d_target
	2176	return_value = dom_write_var( vmea(l)%nc_filename, 'lat_soil', &
	2177	values_realwp_1d = output_values_1d_pointer, &
	2178	bounds_start = (/vmea(l)%start_coord_s/), &
	2179	bounds_end = (/vmea(l)%end_coord_s /) )
	2180
	2181	output_values_1d_target = dum_lon
	2182	output_values_1d_pointer => output_values_1d_target
	2183	return_value = dom_write_var( vmea(l)%nc_filename, 'lon_soil', &
	2184	values_realwp_1d = output_values_1d_pointer, &
	2185	bounds_start = (/vmea(l)%start_coord_s/), &
	2186	bounds_end = (/vmea(l)%end_coord_s /) )
	2187	DEALLOCATE( dum_lat )
	2188	DEALLOCATE( dum_lon )
	2189	!
[4400]	2190	!-- Output of relative height coordinate.
[4498]	2191	!-- Before this is output, first define the relative orographie height and add this to z.
[4400]	2192	ALLOCATE( oro_rel(1:vmea(l)%ns_soil) )
[4406]	2193	DO n = 1, vmea(l)%ns_soil
[4400]	2194	oro_rel(n) = zw(topo_top_ind(vmea(l)%j_soil(n),vmea(l)%i_soil(n),3))
	2195	ENDDO
	2196
	2197	output_values_1d_target = vmea(l)%depth(1:vmea(l)%ns_soil) + oro_rel(:)
	2198	output_values_1d_pointer => output_values_1d_target
[4504]	2199	return_value = dom_write_var( vmea(l)%nc_filename, 'z_soil', &
	2200	values_realwp_1d = output_values_1d_pointer, &
	2201	bounds_start = (/vmea(l)%start_coord_s/), &
	2202	bounds_end = (/vmea(l)%end_coord_s /) )
[4400]	2203	!
[4498]	2204	!-- Write surface altitude for the station. Note, since z is already a relative observation
	2205	!-- height, station_h must be zero, in order to obtain the observation level.
[4400]	2206	output_values_1d_target = oro_rel(:)
	2207	output_values_1d_pointer => output_values_1d_target
[4504]	2208	return_value = dom_write_var( vmea(l)%nc_filename, 'station_h_soil', &
	2209	values_realwp_1d = output_values_1d_pointer, &
	2210	bounds_start = (/vmea(l)%start_coord_s/), &
	2211	bounds_end = (/vmea(l)%end_coord_s /) )
[4400]	2212
	2213	DEALLOCATE( oro_rel )
	2214	DEALLOCATE( output_values_1d_target )
	2215	!
[4421]	2216	!-- Write station name
	2217	ALLOCATE ( station_name(vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
[4498]	2218	ALLOCATE ( output_values_2d_char_target(vmea(l)%start_coord_s:vmea(l)%end_coord_s, &
[4421]	2219	1:maximum_name_length) )
[4408]	2220
[4421]	2221	DO n = vmea(l)%start_coord_s, vmea(l)%end_coord_s
	2222	station_name(n) = REPEAT( ' ', maximum_name_length )
	2223	WRITE( station_name(n), '(A,I10.10)') "station", n
	2224	DO nn = 1, maximum_name_length
	2225	output_values_2d_char_target(n,nn) = station_name(n)(nn:nn)
	2226	ENDDO
	2227	ENDDO
	2228	output_values_2d_char_pointer => output_values_2d_char_target
	2229
[4504]	2230	return_value = dom_write_var( vmea(l)%nc_filename, 'station_name_soil', &
	2231	values_char_2d = output_values_2d_char_pointer, &
	2232	bounds_start = (/ 1, vmea(l)%start_coord_s /), &
	2233	bounds_end = (/ maximum_name_length, &
	2234	vmea(l)%end_coord_s /) )
[4421]	2235
	2236	DEALLOCATE( station_name )
	2237	DEALLOCATE( output_values_2d_char_target )
	2238
[4400]	2239	ENDIF
	2240
	2241	ENDDO ! loop over sites
	2242
	2243	initial_write_coordinates = .TRUE.
	2244	ENDIF
	2245	!
	2246	!-- Loop over all sites.
	2247	DO l = 1, vmea_general%nvm
	2248	!
	2249	!-- Skip if no observations were taken.
	2250	IF ( vmea(l)%ns_tot == 0 .AND. vmea(l)%ns_soil_tot == 0 ) CYCLE
	2251	!
	2252	!-- Determine time index in file.
	2253	t_ind = vmea(l)%file_time_index + 1
	2254	!
	2255	!-- Write output variables. Distinguish between atmosphere and soil variables.
	2256	DO n = 1, vmea(l)%nmeas
[4498]	2257	IF ( vmea(l)%soil_sampling .AND. &
[4400]	2258	ANY( TRIM( vmea(l)%var_atts(n)%name) == soil_vars ) ) THEN
	2259	!
	2260	!-- Write time coordinate to file
	2261	variable_name = 'time_soil'
	2262	ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_s:vmea(l)%end_coord_s) )
	2263	output_values_2d_target(t_ind,:) = time_since_reference_point
	2264	output_values_2d_pointer => output_values_2d_target
	2265
[4504]	2266	return_value = dom_write_var( vmea(l)%nc_filename, variable_name, &
	2267	values_realwp_2d = output_values_2d_pointer, &
[4498]	2268	bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
[4400]	2269	bounds_end = (/vmea(l)%end_coord_s, t_ind /) )
	2270
	2271	variable_name = TRIM( vmea(l)%var_atts(n)%name )
	2272	output_values_2d_target(t_ind,:) = vmea(l)%measured_vars_soil(:,n)
	2273	output_values_2d_pointer => output_values_2d_target
[4504]	2274	return_value = dom_write_var( vmea(l)%nc_filename, variable_name, &
	2275	values_realwp_2d = output_values_2d_pointer, &
	2276	bounds_start = (/vmea(l)%start_coord_s, t_ind/), &
	2277	bounds_end = (/vmea(l)%end_coord_s, t_ind /) )
[4400]	2278	DEALLOCATE( output_values_2d_target )
	2279	ELSE
	2280	!
	2281	!-- Write time coordinate to file
	2282	variable_name = 'time'
	2283	ALLOCATE( output_values_2d_target(t_ind:t_ind,vmea(l)%start_coord_a:vmea(l)%end_coord_a) )
	2284	output_values_2d_target(t_ind,:) = time_since_reference_point
	2285	output_values_2d_pointer => output_values_2d_target
	2286
[4504]	2287	return_value = dom_write_var( vmea(l)%nc_filename, variable_name, &
[4498]	2288	values_realwp_2d = output_values_2d_pointer, &
	2289	bounds_start = (/vmea(l)%start_coord_a, t_ind/), &
[4400]	2290	bounds_end = (/vmea(l)%end_coord_a, t_ind/) )
	2291
	2292	variable_name = TRIM( vmea(l)%var_atts(n)%name )
	2293
	2294	output_values_2d_target(t_ind,:) = vmea(l)%measured_vars(:,n)
	2295	output_values_2d_pointer => output_values_2d_target
[4504]	2296	return_value = dom_write_var( vmea(l)%nc_filename, variable_name, &
	2297	values_realwp_2d = output_values_2d_pointer, &
	2298	bounds_start = (/ vmea(l)%start_coord_a, t_ind /), &
	2299	bounds_end = (/ vmea(l)%end_coord_a, t_ind /) )
[4400]	2300
	2301	DEALLOCATE( output_values_2d_target )
	2302	ENDIF
	2303	ENDDO
	2304	!
	2305	!-- Update number of written time indices
	2306	vmea(l)%file_time_index = t_ind
	2307
	2308	ENDDO ! loop over sites
	2309
[4438]	2310	CALL cpu_log( log_point_s(26), 'VM output', 'stop' )
[4400]	2311
[4438]	2312
[4498]	2313	END SUBROUTINE vm_data_output
[4400]	2314
[4498]	2315	!--------------------------------------------------------------------------------------------------!
[3704]	2316	! Description:
	2317	! ------------
[3434]	2318	!> Sampling of the actual quantities along the observation coordinates
[4498]	2319	!--------------------------------------------------------------------------------------------------!
	2320	SUBROUTINE vm_sampling
[3434]	2321
[4498]	2322	USE radiation_model_mod, &
[4400]	2323	ONLY: radiation
[3522]	2324
[4498]	2325	USE surface_mod, &
	2326	ONLY: surf_def_h, &
	2327	surf_lsm_h, &
[4400]	2328	surf_usm_h
	2329
[3704]	2330	INTEGER(iwp) :: i !< grid index in x-direction
	2331	INTEGER(iwp) :: j !< grid index in y-direction
	2332	INTEGER(iwp) :: k !< grid index in z-direction
	2333	INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
	2334	INTEGER(iwp) :: l !< running index over the number of stations
	2335	INTEGER(iwp) :: m !< running index over all virtual observation coordinates
	2336	INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate
	2337	INTEGER(iwp) :: n !< running index over all measured variables at a station
	2338	INTEGER(iwp) :: nn !< running index over the number of chemcal species
[4400]	2339
[4498]	2340	LOGICAL :: match_lsm !< flag indicating natural-type surface
	2341	LOGICAL :: match_usm !< flag indicating urban-type surface
[4400]	2342
[4498]	2343	REAL(wp) :: e_s !< saturation water vapor pressure
	2344	REAL(wp) :: q_s !< saturation mixing ratio
	2345	REAL(wp) :: q_wv !< mixing ratio
[4438]	2346
	2347	CALL cpu_log( log_point_s(27), 'VM sampling', 'start' )
[3434]	2348	!
[4400]	2349	!-- Loop over all sites.
[3704]	2350	DO l = 1, vmea_general%nvm
[3434]	2351	!
[3704]	2352	!-- At the beginning, set _FillValues
[4498]	2353	IF ( ALLOCATED( vmea(l)%measured_vars ) ) vmea(l)%measured_vars = vmea(l)%fillout
	2354	IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) vmea(l)%measured_vars_soil = vmea(l)%fillout
[3704]	2355	!
[4400]	2356	!-- Loop over all variables measured at this site.
[3833]	2357	DO n = 1, vmea(l)%nmeas
[4400]	2358
	2359	SELECT CASE ( TRIM( vmea(l)%var_atts(n)%name ) )
	2360
	2361	CASE ( 'theta' ) ! potential temperature
[3522]	2362	IF ( .NOT. neutral ) THEN
	2363	DO m = 1, vmea(l)%ns
	2364	k = vmea(l)%k(m)
	2365	j = vmea(l)%j(m)
	2366	i = vmea(l)%i(m)
[3704]	2367	vmea(l)%measured_vars(m,n) = pt(k,j,i)
[3522]	2368	ENDDO
	2369	ENDIF
[4400]	2370
	2371	CASE ( 'ta' ) ! absolute temperature
[3522]	2372	IF ( .NOT. neutral ) THEN
	2373	DO m = 1, vmea(l)%ns
	2374	k = vmea(l)%k(m)
	2375	j = vmea(l)%j(m)
	2376	i = vmea(l)%i(m)
[4498]	2377	vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k ) - degc_to_k
[3522]	2378	ENDDO
	2379	ENDIF
[4400]	2380
[3704]	2381	CASE ( 't_va' )
[4400]	2382
	2383	CASE ( 'hus' ) ! mixing ratio
[3522]	2384	IF ( humidity ) THEN
	2385	DO m = 1, vmea(l)%ns
	2386	k = vmea(l)%k(m)
	2387	j = vmea(l)%j(m)
	2388	i = vmea(l)%i(m)
[3704]	2389	vmea(l)%measured_vars(m,n) = q(k,j,i)
[3522]	2390	ENDDO
	2391	ENDIF
[4400]	2392
	2393	CASE ( 'haa' ) ! absolute humidity
	2394	IF ( humidity ) THEN
	2395	DO m = 1, vmea(l)%ns
	2396	k = vmea(l)%k(m)
	2397	j = vmea(l)%j(m)
	2398	i = vmea(l)%i(m)
[4498]	2399	vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
[4400]	2400	ENDDO
	2401	ENDIF
	2402
	2403	CASE ( 'pwv' ) ! water vapor partial pressure
	2404	IF ( humidity ) THEN
	2405	! DO m = 1, vmea(l)%ns
	2406	! k = vmea(l)%k(m)
	2407	! j = vmea(l)%j(m)
	2408	! i = vmea(l)%i(m)
[4498]	2409	! vmea(l)%measured_vars(m,n) = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) &
	2410	! * rho_air(k)
[4400]	2411	! ENDDO
	2412	ENDIF
	2413
	2414	CASE ( 'hur' ) ! relative humidity
	2415	IF ( humidity ) THEN
	2416	DO m = 1, vmea(l)%ns
	2417	k = vmea(l)%k(m)
	2418	j = vmea(l)%j(m)
	2419	i = vmea(l)%i(m)
	2420	!
[4498]	2421	!-- Calculate actual temperature, water vapor saturation pressure and, based on
	2422	!-- this, the saturation mixing ratio.
[4400]	2423	e_s = magnus( exner(k) * pt(k,j,i) )
	2424	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
	2425	q_wv = ( q(k,j,i) / ( 1.0_wp - q(k,j,i) ) ) * rho_air(k)
	2426
	2427	vmea(l)%measured_vars(m,n) = q_wv / ( q_s + 1E-10_wp )
	2428	ENDDO
	2429	ENDIF
	2430
	2431	CASE ( 'u', 'ua' ) ! u-component
[3522]	2432	DO m = 1, vmea(l)%ns
	2433	k = vmea(l)%k(m)
	2434	j = vmea(l)%j(m)
	2435	i = vmea(l)%i(m)
[3704]	2436	vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
[3522]	2437	ENDDO
[4400]	2438
	2439	CASE ( 'v', 'va' ) ! v-component
[3522]	2440	DO m = 1, vmea(l)%ns
	2441	k = vmea(l)%k(m)
	2442	j = vmea(l)%j(m)
	2443	i = vmea(l)%i(m)
[3704]	2444	vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
[3522]	2445	ENDDO
[4400]	2446
	2447	CASE ( 'w' ) ! w-component
[3522]	2448	DO m = 1, vmea(l)%ns
[4400]	2449	k = MAX ( 1, vmea(l)%k(m) )
[3522]	2450	j = vmea(l)%j(m)
	2451	i = vmea(l)%i(m)
[3704]	2452	vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[3522]	2453	ENDDO
[4400]	2454
	2455	CASE ( 'wspeed' ) ! horizontal wind speed
[3522]	2456	DO m = 1, vmea(l)%ns
	2457	k = vmea(l)%k(m)
	2458	j = vmea(l)%j(m)
	2459	i = vmea(l)%i(m)
[4504]	2460	vmea(l)%measured_vars(m,n) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 &
	2461	+ ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
[3522]	2462	)
	2463	ENDDO
[4400]	2464
	2465	CASE ( 'wdir' ) ! wind direction
[3522]	2466	DO m = 1, vmea(l)%ns
	2467	k = vmea(l)%k(m)
	2468	j = vmea(l)%j(m)
	2469	i = vmea(l)%i(m)
[4400]	2470
[4498]	2471	vmea(l)%measured_vars(m,n) = 180.0_wp + 180.0_wp / pi * ATAN2( &
[4504]	2472	0.5_wp * ( v(k,j,i) + v(k,j+1,i) ), &
	2473	0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) &
	2474	)
[3522]	2475	ENDDO
[4400]	2476
[3704]	2477	CASE ( 'utheta' )
	2478	DO m = 1, vmea(l)%ns
	2479	k = vmea(l)%k(m)
	2480	j = vmea(l)%j(m)
	2481	i = vmea(l)%i(m)
[4498]	2482	vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * pt(k,j,i)
[3704]	2483	ENDDO
[4400]	2484
[3704]	2485	CASE ( 'vtheta' )
	2486	DO m = 1, vmea(l)%ns
	2487	k = vmea(l)%k(m)
	2488	j = vmea(l)%j(m)
	2489	i = vmea(l)%i(m)
[4498]	2490	vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * pt(k,j,i)
[3704]	2491	ENDDO
[4400]	2492
[3704]	2493	CASE ( 'wtheta' )
	2494	DO m = 1, vmea(l)%ns
	2495	k = MAX ( 1, vmea(l)%k(m) )
	2496	j = vmea(l)%j(m)
	2497	i = vmea(l)%i(m)
[4498]	2498	vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * pt(k,j,i)
[3704]	2499	ENDDO
[4400]	2500
	2501	CASE ( 'uqv' )
	2502	IF ( humidity ) THEN
	2503	DO m = 1, vmea(l)%ns
	2504	k = vmea(l)%k(m)
	2505	j = vmea(l)%j(m)
	2506	i = vmea(l)%i(m)
[4498]	2507	vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) * q(k,j,i)
[4400]	2508	ENDDO
	2509	ENDIF
	2510
	2511	CASE ( 'vqv' )
	2512	IF ( humidity ) THEN
	2513	DO m = 1, vmea(l)%ns
	2514	k = vmea(l)%k(m)
	2515	j = vmea(l)%j(m)
	2516	i = vmea(l)%i(m)
[4498]	2517	vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) * q(k,j,i)
[4400]	2518	ENDDO
	2519	ENDIF
	2520
	2521	CASE ( 'wqv' )
	2522	IF ( humidity ) THEN
	2523	DO m = 1, vmea(l)%ns
	2524	k = MAX ( 1, vmea(l)%k(m) )
	2525	j = vmea(l)%j(m)
	2526	i = vmea(l)%i(m)
[4498]	2527	vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k-1,j,i) + w(k,j,i) ) * q(k,j,i)
[4400]	2528	ENDDO
	2529	ENDIF
	2530
[3704]	2531	CASE ( 'uw' )
	2532	DO m = 1, vmea(l)%ns
	2533	k = MAX ( 1, vmea(l)%k(m) )
	2534	j = vmea(l)%j(m)
	2535	i = vmea(l)%i(m)
[4498]	2536	vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) * &
	2537	( u(k,j,i) + u(k,j,i+1) )
[3704]	2538	ENDDO
[4400]	2539
[3704]	2540	CASE ( 'vw' )
	2541	DO m = 1, vmea(l)%ns
	2542	k = MAX ( 1, vmea(l)%k(m) )
	2543	j = vmea(l)%j(m)
	2544	i = vmea(l)%i(m)
[4498]	2545	vmea(l)%measured_vars(m,n) = 0.25_wp * ( w(k-1,j,i) + w(k,j,i) ) * &
	2546	( v(k,j,i) + v(k,j+1,i) )
[3704]	2547	ENDDO
[4400]	2548
[3704]	2549	CASE ( 'uv' )
	2550	DO m = 1, vmea(l)%ns
[4400]	2551	k = vmea(l)%k(m)
[3704]	2552	j = vmea(l)%j(m)
	2553	i = vmea(l)%i(m)
[4498]	2554	vmea(l)%measured_vars(m,n) = 0.25_wp * ( u(k,j,i) + u(k,j,i+1) ) * &
	2555	( v(k,j,i) + v(k,j+1,i) )
[3704]	2556	ENDDO
[3522]	2557	!
[4498]	2558	!-- Chemistry variables. List of variables that may need extension. Note, gas species in
	2559	!-- PALM are in ppm and no distinction is made between mole-fraction and concentration
	2560	!-- quantities (all are output in ppm so far).
	2561	CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfno', 'mfno2', 'mcno', 'mcno2', 'tro3' )
[3704]	2562	IF ( air_chemistry ) THEN
	2563	!
[4400]	2564	!-- First, search for the measured variable in the chem_vars
	2565	!-- list, in order to get the internal name of the variable.
[3704]	2566	DO nn = 1, UBOUND( chem_vars, 2 )
[4498]	2567	IF ( TRIM( vmea(l)%var_atts(n)%name ) == &
[3704]	2568	TRIM( chem_vars(0,nn) ) ) ind_chem = nn
	2569	ENDDO
	2570	!
[4498]	2571	!-- Run loop over all chemical species, if the measured variable matches the interal
	2572	!-- name, sample the variable. Note, nvar as a chemistry-module variable.
[4400]	2573	DO nn = 1, nvar
[4498]	2574	IF ( TRIM( chem_vars(1,ind_chem) ) == TRIM( chem_species(nn)%name ) ) THEN
[4400]	2575	DO m = 1, vmea(l)%ns
[3522]	2576	k = vmea(l)%k(m)
	2577	j = vmea(l)%j(m)
[4400]	2578	i = vmea(l)%i(m)
[4498]	2579	vmea(l)%measured_vars(m,n) = chem_species(nn)%conc(k,j,i)
[3522]	2580	ENDDO
	2581	ENDIF
	2582	ENDDO
	2583	ENDIF
[4400]	2584
	2585	CASE ( 'us' ) ! friction velocity
[3522]	2586	DO m = 1, vmea(l)%ns
	2587	!
[4498]	2588	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2589	!-- limit the indices.
[4400]	2590	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2591	j = MERGE( j , nyn, j < nyn )
	2592	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2593	i = MERGE( i , nxr, i < nxr )
	2594
[4498]	2595	DO mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
[3704]	2596	vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
[3522]	2597	ENDDO
[4671]	2598	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2599	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%us(mm)
[3522]	2600	ENDDO
[4671]	2601	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2602	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%us(mm)
[3522]	2603	ENDDO
	2604	ENDDO
[4400]	2605
	2606	CASE ( 'thetas' ) ! scaling parameter temperature
[3522]	2607	DO m = 1, vmea(l)%ns
	2608	!
[4498]	2609	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2610	!- limit the indices.
[4400]	2611	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2612	j = MERGE( j , nyn, j < nyn )
	2613	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2614	i = MERGE( i , nxr, i < nxr )
	2615
[4498]	2616	DO mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
[3704]	2617	vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
[3522]	2618	ENDDO
[4671]	2619	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2620	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%ts(mm)
[3522]	2621	ENDDO
[4671]	2622	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2623	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%ts(mm)
[3522]	2624	ENDDO
	2625	ENDDO
[4400]	2626
	2627	CASE ( 'hfls' ) ! surface latent heat flux
[3522]	2628	DO m = 1, vmea(l)%ns
	2629	!
[4498]	2630	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2631	!-- limit the indices.
[4400]	2632	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2633	j = MERGE( j , nyn, j < nyn )
	2634	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2635	i = MERGE( i , nxr, i < nxr )
	2636
[4498]	2637	DO mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
[3704]	2638	vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
[3522]	2639	ENDDO
[4671]	2640	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2641	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%qsws(mm)
[3522]	2642	ENDDO
[4671]	2643	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2644	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%qsws(mm)
[3522]	2645	ENDDO
	2646	ENDDO
[4400]	2647
	2648	CASE ( 'hfss' ) ! surface sensible heat flux
[3522]	2649	DO m = 1, vmea(l)%ns
	2650	!
[4498]	2651	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2652	!-- limit the indices.
[4400]	2653	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2654	j = MERGE( j , nyn, j < nyn )
	2655	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2656	i = MERGE( i , nxr, i < nxr )
	2657
[4498]	2658	DO mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
[3704]	2659	vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
[3522]	2660	ENDDO
[4671]	2661	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2662	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%shf(mm)
[3522]	2663	ENDDO
[4671]	2664	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2665	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%shf(mm)
[3522]	2666	ENDDO
	2667	ENDDO
[4400]	2668
	2669	CASE ( 'hfdg' ) ! ground heat flux
	2670	DO m = 1, vmea(l)%ns
	2671	!
[4498]	2672	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2673	!-- limit the indices.
[4400]	2674	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2675	j = MERGE( j , nyn, j < nyn )
	2676	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2677	i = MERGE( i , nxr, i < nxr )
	2678
[4671]	2679	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2680	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%ghf(mm)
[4400]	2681	ENDDO
	2682	ENDDO
	2683
	2684	CASE ( 'lwcs' ) ! liquid water of soil layer
	2685	! DO m = 1, vmea(l)%ns
	2686	! !
[4498]	2687	! !-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2688	! !-- limit the indices.
[4400]	2689	! j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2690	! j = MERGE( j , nyn, j < nyn )
	2691	! i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2692	! i = MERGE( i , nxr, i < nxr )
[4408]	2693	!
[4671]	2694	! DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
[4400]	2695	! vmea(l)%measured_vars(m,n) = ?
	2696	! ENDDO
	2697	! ENDDO
	2698
	2699	CASE ( 'rnds' ) ! surface net radiation
[3522]	2700	IF ( radiation ) THEN
	2701	DO m = 1, vmea(l)%ns
	2702	!
[4498]	2703	!-- Surface data is only available on inner subdomains, not on ghost points.
	2704	!-- Hence, limit the indices.
[4400]	2705	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2706	j = MERGE( j , nyn, j < nyn )
	2707	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2708	i = MERGE( i , nxr, i < nxr )
	2709
[4671]	2710	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2711	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%rad_net(mm)
[3522]	2712	ENDDO
[4671]	2713	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2714	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%rad_net(mm)
[3522]	2715	ENDDO
	2716	ENDDO
	2717	ENDIF
[4400]	2718
	2719	CASE ( 'rsus' ) ! surface shortwave out
[3522]	2720	IF ( radiation ) THEN
	2721	DO m = 1, vmea(l)%ns
	2722	!
[4498]	2723	!-- Surface data is only available on inner subdomains, not on ghost points.
	2724	!-- Hence, limit the indices.
[4400]	2725	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2726	j = MERGE( j , nyn, j < nyn )
	2727	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2728	i = MERGE( i , nxr, i < nxr )
	2729
[4671]	2730	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2731	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%rad_sw_out(mm)
[3522]	2732	ENDDO
[4671]	2733	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2734	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%rad_sw_out(mm)
[3522]	2735	ENDDO
	2736	ENDDO
	2737	ENDIF
[4400]	2738
	2739	CASE ( 'rsds' ) ! surface shortwave in
[3522]	2740	IF ( radiation ) THEN
	2741	DO m = 1, vmea(l)%ns
	2742	!
[4498]	2743	!-- Surface data is only available on inner subdomains, not on ghost points.
	2744	!-- Hence, limit the indices.
[4400]	2745	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2746	j = MERGE( j , nyn, j < nyn )
	2747	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2748	i = MERGE( i , nxr, i < nxr )
	2749
[4671]	2750	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2751	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%rad_sw_in(mm)
[3522]	2752	ENDDO
[4671]	2753	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2754	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%rad_sw_in(mm)
[3522]	2755	ENDDO
	2756	ENDDO
	2757	ENDIF
[4400]	2758
	2759	CASE ( 'rlus' ) ! surface longwave out
[3522]	2760	IF ( radiation ) THEN
	2761	DO m = 1, vmea(l)%ns
	2762	!
[4498]	2763	!-- Surface data is only available on inner subdomains, not on ghost points.
	2764	!-- Hence, limit the indices.
[4400]	2765	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2766	j = MERGE( j , nyn, j < nyn )
	2767	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2768	i = MERGE( i , nxr, i < nxr )
	2769
[4671]	2770	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2771	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%rad_lw_out(mm)
[3522]	2772	ENDDO
[4671]	2773	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2774	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%rad_lw_out(mm)
[3522]	2775	ENDDO
	2776	ENDDO
	2777	ENDIF
[4400]	2778
	2779	CASE ( 'rlds' ) ! surface longwave in
[3522]	2780	IF ( radiation ) THEN
	2781	DO m = 1, vmea(l)%ns
	2782	!
[4498]	2783	!-- Surface data is only available on inner subdomains, not on ghost points.
	2784	!-- Hence, limit the indices.
[4400]	2785	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2786	j = MERGE( j , nyn, j < nyn )
	2787	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2788	i = MERGE( i , nxr, i < nxr )
	2789
[4671]	2790	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2791	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%rad_lw_in(mm)
[3522]	2792	ENDDO
[4671]	2793	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2794	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%rad_lw_in(mm)
[3522]	2795	ENDDO
	2796	ENDDO
	2797	ENDIF
[4400]	2798
	2799	CASE ( 'rsd' ) ! shortwave in
[3704]	2800	IF ( radiation ) THEN
[4400]	2801	IF ( radiation_scheme /= 'rrtmg' ) THEN
	2802	DO m = 1, vmea(l)%ns
	2803	k = 0
	2804	j = vmea(l)%j(m)
	2805	i = vmea(l)%i(m)
	2806	vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
	2807	ENDDO
	2808	ELSE
	2809	DO m = 1, vmea(l)%ns
	2810	k = vmea(l)%k(m)
	2811	j = vmea(l)%j(m)
	2812	i = vmea(l)%i(m)
	2813	vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
	2814	ENDDO
	2815	ENDIF
[3704]	2816	ENDIF
[4400]	2817
	2818	CASE ( 'rsu' ) ! shortwave out
[3704]	2819	IF ( radiation ) THEN
[4400]	2820	IF ( radiation_scheme /= 'rrtmg' ) THEN
	2821	DO m = 1, vmea(l)%ns
	2822	k = 0
	2823	j = vmea(l)%j(m)
	2824	i = vmea(l)%i(m)
	2825	vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
	2826	ENDDO
	2827	ELSE
	2828	DO m = 1, vmea(l)%ns
	2829	k = vmea(l)%k(m)
	2830	j = vmea(l)%j(m)
	2831	i = vmea(l)%i(m)
	2832	vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
	2833	ENDDO
	2834	ENDIF
[3704]	2835	ENDIF
[4400]	2836
	2837	CASE ( 'rlu' ) ! longwave out
[3704]	2838	IF ( radiation ) THEN
[4400]	2839	IF ( radiation_scheme /= 'rrtmg' ) THEN
	2840	DO m = 1, vmea(l)%ns
	2841	k = 0
	2842	j = vmea(l)%j(m)
	2843	i = vmea(l)%i(m)
	2844	vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
	2845	ENDDO
	2846	ELSE
	2847	DO m = 1, vmea(l)%ns
	2848	k = vmea(l)%k(m)
	2849	j = vmea(l)%j(m)
	2850	i = vmea(l)%i(m)
	2851	vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
	2852	ENDDO
	2853	ENDIF
[3704]	2854	ENDIF
[4400]	2855
	2856	CASE ( 'rld' ) ! longwave in
[3704]	2857	IF ( radiation ) THEN
[4400]	2858	IF ( radiation_scheme /= 'rrtmg' ) THEN
	2859	DO m = 1, vmea(l)%ns
	2860	k = 0
	2861	!
[4498]	2862	!-- Surface data is only available on inner subdomains, not on ghost points.
	2863	!-- Hence, limit the indices.
[4400]	2864	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2865	j = MERGE( j , nyn, j < nyn )
	2866	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2867	i = MERGE( i , nxr, i < nxr )
	2868
	2869	vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
	2870	ENDDO
	2871	ELSE
	2872	DO m = 1, vmea(l)%ns
	2873	k = vmea(l)%k(m)
	2874	j = vmea(l)%j(m)
	2875	i = vmea(l)%i(m)
	2876	vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
	2877	ENDDO
	2878	ENDIF
[3704]	2879	ENDIF
[4400]	2880
	2881	CASE ( 'rsddif' ) ! shortwave in, diffuse part
[3704]	2882	IF ( radiation ) THEN
	2883	DO m = 1, vmea(l)%ns
	2884	j = vmea(l)%j(m)
	2885	i = vmea(l)%i(m)
[4400]	2886
[3704]	2887	vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
	2888	ENDDO
	2889	ENDIF
[4400]	2890
	2891	CASE ( 't_soil' ) ! soil and wall temperature
[3704]	2892	DO m = 1, vmea(l)%ns_soil
[4400]	2893	j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
	2894	j = MERGE( j , nyn, j < nyn )
	2895	i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
	2896	i = MERGE( i , nxr, i < nxr )
[3704]	2897	k = vmea(l)%k_soil(m)
[4400]	2898
[4671]	2899	match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
	2900	match_usm = surf_usm_h(0)%start_index(j,i) <= surf_usm_h(0)%end_index(j,i)
[4400]	2901
[3704]	2902	IF ( match_lsm ) THEN
[4671]	2903	mm = surf_lsm_h(0)%start_index(j,i)
	2904	vmea(l)%measured_vars_soil(m,n) = t_soil_h(0)%var_2d(k,mm)
[3704]	2905	ENDIF
[4400]	2906
[3704]	2907	IF ( match_usm ) THEN
[4671]	2908	mm = surf_usm_h(0)%start_index(j,i)
	2909	vmea(l)%measured_vars_soil(m,n) = t_wall_h(0)%val(k,mm)
[3704]	2910	ENDIF
	2911	ENDDO
[4400]	2912
	2913	CASE ( 'm_soil' ) ! soil moisture
[3704]	2914	DO m = 1, vmea(l)%ns_soil
[4400]	2915	j = MERGE( vmea(l)%j_soil(m), nys, vmea(l)%j_soil(m) > nys )
	2916	j = MERGE( j , nyn, j < nyn )
	2917	i = MERGE( vmea(l)%i_soil(m), nxl, vmea(l)%i_soil(m) > nxl )
	2918	i = MERGE( i , nxr, i < nxr )
[3704]	2919	k = vmea(l)%k_soil(m)
[4400]	2920
[4671]	2921	match_lsm = surf_lsm_h(0)%start_index(j,i) <= surf_lsm_h(0)%end_index(j,i)
[4400]	2922
[3704]	2923	IF ( match_lsm ) THEN
[4671]	2924	mm = surf_lsm_h(0)%start_index(j,i)
	2925	vmea(l)%measured_vars_soil(m,n) = m_soil_h(0)%var_2d(k,mm)
[3704]	2926	ENDIF
[4400]	2927
[3704]	2928	ENDDO
[4400]	2929
	2930	CASE ( 'ts' ) ! surface temperature
	2931	DO m = 1, vmea(l)%ns
[3522]	2932	!
[4498]	2933	!-- Surface data is only available on inner subdomains, not on ghost points. Hence,
	2934	!-- limit the indices.
[4400]	2935	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2936	j = MERGE( j , nyn, j < nyn )
	2937	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2938	i = MERGE( i , nxr, i < nxr )
	2939
[4498]	2940	DO mm = surf_def_h(0)%start_index(j,i), surf_def_h(0)%end_index(j,i)
[4400]	2941	vmea(l)%measured_vars(m,n) = surf_def_h(0)%pt_surface(mm)
	2942	ENDDO
[4671]	2943	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2944	vmea(l)%measured_vars(m,n) = surf_lsm_h(0)%pt_surface(mm)
[4400]	2945	ENDDO
[4671]	2946	DO mm = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
	2947	vmea(l)%measured_vars(m,n) = surf_usm_h(0)%pt_surface(mm)
[4400]	2948	ENDDO
	2949	ENDDO
	2950
	2951	CASE ( 'lwp' ) ! liquid water path
	2952	IF ( ASSOCIATED( ql ) ) THEN
	2953	DO m = 1, vmea(l)%ns
	2954	j = vmea(l)%j(m)
	2955	i = vmea(l)%i(m)
	2956
[4498]	2957	vmea(l)%measured_vars(m,n) = SUM( ql(nzb:nzt,j,i) * dzw(1:nzt+1) ) &
	2958	* rho_surface
[4400]	2959	ENDDO
	2960	ENDIF
	2961
	2962	CASE ( 'ps' ) ! surface pressure
	2963	vmea(l)%measured_vars(:,n) = surface_pressure
	2964
	2965	CASE ( 'pswrtg' ) ! platform speed above ground
	2966	vmea(l)%measured_vars(:,n) = 0.0_wp
	2967
	2968	CASE ( 'pswrta' ) ! platform speed in air
	2969	vmea(l)%measured_vars(:,n) = 0.0_wp
	2970
	2971	CASE ( 't_lw' ) ! water temperature
	2972	DO m = 1, vmea(l)%ns
	2973	!
	2974	!-- Surface data is only available on inner subdomains, not
	2975	!-- on ghost points. Hence, limit the indices.
	2976	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) > nys )
	2977	j = MERGE( j , nyn, j < nyn )
	2978	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) > nxl )
	2979	i = MERGE( i , nxr, i < nxr )
	2980
[4671]	2981	DO mm = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
	2982	IF ( surf_lsm_h(0)%water_surface(m) ) &
	2983	vmea(l)%measured_vars(m,n) = t_soil_h(0)%var_2d(nzt,m)
[4400]	2984	ENDDO
	2985
	2986	ENDDO
	2987	!
[3522]	2988	!-- More will follow ...
[4641]	2989	CASE ( 'ncaa' )
[3704]	2990	!
	2991	!-- No match found - just set a fill value
	2992	CASE DEFAULT
	2993	vmea(l)%measured_vars(:,n) = vmea(l)%fillout
[3522]	2994	END SELECT
	2995
[3494]	2996	ENDDO
[3434]	2997
	2998	ENDDO
[4400]	2999
[4438]	3000	CALL cpu_log( log_point_s(27), 'VM sampling', 'stop' )
	3001
[4498]	3002	END SUBROUTINE vm_sampling
[3434]	3003
[4400]	3004
[3471]	3005	END MODULE virtual_measurement_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |