Home

Context Navigation

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

Last change on this file since 4406 was 4406, checked in by knoop, 4 years ago
Bugix: removed oro_rel wrong loop bounds and removed unnecessary restart method
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 164.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |