Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 76.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |