Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |