Home

Context Navigation

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

Last change on this file since 4642 was 4642, checked in by suehring, 4 years ago
Do not set attribute bounds for time variable, as it refers to time_bounds which is not defined for non-aggregated quantities (according to data standard)
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 154.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |