Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |