Home

Context Navigation

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

Last change on this file since 3912 was 3910, checked in by suehring, 6 years ago
Bugfix in rotation of UTM coordinates
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 74.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |