Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |