Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |