Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |