Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |