Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |