Home

Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

virtual_measurement_mod.f90

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

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

Context Navigation

source: palm/trunk/SOURCE/virtual_measurement_mod.f90

Download in other formats: