Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |