Home

Context Navigation

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

Last change on this file since 4332 was 4329, checked in by motisi, 5 years ago
Renamed wall_flags_0 to wall_flags_static_0
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 76.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |