Home

Context Navigation

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

Last change on this file since 4762 was 4752, checked in by suehring, 4 years ago
virtual measurements: Remove unnecessary output and merge output for quantities that are represented by the same variable in PALM, e.g. surface temperature and brightness temperature
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 157.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |