Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |