Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |