Home

Context Navigation

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

Last change on this file since 4403 was 4400, checked in by suehring, 5 years ago
Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 166.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |