Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |