Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |