Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |