Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |