Home

Context Navigation

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

Last change on this file since 4377 was 4346, checked in by motisi, 5 years ago
Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 76.8 KB

Line
1	!> @virtual_measurement_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: virtual_measurement_mod.f90 4346 2019-12-18 11:55:56Z gronemeier $
27	! Introduction of wall_flags_total_0, which currently sets bits based on static
28	! topography information used in wall_flags_static_0
29	!
30	! 4329 2019-12-10 15:46:36Z motisi
31	! Renamed wall_flags_0 to wall_flags_static_0
32	!
33	! 4226 2019-09-10 17:03:24Z suehring
34	! Netcdf input routine for dimension length renamed
35	!
36	! 4182 2019-08-22 15:20:23Z scharf
37	! Corrected "Former revisions" section
38	!
39	! 4168 2019-08-16 13:50:17Z suehring
40	! Replace function get_topography_top_index by topo_top_ind
41	!
42	! 3988 2019-05-22 11:32:37Z kanani
43	! Add variables to enable steering of output interval for virtual measurements
44	!
45	! 3913 2019-04-17 15:12:28Z gronemeier
46	! Bugfix: rotate positions of measurements before writing them into file
47	!
48	! 3910 2019-04-17 11:46:56Z suehring
49	! Bugfix in rotation of UTM coordinates
50	!
51	! 3904 2019-04-16 18:22:51Z gronemeier
52	! Rotate coordinates of stations by given rotation_angle
53	!
54	! 3876 2019-04-08 18:41:49Z knoop
55	! Remove print statement
56	!
57	! 3854 2019-04-02 16:59:33Z suehring
58	! renamed nvar to nmeas, replaced USE chem_modules by USE chem_gasphase_mod and
59	! nspec by nvar
60	!
61	! 3766 2019-02-26 16:23:41Z raasch
62	! unused variables removed
63	!
64	! 3718 2019-02-06 11:08:28Z suehring
65	! Adjust variable name connections between UC2 and chemistry variables
66	!
67	! 3717 2019-02-05 17:21:16Z suehring
68	! Additional check + error numbers adjusted
69	!
70	! 3706 2019-01-29 20:02:26Z suehring
71	! unused variables removed
72	!
73	! 3705 2019-01-29 19:56:39Z suehring
74	! - initialization revised
75	! - binary data output
76	! - list of allowed variables extended
77	!
78	! 3704 2019-01-29 19:51:41Z suehring
79	! Sampling of variables
80	!
81	! 3473 2018-10-30 20:50:15Z suehring
82	! Initial revision
83	!
84	! Authors:
85	! --------
86	! @author Matthias Suehring
87	!
88	! Description:
89	! ------------
90	!> The module acts as an interface between 'real-world' observations and
91	!> model simulations. Virtual measurements will be taken in the model at the
92	!> coordinates representative for the 'real-world' observation coordinates.
93	!> More precisely, coordinates and measured quanties will be read from a
94	!> NetCDF file which contains all required information. In the model,
95	!> the same quantities (as long as all the required components are switched-on)
96	!> will be sampled at the respective positions and output into an extra file,
97	!> which allows for straight-forward comparison of model results with
98	!> observations.
99	!>
100	!> @todo list_of_allowed variables needs careful checking
101	!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
102	!> the (UC)2 standard
103	!> @note Fluxes are not processed
104	!------------------------------------------------------------------------------!
105	MODULE virtual_measurement_mod
106
107	USE arrays_3d, &
108	ONLY: q, pt, u, v, w, zu, zw
109
110	USE basic_constants_and_equations_mod, &
111	ONLY: pi
112
113	USE chem_gasphase_mod, &
114	ONLY: nvar
115
116	USE chem_modules, &
117	ONLY: chem_species
118
119	USE control_parameters, &
120	ONLY: air_chemistry, dz, humidity, io_blocks, io_group, neutral, &
121	message_string, time_since_reference_point, virtual_measurement
122
123	USE cpulog, &
124	ONLY: cpu_log, log_point
125
126	USE grid_variables, &
127	ONLY: dx, dy
128
129	USE indices, &
130	ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind, &
131	wall_flags_total_0
132
133	USE kinds
134
135	USE netcdf_data_input_mod, &
136	ONLY: init_model
137
138	USE pegrid
139
140	USE surface_mod, &
141	ONLY: surf_lsm_h, surf_usm_h
142
143	USE land_surface_model_mod, &
144	ONLY: nzb_soil, nzs, nzt_soil, zs, t_soil_h, m_soil_h
145
146	USE radiation_model_mod
147
148	USE urban_surface_mod, &
149	ONLY: nzb_wall, nzt_wall, t_wall_h
150
151
152	IMPLICIT NONE
153
154	TYPE virt_general
155	INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements
156	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
157	END TYPE virt_general
158
159	TYPE virt_mea
160
161	CHARACTER(LEN=100) :: feature_type !< type of the measurement
162	CHARACTER(LEN=100) :: filename_original !< name of the original file
163	CHARACTER(LEN=100) :: site !< name of the measurement site
164
165	CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables
166
167	INTEGER(iwp) :: ns = 0 !< number of observation coordinates on subdomain, for atmospheric measurements
168	INTEGER(iwp) :: ns_tot = 0 !< total number of observation coordinates, for atmospheric measurements
169	INTEGER(iwp) :: ntraj !< number of trajectories of a measurement
170	INTEGER(iwp) :: nmeas !< number of measured variables (atmosphere + soil)
171
172	INTEGER(iwp) :: ns_soil = 0 !< number of observation coordinates on subdomain, for soil measurements
173	INTEGER(iwp) :: ns_soil_tot = 0 !< total number of observation coordinates, for soil measurements
174
175	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues
176
177	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
178	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
179	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
180
181	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i_soil !< grid index for measurement position in x-direction
182	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j_soil !< grid index for measurement position in y-direction
183	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_soil !< grid index for measurement position in k-direction
184
185	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
186	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
187	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
188	LOGICAL :: soil_sampling = .FALSE. !< flag indicating that soil state variables were sampled
189
190	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
191	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
192	REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values
193	REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building
194	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
195	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
196
197	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_ag !< measurement height above ground level
198	REAL(wp), DIMENSION(:), ALLOCATABLE :: depth !< measurement depth in soil
199
200	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
201	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars_soil !< measured variables
202
203	END TYPE virt_mea
204
205	CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting
206	CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type
207
208	! This need to generalized
209	CHARACTER(LEN=8) :: char_filename = "filename" !< attribute name for filename
210	CHARACTER(LEN=11) :: char_soil = "soil_sample" !< attribute name for soil sampling indication
211	CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue
212	CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names
213	CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing
214	CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations
215	CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x
216	CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y
217	CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name
218	CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable
219	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
220	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
221	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
222
223	CHARACTER(LEN=6), DIMENSION(1:5) :: soil_vars = (/ & !< list of soil variables
224	't_soil', &
225	'm_soil', &
226	'lwc ', &
227	'lwcs ', &
228	'smp ' /)
229
230	CHARACTER(LEN=10), DIMENSION(0:1,1:8) :: chem_vars = RESHAPE( (/ &
231	'mcpm1 ', 'PM1 ', &
232	'mcpm2p5 ', 'PM2.5 ', &
233	'mcpm25 ', 'PM25 ', &
234	'mcpm10 ', 'PM10 ', &
235	'mfno2 ', 'NO2 ', &
236	'mfno ', 'NO ', &
237	'tro3 ', 'O3 ', &
238	'mfco ', 'CO ' &
239	/), (/ 2, 8 /) )
240	!
241	!-- MS: List requires careful revision!
242	CHARACTER(LEN=10), DIMENSION(1:54), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM
243	(/ 'hfls ', & ! surface latent heat flux (W/m2)
244	'hfss ', & ! surface sensible heat flux (W/m2)
245	'hur ', & ! relative humidity (-)
246	'hus ', & ! specific humidity (g/kg)
247	'haa ', & ! absolute atmospheric humidity (kg/m3)
248	'mcpm1 ', & ! mass concentration of PM1 (kg/m3)
249	'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3)
250	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
251	'mcco ', & ! mass concentration of CO (kg/m3)
252	'mcco2 ', & ! mass concentration of CO2 (kg/m3)
253	'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3)
254	'ncaa ', & ! number concentation of particles (1/m3)
255	'mfco ', & ! mole fraction of CO (mol/mol)
256	'mfco2 ', & ! mole fraction of CO2 (mol/mol)
257	'mfch4 ', & ! mole fraction of methane (mol/mol)
258	'mfnh3 ', & ! mole fraction of amonia (mol/mol)
259	'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol)
260	'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol)
261	'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol)
262	'mfh20 ', & ! mole fraction of water (mol/mol)
263	'plev ', & ! ? air pressure - hydrostaic + perturbation?
264	'rlds ', & ! surface downward longwave flux (W/m2)
265	'rlus ', & ! surface upward longwave flux (W/m2)
266	'rsds ', & ! surface downward shortwave flux (W/m2)
267	'rsus ', & ! surface upward shortwave flux (W/m2)
268	'ta ', & ! air temperature (degree C)
269	't_va ', & ! virtual accoustic temperature (K)
270	'theta ', & ! potential temperature (K)
271	'tro3 ', & ! mole fraction of ozone air (mol/mol)
272	'ts ', & ! scaling parameter of temperature (K)
273	'wspeed ', & ! ? wind speed - horizontal?
274	'wdir ', & ! wind direction
275	'us ', & ! friction velocity
276	'msoil ', & ! ? soil moisture - which depth?
277	'tsoil ', & ! ? soil temperature - which depth?
278	'u ', & ! u-component
279	'utheta ', & ! total eastward kinematic heat flux
280	'ua ', & ! eastward wind (is there any difference to u?)
281	'v ', & ! v-component
282	'vtheta ', & ! total northward kinematic heat flux
283	'va ', & ! northward wind (is there any difference to v?)
284	'w ', & ! w-component
285	'wtheta ', & ! total vertical kinematic heat flux
286	'rld ', & ! downward longwave radiative flux (W/m2)
287	'rlu ', & ! upnward longwave radiative flux (W/m2)
288	'rsd ', & ! downward shortwave radiative flux (W/m2)
289	'rsu ', & ! upward shortwave radiative flux (W/m2)
290	'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2)
291	'rnds ', & ! surface net downward radiative flux (W/m2)
292	't_soil ', &
293	'm_soil ', &
294	'lwc ', &
295	'lwcs ', &
296	'smp ' &
297	/)
298
299
300	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
301	LOGICAL :: init = .TRUE. !< flag indicating initialization of data output
302	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
303
304	REAL(wp) :: dt_virtual_measurement = 0.0_wp !< namelist parameter
305	REAL(wp) :: time_virtual_measurement = 0.0_wp !< time since last 3d output
306	REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start
307
308	TYPE( virt_general ) :: vmea_general !< data structure which encompass general variables
309	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure
310
311	INTERFACE vm_check_parameters
312	MODULE PROCEDURE vm_check_parameters
313	END INTERFACE vm_check_parameters
314
315	INTERFACE vm_data_output
316	MODULE PROCEDURE vm_data_output
317	END INTERFACE vm_data_output
318
319	INTERFACE vm_init
320	MODULE PROCEDURE vm_init
321	END INTERFACE vm_init
322
323	INTERFACE vm_last_actions
324	MODULE PROCEDURE vm_last_actions
325	END INTERFACE vm_last_actions
326
327	INTERFACE vm_parin
328	MODULE PROCEDURE vm_parin
329	END INTERFACE vm_parin
330
331	INTERFACE vm_sampling
332	MODULE PROCEDURE vm_sampling
333	END INTERFACE vm_sampling
334
335	SAVE
336
337	PRIVATE
338
339	!
340	!-- Public interfaces
341	PUBLIC vm_check_parameters, vm_data_output, vm_init, vm_last_actions, &
342	vm_parin, vm_sampling
343
344	!
345	!-- Public variables
346	PUBLIC dt_virtual_measurement, time_virtual_measurement, vmea, vmea_general, vm_time_start
347
348	CONTAINS
349
350
351	!------------------------------------------------------------------------------!
352	! Description:
353	! ------------
354	!> Check parameters for virtual measurement module
355	!------------------------------------------------------------------------------!
356	SUBROUTINE vm_check_parameters
357
358	USE control_parameters, &
359	ONLY: message_string, virtual_measurement
360
361	USE netcdf_data_input_mod, &
362	ONLY: input_pids_static, input_pids_vm
363
364	IMPLICIT NONE
365
366	!
367	!-- Virtual measurements require a setup file.
368	IF ( virtual_measurement .AND. .NOT. input_pids_vm ) THEN
369	message_string = 'If virtual measurements are taken, a setup input ' // &
370	'file for the site locations is mandatory.'
371	CALL message( 'vm_check_parameters', 'PA0533', 1, 2, 0, 6, 0 )
372	ENDIF
373	!
374	!-- In case virtual measurements are taken, a static input file is required.
375	!-- This is because UTM coordinates for the PALM domain origin are required
376	!-- for correct mapping of the measurements.
377	!-- ToDo: Revise this later and remove this requirement.
378	IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN
379	message_string = 'If virtual measurements are taken, a static input ' //&
380	'file is mandatory.'
381	CALL message( 'vm_check_parameters', 'PA0534', 1, 2, 0, 6, 0 )
382	ENDIF
383
384	END SUBROUTINE vm_check_parameters
385
386	!------------------------------------------------------------------------------!
387	! Description:
388	! ------------
389	!> Read namelist for the virtual measurement module
390	!------------------------------------------------------------------------------!
391	SUBROUTINE vm_parin
392
393	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
394
395	NAMELIST /virtual_measurement_parameters/ dt_virtual_measurement, &
396	use_virtual_measurement, &
397	vm_time_start
398
399	line = ' '
400
401	!
402	!-- Try to find stg package
403	REWIND ( 11 )
404	line = ' '
405	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
406	READ ( 11, '(A)', END=20 ) line
407	ENDDO
408	BACKSPACE ( 11 )
409
410	!
411	!-- Read namelist
412	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
413
414	!
415	!-- Set flag that indicates that the virtual measurement module is switched on
416	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
417
418	GOTO 20
419
420	10 BACKSPACE( 11 )
421	READ( 11 , '(A)') line
422	CALL parin_fail_message( 'virtual_measurement_parameters', line )
423
424	20 CONTINUE
425
426	END SUBROUTINE vm_parin
427
428
429	!------------------------------------------------------------------------------!
430	! Description:
431	! ------------
432	!> Initialize virtual measurements: read coordiante arrays and measured
433	!> variables, set indicies indicating the measurement points, read further
434	!> attributes, etc..
435	!------------------------------------------------------------------------------!
436	SUBROUTINE vm_init
437
438	USE arrays_3d, &
439	ONLY: zu, zw
440
441	USE grid_variables, &
442	ONLY: ddx, ddy, dx, dy
443
444	USE indices, &
445	ONLY: nxl, nxr, nyn, nys
446
447	USE netcdf_data_input_mod, &
448	ONLY: get_dimension_length, &
449	init_model, &
450	input_file_vm, &
451	netcdf_data_input_att, &
452	netcdf_data_input_var
453
454	IMPLICIT NONE
455
456	CHARACTER(LEN=5) :: dum !< dummy string indicate station id
457	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
458	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
459
460	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
461	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
462	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
463	INTEGER(iwp) :: j !< grid index of observation point in x-direction
464	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
465	INTEGER(iwp) :: k !< grid index of observation point in x-direction
466	INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate
467	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
468	INTEGER(iwp) :: ksurf !< topography top index
469	INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate
470	INTEGER(iwp) :: l !< running index over all stations
471	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
472	INTEGER(iwp) :: ll !< running index over all measured variables in file
473	INTEGER(iwp) :: lll !< running index over all allowed variables
474	INTEGER(iwp) :: n !< running index over trajectory coordinates
475	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
476	INTEGER(iwp) :: t !< running index over number of trajectories
477	INTEGER(iwp) :: m
478
479	INTEGER(KIND=1):: soil_dum
480
481	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ns_all !< dummy array used to sum-up the number of observation coordinates
482
483	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions
484
485	! LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism
486	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
487
488	REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM
489	REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM
490	REAL(wp) :: fill_zag !< _FillValue for height coordinate
491
492	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable
493	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable
494	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm_tmp !< EUTM coordinate before rotation
495	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm_tmp !< NUTM coordinate before rotation
496	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable
497	!
498	!-- Obtain number of sites. Also, pass the 'open' string, in order to initially
499	!-- open the measurement driver.
500	CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, &
501	vmea_general%id_vm, input_file_vm, &
502	global_attribute, 'open', '' )
503	!
504	!-- Allocate data structure which encompass all required information, such as
505	!-- grid points indicies, absolute UTM coordinates, the measured quantities,
506	!-- etc. .
507	ALLOCATE( vmea(1:vmea_general%nvm) )
508	!
509	!-- Allocate flag array. This dummy array is used to identify grid points
510	!-- where virtual measurements should be taken. Please note, at least one
511	!-- ghost point is required, in order to include also the surrounding
512	!-- grid points of the original coordinate.
513	ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
514	meas_flag = 0
515	!
516	!-- Loop over all sites.
517	DO l = 1, vmea_general%nvm
518	!
519	!-- Determine suffix which contains the ID, ordered according to the number
520	!-- of measurements.
521	IF( l < 10 ) THEN
522	WRITE( dum, '(I1)') l
523	ELSEIF( l < 100 ) THEN
524	WRITE( dum, '(I2)') l
525	ELSEIF( l < 1000 ) THEN
526	WRITE( dum, '(I3)') l
527	ELSEIF( l < 10000 ) THEN
528	WRITE( dum, '(I4)') l
529	ELSEIF( l < 100000 ) THEN
530	WRITE( dum, '(I5)') l
531	ENDIF
532	!
533	!-- Read site coordinates (UTM).
534	CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx // &
535	TRIM( dum ), vmea_general%id_vm, '', &
536	global_attribute, '', '' )
537	CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy // &
538	TRIM( dum ), vmea_general%id_vm, '', &
539	global_attribute, '', '' )
540	!
541	!-- Read site name
542	CALL netcdf_data_input_att( vmea(l)%site, char_site // TRIM( dum ), &
543	vmea_general%id_vm, '', global_attribute, &
544	'', '' )
545	!
546	!-- Read type of the measurement (trajectory, profile, timeseries).
547	CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature // &
548	TRIM( dum ), vmea_general%id_vm, '', &
549	global_attribute, '', '' )
550	!
551	!-- Read the name of the original file where observational data is stored.
552	CALL netcdf_data_input_att( vmea(l)%filename_original, char_filename // &
553	TRIM( dum ), vmea_general%id_vm, '', &
554	global_attribute, '', '' )
555	!
556	!-- Read a flag which indicates that also soil quantities are take at the
557	!-- respective site (is part of the virtual measurement driver).
558	CALL netcdf_data_input_att( soil_dum, char_soil // TRIM( dum ), &
559	vmea_general%id_vm, '', global_attribute, &
560	'', '' )
561	!
562	!-- Set flag for soil-sampling.
563	IF ( soil_dum == 1 ) vmea(l)%soil_sampling = .TRUE.
564	!
565	!--- Set logicals depending on the type of the measurement
566	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
567	vmea(l)%timseries_profile = .TRUE.
568	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
569	vmea(l)%timseries = .TRUE.
570	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
571	vmea(l)%trajectory = .TRUE.
572	!
573	!-- Give error message in case the type matches non of the pre-defined types.
574	ELSE
575	message_string = 'Attribue featureType = ' // &
576	TRIM( vmea(l)%feature_type ) // &
577	' is not allowed.'
578	CALL message( 'vm_init', 'PA0535', 1, 2, 0, 6, 0 )
579	ENDIF
580	!
581	!-- Read string with all measured variables at this site
582	measured_variables_file = ''
583	CALL netcdf_data_input_var( measured_variables_file, &
584	char_mv // TRIM( dum ), vmea_general%id_vm )
585	!
586	!-- Count the number of measured variables. Only count variables that match
587	!-- with the allowed variables.
588	!-- Please note, for some NetCDF interal reasons characters end with a NULL,
589	!-- i.e. also empty characters contain a NULL. Therefore, check the strings
590	!-- for a NULL to get the correct character length in order to compare
591	!-- them with the list of allowed variables.
592	vmea(l)%nmeas = 0
593	DO ll = 1, SIZE( measured_variables_file )
594	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
595	measured_variables_file(ll)(1:1) /= ' ') THEN
596	!
597	!-- Obtain character length of the character
598	len_char = 1
599	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
600	.AND. measured_variables_file(ll)(len_char:len_char) /= ' ' )
601	len_char = len_char + 1
602	ENDDO
603	len_char = len_char - 1
604	!
605	!-- Now, compare the measured variable with the list of allowed
606	!-- variables.
607	DO lll= 1, SIZE( list_allowed_variables )
608	IF ( measured_variables_file(ll)(1:len_char) == &
609	TRIM( list_allowed_variables(lll) ) ) THEN
610	vmea(l)%nmeas = vmea(l)%nmeas + 1
611	measured_variables(vmea(l)%nmeas) = &
612	measured_variables_file(ll)(1:len_char)
613	ENDIF
614	ENDDO
615	ENDIF
616	ENDDO
617
618
619	!
620	!-- Allocate array for the measured variables names for the respective site.
621	ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nmeas) )
622
623	DO ll = 1, vmea(l)%nmeas
624	vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
625	ENDDO
626	!
627	!-- In case of chemistry, check if species is considered in the modelled
628	!-- chemistry mechanism.
629	! IF ( air_chemistry ) THEN
630	! DO ll = 1, vmea(l)%nmeas
631	! chem_include = .FALSE.
632	! DO n = 1, nvar
633	! IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == &
634	! TRIM( chem_species(n)%name ) ) chem_include = .TRUE.
635	! ENDDO
636	! !
637	! !-- Revise this. It should only check for chemistry variables and not for all!
638	! IF ( .NOT. chem_include ) THEN
639	! message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // &
640	! ' is not considered in the modelled ' // &
641	! 'chemistry mechanism'
642	! CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
643	! ENDIF
644	! ENDDO
645	! ENDIF
646	!
647	!-- Read the UTM coordinates for the actual site. Based on the coordinates,
648	!-- define the grid-index space on each subdomain where virtual measurements
649	!-- should be taken. Note, the entire coordinate arrays will not be stored
650	!-- as this would exceed memory requirements, particularly for trajectory
651	!-- measurements.
652	IF ( vmea(l)%nmeas > 0 ) THEN
653	!
654	!-- For stationary measurements UTM coordinates are just one value and
655	!-- its dimension is "station", while for mobile measurements UTM
656	!-- coordinates are arrays depending on the number of trajectories and
657	!-- time, according to (UC)2 standard. First, inquire dimension length
658	!-- of the UTM coordinates.
659	IF ( vmea(l)%trajectory ) THEN
660	!
661	!-- For non-stationary measurements read the number of trajectories
662	!-- and the number of time coordinates.
663	CALL get_dimension_length( vmea_general%id_vm, &
664	vmea(l)%ntraj, &
665	"traj" // &
666	TRIM( dum ) )
667	CALL get_dimension_length( vmea_general%id_vm, &
668	dim_ntime, &
669	"ntime" // &
670	TRIM( dum ) )
671	!
672	!-- For stationary measurements the dimension for UTM and time
673	!-- coordinates is 1.
674	ELSE
675	vmea(l)%ntraj = 1
676	dim_ntime = 1
677	ENDIF
678	!
679	!- Allocate array which defines individual time frame for each
680	!-- trajectory or station.
681	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
682	!
683	!-- Allocate temporary arrays for UTM and height coordinates. Note,
684	!-- on file UTM coordinates might be 1D or 2D variables
685	ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
686	ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
687	ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) )
688
689	ALLOCATE( e_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
690	ALLOCATE( n_utm_tmp(1:vmea(l)%ntraj,1:dim_ntime) )
691	!
692	!-- Read _FillValue attributes of the coordinate dimensions.
693	CALL netcdf_data_input_att( fill_eutm, char_fillvalue, &
694	vmea_general%id_vm, '', &
695	.NOT. global_attribute, '', &
696	char_eutm // TRIM( dum ) )
697	CALL netcdf_data_input_att( fill_nutm, char_fillvalue, &
698	vmea_general%id_vm, '', &
699	.NOT. global_attribute, '', &
700	char_nutm // TRIM( dum ) )
701	CALL netcdf_data_input_att( fill_zag, char_fillvalue, &
702	vmea_general%id_vm, '', &
703	.NOT. global_attribute, '', &
704	char_zag // TRIM( dum ) )
705	!
706	!-- Read UTM and height coordinates coordinates for all trajectories and
707	!-- times.
708	IF ( vmea(l)%trajectory ) THEN
709	CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), &
710	vmea_general%id_vm, &
711	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
712	CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), &
713	vmea_general%id_vm, &
714	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
715	CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), &
716	vmea_general%id_vm, &
717	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
718	ELSE
719	CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), &
720	vmea_general%id_vm )
721	CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), &
722	vmea_general%id_vm )
723	CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), &
724	vmea_general%id_vm )
725	ENDIF
726	!
727	!-- Based on UTM coordinates, check if the measurement station or parts
728	!-- of the trajectory is on subdomain. This case, setup grid index space
729	!-- sample these quantities.
730	meas_flag = 0
731	DO t = 1, vmea(l)%ntraj
732	!
733	!-- First, compute relative x- and y-coordinates with respect to the
734	!-- lower-left origin of the model domain, which is the difference
735	!-- between UTM coordinates. Note, if the origin is not correct, the
736	!-- virtual sites will be misplaced.
737	e_utm_tmp(t,1:dim_ntime) = e_utm(t,1:dim_ntime) - init_model%origin_x
738	n_utm_tmp(t,1:dim_ntime) = n_utm(t,1:dim_ntime) - init_model%origin_y
739	e_utm(t,1:dim_ntime) = COS( init_model%rotation_angle * pi / 180.0_wp ) &
740	* e_utm_tmp(t,1:dim_ntime) &
741	- SIN( init_model%rotation_angle * pi / 180.0_wp ) &
742	* n_utm_tmp(t,1:dim_ntime)
743	n_utm(t,1:dim_ntime) = SIN( init_model%rotation_angle * pi / 180.0_wp ) &
744	* e_utm_tmp(t,1:dim_ntime) &
745	+ COS( init_model%rotation_angle * pi / 180.0_wp ) &
746	* n_utm_tmp(t,1:dim_ntime)
747	!
748	!-- Determine the individual time coordinate length for each station and
749	!-- trajectory. This is required as several stations and trajectories
750	!-- are merged into one file but they do not have the same number of
751	!-- points in time, hence, missing values may occur and cannot be
752	!-- processed further. This is actually a work-around for the specific
753	!-- (UC)2 dataset, but it won't harm in anyway.
754	vmea(l)%dim_t(t) = 0
755	DO n = 1, dim_ntime
756	IF ( e_utm(t,n) /= fill_eutm .AND. &
757	n_utm(t,n) /= fill_nutm .AND. &
758	z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n
759	ENDDO
760	!
761	!-- Compute grid indices relative to origin and check if these are
762	!-- on the subdomain. Note, virtual measurements will be taken also
763	!-- at grid points surrounding the station, hence, check also for
764	!-- these grid points.
765	DO n = 1, vmea(l)%dim_t(t)
766	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
767	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
768	!
769	!-- Is the observation point on subdomain?
770	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
771	js >= nys .AND. js <= nyn )
772	!
773	!-- Check if observation coordinate is on subdomain
774	IF ( on_pe ) THEN
775	!
776	!-- Determine vertical index which correspond to the observation
777	!-- height.
778	ksurf = topo_top_ind(js,is,0)
779	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
780	!
781	!-- Set mask array at the observation coordinates. Also, flag the
782	!-- surrounding coordinate points, but first check whether the
783	!-- surrounding coordinate points are on the subdomain.
784	kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 >= ksurf )
785	ku = MERGE( ks+1, ks, ks+1 < nzt+1 )
786
787	DO i = is-1, is+1
788	DO j = js-1, js+1
789	DO k = kl, ku
790	meas_flag(k,j,i) = MERGE( &
791	IBSET( meas_flag(k,j,i), 0 ), &
792	0, &
793	BTEST( wall_flags_total_0(k,j,i), 0 ) &
794	)
795	ENDDO
796	ENDDO
797	ENDDO
798	ENDIF
799	ENDDO
800
801	ENDDO
802	!
803	!-- Based on the flag array count the number of of sampling coordinates.
804	!-- Please note, sampling coordinates in atmosphere and soil may be
805	!-- different, as within the soil all levels will be measured.
806	!-- Hence, count individually. Start with atmoshere.
807	ns = 0
808	DO i = nxl-1, nxr+1
809	DO j = nys-1, nyn+1
810	DO k = nzb, nzt+1
811	ns = ns + MERGE( 1, 0, BTEST( meas_flag(k,j,i), 0 ) )
812	ENDDO
813	ENDDO
814	ENDDO
815
816	!
817	!-- Store number of observation points on subdomain and allocate index
818	!-- arrays as well as array containing height information.
819	vmea(l)%ns = ns
820
821	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
822	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
823	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
824	ALLOCATE( vmea(l)%z_ag(1:vmea(l)%ns) )
825	!
826	!-- Based on the flag array store the grid indices which correspond to
827	!-- the observation coordinates.
828	ns = 0
829	DO i = nxl-1, nxr+1
830	DO j = nys-1, nyn+1
831	DO k = nzb, nzt+1
832	IF ( BTEST( meas_flag(k,j,i), 0 ) ) THEN
833	ns = ns + 1
834	vmea(l)%i(ns) = i
835	vmea(l)%j(ns) = j
836	vmea(l)%k(ns) = k
837	vmea(l)%z_ag(ns) = zu(k) - zw(topo_top_ind(j,i,0))
838	ENDIF
839	ENDDO
840	ENDDO
841	ENDDO
842	!
843	!-- Same for the soil. Based on the flag array, count the number of
844	!-- sampling coordinates in soil. Sample at all soil levels in this case.
845	IF ( vmea(l)%soil_sampling ) THEN
846	DO i = nxl, nxr
847	DO j = nys, nyn
848	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
849	IF ( surf_lsm_h%start_index(j,i) <= &
850	surf_lsm_h%end_index(j,i) ) THEN
851	vmea(l)%ns_soil = vmea(l)%ns_soil + &
852	nzt_soil - nzb_soil + 1
853	ENDIF
854	IF ( surf_usm_h%start_index(j,i) <= &
855	surf_usm_h%end_index(j,i) ) THEN
856	vmea(l)%ns_soil = vmea(l)%ns_soil + &
857	nzt_wall - nzb_wall + 1
858	ENDIF
859	ENDIF
860	ENDDO
861	ENDDO
862	ENDIF
863	!
864	!-- Allocate index arrays as well as array containing height information
865	!-- for soil.
866	IF ( vmea(l)%soil_sampling ) THEN
867	ALLOCATE( vmea(l)%i_soil(1:vmea(l)%ns_soil) )
868	ALLOCATE( vmea(l)%j_soil(1:vmea(l)%ns_soil) )
869	ALLOCATE( vmea(l)%k_soil(1:vmea(l)%ns_soil) )
870	ALLOCATE( vmea(l)%depth(1:vmea(l)%ns_soil) )
871	ENDIF
872	!
873	!-- For soil, store the grid indices.
874	ns = 0
875	IF ( vmea(l)%soil_sampling ) THEN
876	DO i = nxl, nxr
877	DO j = nys, nyn
878	IF ( ANY( BTEST( meas_flag(:,j,i), 0 ) ) ) THEN
879	IF ( surf_lsm_h%start_index(j,i) <= &
880	surf_lsm_h%end_index(j,i) ) THEN
881	m = surf_lsm_h%start_index(j,i)
882	DO k = nzb_soil, nzt_soil
883	ns = ns + 1
884	vmea(l)%i_soil(ns) = i
885	vmea(l)%j_soil(ns) = j
886	vmea(l)%k_soil(ns) = k
887	vmea(l)%depth(ns) = zs(k)
888	ENDDO
889	ENDIF
890
891	IF ( surf_usm_h%start_index(j,i) <= &
892	surf_usm_h%end_index(j,i) ) THEN
893	m = surf_usm_h%start_index(j,i)
894	DO k = nzb_wall, nzt_wall
895	ns = ns + 1
896	vmea(l)%i_soil(ns) = i
897	vmea(l)%j_soil(ns) = j
898	vmea(l)%k_soil(ns) = k
899	vmea(l)%depth(ns) = surf_usm_h%zw(k,m)
900	ENDDO
901	ENDIF
902	ENDIF
903	ENDDO
904	ENDDO
905	ENDIF
906	!
907	!-- Allocate array to save the sampled values.
908	ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) )
909
910	IF ( vmea(l)%soil_sampling ) &
911	ALLOCATE( vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil, &
912	1:vmea(l)%nmeas) )
913	!
914	!-- Initialize with _FillValues
915	vmea(l)%measured_vars(1:vmea(l)%ns,1:vmea(l)%nmeas) = vmea(l)%fillout
916	IF ( vmea(l)%soil_sampling ) &
917	vmea(l)%measured_vars_soil(1:vmea(l)%ns_soil,1:vmea(l)%nmeas) = &
918	vmea(l)%fillout
919	!
920	!-- Deallocate temporary coordinate arrays
921	IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm )
922	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
923	IF ( ALLOCATED( e_utm_tmp ) ) DEALLOCATE( e_utm_tmp )
924	IF ( ALLOCATED( n_utm_tmp ) ) DEALLOCATE( n_utm_tmp )
925	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
926	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag )
927	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( vmea(l)%dim_t )
928	ENDIF
929	ENDDO
930	!
931	!-- Close input file for virtual measurements. Therefore, just call
932	!-- the read attribute routine with the "close" option.
933	CALL netcdf_data_input_att( vmea_general%nvm, char_numstations, &
934	vmea_general%id_vm, '', &
935	global_attribute, 'close', '' )
936	!
937	!-- Sum-up the number of observation coordiates, for atmosphere first.
938	!-- This is actually only required for data output.
939	ALLOCATE( ns_all(1:vmea_general%nvm) )
940	ns_all = 0
941	#if defined( __parallel )
942	CALL MPI_ALLREDUCE( vmea(:)%ns, ns_all(:), vmea_general%nvm, MPI_INTEGER, &
943	MPI_SUM, comm2d, ierr )
944	#else
945	ns_all(:) = vmea(:)%ns
946	#endif
947	vmea(:)%ns_tot = ns_all(:)
948	!
949	!-- Now for soil
950	ns_all = 0
951	#if defined( __parallel )
952	CALL MPI_ALLREDUCE( vmea(:)%ns_soil, ns_all(:), vmea_general%nvm, &
953	MPI_INTEGER, MPI_SUM, comm2d, ierr )
954	#else
955	ns_all(:) = vmea(:)%ns_soil
956	#endif
957	vmea(:)%ns_soil_tot = ns_all(:)
958
959	DEALLOCATE( ns_all )
960	!
961	!-- Dellocate flag array
962	DEALLOCATE( meas_flag )
963	!
964	!-- Initialize binary data output of virtual measurements.
965	!-- Open binary output file.
966	CALL check_open( 27 )
967	!
968	!-- Output header information.
969	CALL vm_data_output
970
971	END SUBROUTINE vm_init
972
973
974	!------------------------------------------------------------------------------!
975	! Description:
976	! ------------
977	!> Binary data output.
978	!------------------------------------------------------------------------------!
979	SUBROUTINE vm_data_output
980
981	USE pegrid
982
983	IMPLICIT NONE
984
985	INTEGER(iwp) :: i !< running index over IO blocks
986	INTEGER(iwp) :: l !< running index over all stations
987	INTEGER(iwp) :: n !< running index over all measured variables at a station
988	!
989	!-- Header output on each PE
990	IF ( init ) THEN
991
992	DO i = 0, io_blocks-1
993	IF ( i == io_group ) THEN
994	WRITE ( 27 ) 'number of measurements '
995	WRITE ( 27 ) vmea_general%nvm
996
997	DO l = 1, vmea_general%nvm
998	WRITE ( 27 ) 'site '
999	WRITE ( 27 ) vmea(l)%site
1000	WRITE ( 27 ) 'file '
1001	WRITE ( 27 ) vmea(l)%filename_original
1002	WRITE ( 27 ) 'feature_type '
1003	WRITE ( 27 ) vmea(l)%feature_type
1004	WRITE ( 27 ) 'origin_x_obs '
1005	WRITE ( 27 ) vmea(l)%origin_x_obs
1006	WRITE ( 27 ) 'origin_y_obs '
1007	WRITE ( 27 ) vmea(l)%origin_y_obs
1008	WRITE ( 27 ) 'total number of observation points'
1009	WRITE ( 27 ) vmea(l)%ns_tot
1010	WRITE ( 27 ) 'number of measured variables '
1011	WRITE ( 27 ) vmea(l)%nmeas
1012	WRITE ( 27 ) 'variables '
1013	WRITE ( 27 ) vmea(l)%measured_vars_name(:)
1014	WRITE ( 27 ) 'number of observation points '
1015	WRITE ( 27 ) vmea(l)%ns
1016	WRITE ( 27 ) 'E_UTM '
1017	WRITE ( 27 ) init_model%origin_x &
1018	+ REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
1019	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
1020	+ REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
1021	* SIN( init_model%rotation_angle * pi / 180.0_wp )
1022	WRITE ( 27 ) 'N_UTM '
1023	WRITE ( 27 ) init_model%origin_y &
1024	- REAL( vmea(l)%i(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dx &
1025	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
1026	+ REAL( vmea(l)%j(1:vmea(l)%ns) + 0.5_wp, KIND = wp ) * dy &
1027	* COS( init_model%rotation_angle * pi / 180.0_wp )
1028	WRITE ( 27 ) 'Z_AG '
1029	WRITE ( 27 ) vmea(l)%z_ag(1:vmea(l)%ns)
1030	WRITE ( 27 ) 'soil sampling '
1031	WRITE ( 27 ) MERGE( 'yes ', &
1032	'no ', &
1033	vmea(l)%soil_sampling )
1034
1035	IF ( vmea(l)%soil_sampling ) THEN
1036	WRITE ( 27 ) 'total number of soil points '
1037	WRITE ( 27 ) vmea(l)%ns_soil_tot
1038	WRITE ( 27 ) 'number of soil points '
1039	WRITE ( 27 ) vmea(l)%ns_soil
1040	WRITE ( 27 ) 'E_UTM soil '
1041	WRITE ( 27 ) init_model%origin_x &
1042	+ REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
1043	KIND = wp ) * dx &
1044	* COS( init_model%rotation_angle * pi / 180.0_wp ) &
1045	+ REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
1046	KIND = wp ) * dy &
1047	* SIN( init_model%rotation_angle * pi / 180.0_wp )
1048	WRITE ( 27 ) 'N_UTM soil '
1049	WRITE ( 27 ) init_model%origin_y &
1050	- REAL( vmea(l)%i_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
1051	KIND = wp ) * dx &
1052	* SIN( init_model%rotation_angle * pi / 180.0_wp ) &
1053	+ REAL( vmea(l)%j_soil(1:vmea(l)%ns_soil) + 0.5_wp, &
1054	KIND = wp ) * dy &
1055	* COS( init_model%rotation_angle * pi / 180.0_wp )
1056	WRITE ( 27 ) 'DEPTH '
1057	WRITE ( 27 ) vmea(l)%depth(1:vmea(l)%ns_soil)
1058	ENDIF
1059	ENDDO
1060
1061	ENDIF
1062	ENDDO
1063
1064	#if defined( __parallel )
1065	CALL MPI_BARRIER( comm2d, ierr )
1066	#endif
1067	!
1068	!-- After header information is written, set control flag to .FALSE.
1069	init = .FALSE.
1070	!
1071	!-- Data output at each measurement timestep on each PE
1072	ELSE
1073	DO i = 0, io_blocks-1
1074
1075	IF ( i == io_group ) THEN
1076	WRITE( 27 ) 'output time '
1077	WRITE( 27 ) time_since_reference_point
1078	DO l = 1, vmea_general%nvm
1079	!
1080	!-- Skip binary writing if no observation points are defined on PE
1081	IF ( vmea(l)%ns < 1 .AND. vmea(l)%ns_soil < 1) CYCLE
1082	DO n = 1, vmea(l)%nmeas
1083	WRITE( 27 ) vmea(l)%measured_vars_name(n)
1084	IF ( vmea(l)%soil_sampling .AND. &
1085	ANY( TRIM( vmea(l)%measured_vars_name(n)) == &
1086	soil_vars ) ) THEN
1087	WRITE( 27 ) vmea(l)%measured_vars_soil(:,n)
1088	ELSE
1089	WRITE( 27 ) vmea(l)%measured_vars(:,n)
1090	ENDIF
1091	ENDDO
1092
1093	ENDDO
1094	ENDIF
1095	ENDDO
1096	#if defined( __parallel )
1097	CALL MPI_BARRIER( comm2d, ierr )
1098	#endif
1099	ENDIF
1100
1101	END SUBROUTINE vm_data_output
1102
1103
1104	!------------------------------------------------------------------------------!
1105	! Description:
1106	! ------------
1107	!> Write end-of-file statement as last action.
1108	!------------------------------------------------------------------------------!
1109	SUBROUTINE vm_last_actions
1110
1111	USE pegrid
1112
1113	IMPLICIT NONE
1114
1115	INTEGER(iwp) :: i !< running index over IO blocks
1116
1117	DO i = 0, io_blocks-1
1118	IF ( i == io_group ) THEN
1119	WRITE( 27 ) 'EOF '
1120	ENDIF
1121	ENDDO
1122	#if defined( __parallel )
1123	CALL MPI_BARRIER( comm2d, ierr )
1124	#endif
1125	!
1126	!-- Close binary file
1127	CALL close_file( 27 )
1128
1129	END SUBROUTINE vm_last_actions
1130
1131	!------------------------------------------------------------------------------!
1132	! Description:
1133	! ------------
1134	!> Sampling of the actual quantities along the observation coordinates
1135	!------------------------------------------------------------------------------!
1136	SUBROUTINE vm_sampling
1137
1138	USE arrays_3d, &
1139	ONLY: exner, pt, q, u, v, w
1140
1141	USE radiation_model_mod, &
1142	ONLY: radiation
1143
1144	USE surface_mod, &
1145	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
1146
1147	IMPLICIT NONE
1148
1149	INTEGER(iwp) :: i !< grid index in x-direction
1150	INTEGER(iwp) :: j !< grid index in y-direction
1151	INTEGER(iwp) :: k !< grid index in z-direction
1152	INTEGER(iwp) :: ind_chem !< dummy index to identify chemistry variable and translate it from (UC)2 standard to interal naming
1153	INTEGER(iwp) :: l !< running index over the number of stations
1154	INTEGER(iwp) :: m !< running index over all virtual observation coordinates
1155	INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate
1156	INTEGER(iwp) :: n !< running index over all measured variables at a station
1157	INTEGER(iwp) :: nn !< running index over the number of chemcal species
1158
1159	LOGICAL :: match_lsm !< flag indicating natural-type surface
1160	LOGICAL :: match_usm !< flag indicating urban-type surface
1161	!
1162	!-- Loop over all sites.
1163	DO l = 1, vmea_general%nvm
1164	!
1165	!-- At the beginning, set _FillValues
1166	IF ( ALLOCATED( vmea(l)%measured_vars ) ) &
1167	vmea(l)%measured_vars = vmea(l)%fillout
1168	IF ( ALLOCATED( vmea(l)%measured_vars_soil ) ) &
1169	vmea(l)%measured_vars_soil = vmea(l)%fillout
1170	!
1171	!-- Loop over all variables measured at this site.
1172	DO n = 1, vmea(l)%nmeas
1173
1174	SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
1175
1176	CASE ( 'theta' )
1177	IF ( .NOT. neutral ) THEN
1178	DO m = 1, vmea(l)%ns
1179	k = vmea(l)%k(m)
1180	j = vmea(l)%j(m)
1181	i = vmea(l)%i(m)
1182	vmea(l)%measured_vars(m,n) = pt(k,j,i)
1183	ENDDO
1184	ENDIF
1185
1186	CASE ( 'ta' )
1187	IF ( .NOT. neutral ) THEN
1188	DO m = 1, vmea(l)%ns
1189	k = vmea(l)%k(m)
1190	j = vmea(l)%j(m)
1191	i = vmea(l)%i(m)
1192	vmea(l)%measured_vars(m,n) = pt(k,j,i) * exner( k )
1193	ENDDO
1194	ENDIF
1195
1196	CASE ( 't_va' )
1197
1198	CASE ( 'hus', 'haa' )
1199	IF ( humidity ) THEN
1200	DO m = 1, vmea(l)%ns
1201	k = vmea(l)%k(m)
1202	j = vmea(l)%j(m)
1203	i = vmea(l)%i(m)
1204	vmea(l)%measured_vars(m,n) = q(k,j,i)
1205	ENDDO
1206	ENDIF
1207
1208	CASE ( 'u', 'ua' )
1209	DO m = 1, vmea(l)%ns
1210	k = vmea(l)%k(m)
1211	j = vmea(l)%j(m)
1212	i = vmea(l)%i(m)
1213	vmea(l)%measured_vars(m,n) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1214	ENDDO
1215
1216	CASE ( 'v', 'va' )
1217	DO m = 1, vmea(l)%ns
1218	k = vmea(l)%k(m)
1219	j = vmea(l)%j(m)
1220	i = vmea(l)%i(m)
1221	vmea(l)%measured_vars(m,n) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1222	ENDDO
1223
1224	CASE ( 'w' )
1225	DO m = 1, vmea(l)%ns
1226	k = MAX ( 1, vmea(l)%k(m) )
1227	j = vmea(l)%j(m)
1228	i = vmea(l)%i(m)
1229	vmea(l)%measured_vars(m,n) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1230	ENDDO
1231
1232	CASE ( 'wspeed' )
1233	DO m = 1, vmea(l)%ns
1234	k = vmea(l)%k(m)
1235	j = vmea(l)%j(m)
1236	i = vmea(l)%i(m)
1237	vmea(l)%measured_vars(m,n) = SQRT( &
1238	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
1239	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
1240	)
1241	ENDDO
1242
1243	CASE ( 'wdir' )
1244	DO m = 1, vmea(l)%ns
1245	k = vmea(l)%k(m)
1246	j = vmea(l)%j(m)
1247	i = vmea(l)%i(m)
1248
1249	vmea(l)%measured_vars(m,n) = ATAN2( &
1250	- 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), &
1251	- 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) &
1252	) * 180.0_wp / pi
1253	ENDDO
1254
1255	CASE ( 'utheta' )
1256	DO m = 1, vmea(l)%ns
1257	k = vmea(l)%k(m)
1258	j = vmea(l)%j(m)
1259	i = vmea(l)%i(m)
1260	vmea(l)%measured_vars(m,n) = 0.5_wp * &
1261	( u(k,j,i) + u(k,j,i+1) ) * &
1262	pt(k,j,i)
1263	ENDDO
1264
1265	CASE ( 'vtheta' )
1266	DO m = 1, vmea(l)%ns
1267	k = vmea(l)%k(m)
1268	j = vmea(l)%j(m)
1269	i = vmea(l)%i(m)
1270	vmea(l)%measured_vars(m,n) = 0.5_wp * &
1271	( v(k,j,i) + v(k,j+1,i) ) * &
1272	pt(k,j,i)
1273	ENDDO
1274
1275	CASE ( 'wtheta' )
1276	DO m = 1, vmea(l)%ns
1277	k = MAX ( 1, vmea(l)%k(m) )
1278	j = vmea(l)%j(m)
1279	i = vmea(l)%i(m)
1280	vmea(l)%measured_vars(m,n) = 0.5_wp * &
1281	( w(k-1,j,i) + w(k,j,i) ) * &
1282	pt(k,j,i)
1283	ENDDO
1284
1285	CASE ( 'uw' )
1286	DO m = 1, vmea(l)%ns
1287	k = MAX ( 1, vmea(l)%k(m) )
1288	j = vmea(l)%j(m)
1289	i = vmea(l)%i(m)
1290	vmea(l)%measured_vars(m,n) = 0.25_wp * &
1291	( w(k-1,j,i) + w(k,j,i) ) * &
1292	( u(k,j,i) + u(k,j,i+1) )
1293	ENDDO
1294
1295	CASE ( 'vw' )
1296	DO m = 1, vmea(l)%ns
1297	k = MAX ( 1, vmea(l)%k(m) )
1298	j = vmea(l)%j(m)
1299	i = vmea(l)%i(m)
1300	vmea(l)%measured_vars(m,n) = 0.25_wp * &
1301	( w(k-1,j,i) + w(k,j,i) ) * &
1302	( v(k,j,i) + v(k,j+1,i) )
1303	ENDDO
1304
1305	CASE ( 'uv' )
1306	DO m = 1, vmea(l)%ns
1307	k = MAX ( 1, vmea(l)%k(m) )
1308	j = vmea(l)%j(m)
1309	i = vmea(l)%i(m)
1310	vmea(l)%measured_vars(m,n) = 0.25_wp * &
1311	( u(k,j,i) + u(k,j,i+1) ) * &
1312	( v(k,j,i) + v(k,j+1,i) )
1313	ENDDO
1314	!
1315	!-- List of variables may need extension.
1316	CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mfco', 'mfno', 'mfno2', &
1317	'tro3' )
1318	IF ( air_chemistry ) THEN
1319	!
1320	!-- First, search for the measured variable in the chem_vars
1321	!-- list, in order to get the internal name of the variable.
1322	DO nn = 1, UBOUND( chem_vars, 2 )
1323	IF ( TRIM( vmea(l)%measured_vars_name(m) ) == &
1324	TRIM( chem_vars(0,nn) ) ) ind_chem = nn
1325	ENDDO
1326	!
1327	!-- Run loop over all chemical species, if the measured
1328	!-- variable matches the interal name, sample the variable.
1329	DO nn = 1, nvar
1330	IF ( TRIM( chem_vars(1,ind_chem) ) == &
1331	TRIM( chem_species(nn)%name ) ) THEN
1332	DO m = 1, vmea(l)%ns
1333	k = vmea(l)%k(m)
1334	j = vmea(l)%j(m)
1335	i = vmea(l)%i(m)
1336	vmea(l)%measured_vars(m,n) = &
1337	chem_species(nn)%conc(k,j,i)
1338	ENDDO
1339	ENDIF
1340	ENDDO
1341	ENDIF
1342
1343	CASE ( 'us' )
1344	DO m = 1, vmea(l)%ns
1345	!
1346	!-- Surface data is only available on inner subdomains, not
1347	!-- on ghost points. Hence, limit the indices.
1348	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1349	j = MERGE( j , nyn, j > nyn )
1350	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1351	i = MERGE( i , nxr, i > nxr )
1352
1353	DO mm = surf_def_h(0)%start_index(j,i), &
1354	surf_def_h(0)%end_index(j,i)
1355	vmea(l)%measured_vars(m,n) = surf_def_h(0)%us(mm)
1356	ENDDO
1357	DO mm = surf_lsm_h%start_index(j,i), &
1358	surf_lsm_h%end_index(j,i)
1359	vmea(l)%measured_vars(m,n) = surf_lsm_h%us(mm)
1360	ENDDO
1361	DO mm = surf_usm_h%start_index(j,i), &
1362	surf_usm_h%end_index(j,i)
1363	vmea(l)%measured_vars(m,n) = surf_usm_h%us(mm)
1364	ENDDO
1365	ENDDO
1366
1367	CASE ( 'ts' )
1368	DO m = 1, vmea(l)%ns
1369	!
1370	!-- Surface data is only available on inner subdomains, not
1371	!-- on ghost points. Hence, limit the indices.
1372	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1373	j = MERGE( j , nyn, j > nyn )
1374	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1375	i = MERGE( i , nxr, i > nxr )
1376
1377	DO mm = surf_def_h(0)%start_index(j,i), &
1378	surf_def_h(0)%end_index(j,i)
1379	vmea(l)%measured_vars(m,n) = surf_def_h(0)%ts(mm)
1380	ENDDO
1381	DO mm = surf_lsm_h%start_index(j,i), &
1382	surf_lsm_h%end_index(j,i)
1383	vmea(l)%measured_vars(m,n) = surf_lsm_h%ts(mm)
1384	ENDDO
1385	DO mm = surf_usm_h%start_index(j,i), &
1386	surf_usm_h%end_index(j,i)
1387	vmea(l)%measured_vars(m,n) = surf_usm_h%ts(mm)
1388	ENDDO
1389	ENDDO
1390
1391	CASE ( 'hfls' )
1392	DO m = 1, vmea(l)%ns
1393	!
1394	!-- Surface data is only available on inner subdomains, not
1395	!-- on ghost points. Hence, limit the indices.
1396	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1397	j = MERGE( j , nyn, j > nyn )
1398	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1399	i = MERGE( i , nxr, i > nxr )
1400
1401	DO mm = surf_def_h(0)%start_index(j,i), &
1402	surf_def_h(0)%end_index(j,i)
1403	vmea(l)%measured_vars(m,n) = surf_def_h(0)%qsws(mm)
1404	ENDDO
1405	DO mm = surf_lsm_h%start_index(j,i), &
1406	surf_lsm_h%end_index(j,i)
1407	vmea(l)%measured_vars(m,n) = surf_lsm_h%qsws(mm)
1408	ENDDO
1409	DO mm = surf_usm_h%start_index(j,i), &
1410	surf_usm_h%end_index(j,i)
1411	vmea(l)%measured_vars(m,n) = surf_usm_h%qsws(mm)
1412	ENDDO
1413	ENDDO
1414
1415	CASE ( 'hfss' )
1416	DO m = 1, vmea(l)%ns
1417	!
1418	!-- Surface data is only available on inner subdomains, not
1419	!-- on ghost points. Hence, limit the indices.
1420	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1421	j = MERGE( j , nyn, j > nyn )
1422	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1423	i = MERGE( i , nxr, i > nxr )
1424
1425	DO mm = surf_def_h(0)%start_index(j,i), &
1426	surf_def_h(0)%end_index(j,i)
1427	vmea(l)%measured_vars(m,n) = surf_def_h(0)%shf(mm)
1428	ENDDO
1429	DO mm = surf_lsm_h%start_index(j,i), &
1430	surf_lsm_h%end_index(j,i)
1431	vmea(l)%measured_vars(m,n) = surf_lsm_h%shf(mm)
1432	ENDDO
1433	DO mm = surf_usm_h%start_index(j,i), &
1434	surf_usm_h%end_index(j,i)
1435	vmea(l)%measured_vars(m,n) = surf_usm_h%shf(mm)
1436	ENDDO
1437	ENDDO
1438
1439	CASE ( 'rnds' )
1440	IF ( radiation ) THEN
1441	DO m = 1, vmea(l)%ns
1442	!
1443	!-- Surface data is only available on inner subdomains, not
1444	!-- on ghost points. Hence, limit the indices.
1445	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1446	j = MERGE( j , nyn, j > nyn )
1447	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1448	i = MERGE( i , nxr, i > nxr )
1449
1450	DO mm = surf_lsm_h%start_index(j,i), &
1451	surf_lsm_h%end_index(j,i)
1452	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_net(mm)
1453	ENDDO
1454	DO mm = surf_usm_h%start_index(j,i), &
1455	surf_usm_h%end_index(j,i)
1456	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_net(mm)
1457	ENDDO
1458	ENDDO
1459	ENDIF
1460
1461	CASE ( 'rsus' )
1462	IF ( radiation ) THEN
1463	DO m = 1, vmea(l)%ns
1464	!
1465	!-- Surface data is only available on inner subdomains, not
1466	!-- on ghost points. Hence, limit the indices.
1467	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1468	j = MERGE( j , nyn, j > nyn )
1469	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1470	i = MERGE( i , nxr, i > nxr )
1471
1472	DO mm = surf_lsm_h%start_index(j,i), &
1473	surf_lsm_h%end_index(j,i)
1474	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_out(mm)
1475	ENDDO
1476	DO mm = surf_usm_h%start_index(j,i), &
1477	surf_usm_h%end_index(j,i)
1478	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_out(mm)
1479	ENDDO
1480	ENDDO
1481	ENDIF
1482
1483	CASE ( 'rsds' )
1484	IF ( radiation ) THEN
1485	DO m = 1, vmea(l)%ns
1486	!
1487	!-- Surface data is only available on inner subdomains, not
1488	!-- on ghost points. Hence, limit the indices.
1489	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1490	j = MERGE( j , nyn, j > nyn )
1491	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1492	i = MERGE( i , nxr, i > nxr )
1493
1494	DO mm = surf_lsm_h%start_index(j,i), &
1495	surf_lsm_h%end_index(j,i)
1496	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_sw_in(mm)
1497	ENDDO
1498	DO mm = surf_usm_h%start_index(j,i), &
1499	surf_usm_h%end_index(j,i)
1500	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_sw_in(mm)
1501	ENDDO
1502	ENDDO
1503	ENDIF
1504
1505	CASE ( 'rlus' )
1506	IF ( radiation ) THEN
1507	DO m = 1, vmea(l)%ns
1508	!
1509	!-- Surface data is only available on inner subdomains, not
1510	!-- on ghost points. Hence, limit the indices.
1511	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1512	j = MERGE( j , nyn, j > nyn )
1513	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1514	i = MERGE( i , nxr, i > nxr )
1515
1516	DO mm = surf_lsm_h%start_index(j,i), &
1517	surf_lsm_h%end_index(j,i)
1518	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_out(mm)
1519	ENDDO
1520	DO mm = surf_usm_h%start_index(j,i), &
1521	surf_usm_h%end_index(j,i)
1522	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_out(mm)
1523	ENDDO
1524	ENDDO
1525	ENDIF
1526
1527	CASE ( 'rlds' )
1528	IF ( radiation ) THEN
1529	DO m = 1, vmea(l)%ns
1530	!
1531	!-- Surface data is only available on inner subdomains, not
1532	!-- on ghost points. Hence, limit the indices.
1533	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
1534	j = MERGE( j , nyn, j > nyn )
1535	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
1536	i = MERGE( i , nxr, i > nxr )
1537
1538	DO mm = surf_lsm_h%start_index(j,i), &
1539	surf_lsm_h%end_index(j,i)
1540	vmea(l)%measured_vars(m,n) = surf_lsm_h%rad_lw_in(mm)
1541	ENDDO
1542	DO mm = surf_usm_h%start_index(j,i), &
1543	surf_usm_h%end_index(j,i)
1544	vmea(l)%measured_vars(m,n) = surf_usm_h%rad_lw_in(mm)
1545	ENDDO
1546	ENDDO
1547	ENDIF
1548
1549	CASE ( 'rsd' )
1550	IF ( radiation ) THEN
1551	DO m = 1, vmea(l)%ns
1552	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
1553	j = vmea(l)%j(m)
1554	i = vmea(l)%i(m)
1555
1556	vmea(l)%measured_vars(m,n) = rad_sw_in(k,j,i)
1557	ENDDO
1558	ENDIF
1559
1560	CASE ( 'rsu' )
1561	IF ( radiation ) THEN
1562	DO m = 1, vmea(l)%ns
1563	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
1564	j = vmea(l)%j(m)
1565	i = vmea(l)%i(m)
1566
1567	vmea(l)%measured_vars(m,n) = rad_sw_out(k,j,i)
1568	ENDDO
1569	ENDIF
1570
1571	CASE ( 'rlu' )
1572	IF ( radiation ) THEN
1573	DO m = 1, vmea(l)%ns
1574	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
1575	j = vmea(l)%j(m)
1576	i = vmea(l)%i(m)
1577
1578	vmea(l)%measured_vars(m,n) = rad_lw_out(k,j,i)
1579	ENDDO
1580	ENDIF
1581
1582	CASE ( 'rld' )
1583	IF ( radiation ) THEN
1584	DO m = 1, vmea(l)%ns
1585	k = MERGE( 0, vmea(l)%k(m), radiation_scheme /= 'rrtmg' )
1586	j = vmea(l)%j(m)
1587	i = vmea(l)%i(m)
1588
1589	vmea(l)%measured_vars(m,n) = rad_lw_in(k,j,i)
1590	ENDDO
1591	ENDIF
1592
1593	CASE ( 'rsddif' )
1594	IF ( radiation ) THEN
1595	DO m = 1, vmea(l)%ns
1596	j = vmea(l)%j(m)
1597	i = vmea(l)%i(m)
1598
1599	vmea(l)%measured_vars(m,n) = rad_sw_in_diff(j,i)
1600	ENDDO
1601	ENDIF
1602
1603	CASE ( 't_soil' )
1604	DO m = 1, vmea(l)%ns_soil
1605	i = vmea(l)%i_soil(m)
1606	j = vmea(l)%j_soil(m)
1607	k = vmea(l)%k_soil(m)
1608
1609	match_lsm = surf_lsm_h%start_index(j,i) <= &
1610	surf_lsm_h%end_index(j,i)
1611	match_usm = surf_usm_h%start_index(j,i) <= &
1612	surf_usm_h%end_index(j,i)
1613
1614	IF ( match_lsm ) THEN
1615	mm = surf_lsm_h%start_index(j,i)
1616	vmea(l)%measured_vars_soil(m,n) = t_soil_h%var_2d(k,mm)
1617	ENDIF
1618
1619	IF ( match_usm ) THEN
1620	mm = surf_usm_h%start_index(j,i)
1621	vmea(l)%measured_vars_soil(m,n) = t_wall_h(k,mm)
1622	ENDIF
1623	ENDDO
1624
1625	CASE ( 'm_soil' )
1626	DO m = 1, vmea(l)%ns_soil
1627	i = vmea(l)%i_soil(m)
1628	j = vmea(l)%j_soil(m)
1629	k = vmea(l)%k_soil(m)
1630
1631	match_lsm = surf_lsm_h%start_index(j,i) <= &
1632	surf_lsm_h%end_index(j,i)
1633
1634	IF ( match_lsm ) THEN
1635	mm = surf_lsm_h%start_index(j,i)
1636	vmea(l)%measured_vars_soil(m,n) = m_soil_h%var_2d(k,mm)
1637	ENDIF
1638
1639	ENDDO
1640	!
1641	!-- More will follow ...
1642
1643	!
1644	!-- No match found - just set a fill value
1645	CASE DEFAULT
1646	vmea(l)%measured_vars(:,n) = vmea(l)%fillout
1647	END SELECT
1648
1649	ENDDO
1650
1651	ENDDO
1652
1653	END SUBROUTINE vm_sampling
1654
1655
1656	END MODULE virtual_measurement_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |