Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |