Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |