Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |