Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |