Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |