Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |