Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |