Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |