Home

Context Navigation

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

Last change on this file since 3890 was 3876, checked in by knoop, 6 years ago
Moved "photolysis_scheme", "chem_species" and "phot_frequen" to chem_modules
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 73.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |