Home

Context Navigation

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

Last change on this file since 4191 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 76.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |