Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |