Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |