Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |