Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |