Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |