Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |