Home

Context Navigation

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

Last change on this file since 4652 was 4645, checked in by suehring, 4 years ago
Bugfix in output of E_UTM_soil coordinate
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 154.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |