Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |