Home

Context Navigation

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

Last change on this file since 3619 was 3522, checked in by suehring, 6 years ago
Sampling of variables extended in virtual measurement module
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 48.6 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	! -----------------
	22	!
[3473]	23	!
[3434]	24	! Former revisions:
	25	! -----------------
[3472]	26	! $Id: virtual_measurement_mod.f90 3522 2018-11-13 12:14:36Z suehring $
[3522]	27	! Sampling of variables
	28	!
	29	! 3494 2018-11-06 14:51:27Z suehring
[3494]	30	! Bugfixing
	31	!
	32	! 3473 2018-10-30 20:50:15Z suehring
[3473]	33	! Initial revision
[3434]	34	!
	35	! Authors:
	36	! --------
[3522]	37	! @author Matthias Suehring
[3434]	38	!
	39	! Description:
	40	! ------------
[3471]	41	!> The module acts as an interface between 'real-world' observations and
	42	!> model simulations. Virtual measurements will be taken in the model at the
	43	!> coordinates representative for the 'real-world' measurement positions.
	44	!> More precisely, coordinates and measured quanties will be read from a
	45	!> NetCDF file which contains all required information. In the model,
	46	!> the same quantities (as long as all the required components are switched-on)
	47	!> will be sampled at the respective positions and output into an extra file,
	48	!> which allows for straight-forward comparison of model results with
	49	!> observations.
[3522]	50	!>
	51	!> @todo list_of_allowed variables needs careful checking
	52	!> @todo output (binary or NetCDF) needs to be implemented
	53	!> @todo clean-up anything from current test modus
	54	!> @todo Check if sign of surface fluxes for heat, radiation, etc., follows
	55	!> the (UC)2 standard
	56	!> @note Fluxes are not processed
[3434]	57	!------------------------------------------------------------------------------!
[3471]	58	MODULE virtual_measurement_mod
[3434]	59
	60	USE arrays_3d, &
	61	ONLY: q, pt, u, v, w, zu, zw
	62
[3522]	63	USE chem_modules, &
	64	ONLY: nspec
	65
	66	USE chemistry_model_mod, &
	67	ONLY: chem_species
	68
[3434]	69	USE control_parameters, &
[3522]	70	ONLY: air_chemistry, dz, humidity, neutral, message_string, &
	71	virtual_measurement
[3434]	72
	73	USE cpulog, &
	74	ONLY: cpu_log, log_point
	75
	76	USE grid_variables, &
	77	ONLY: dx, dy
	78
	79	USE indices, &
[3522]	80	ONLY: nzb, nzt, nxl, nxr, nys, nyn, nx, ny
[3434]	81
	82	USE kinds
	83
	84
	85	IMPLICIT NONE
	86
	87	TYPE virt_mea
	88
	89	CHARACTER(LEN=100) :: feature_type !< type of the measurement
	90	CHARACTER(LEN=100) :: site !< name of the measurement site
	91
	92	CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables
	93
[3522]	94	INTEGER(iwp) :: ns = 0 !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
	95	INTEGER(iwp) :: ntraj !< number of trajectories of a measurement
	96	INTEGER(iwp) :: nvar !< number of measured variables
[3434]	97
	98	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues
	99
	100	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
	101	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
	102	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
	103
	104	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
	105	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
	106	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
	107
	108	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
	109	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
	110	REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values
	111	REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building
	112	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
	113	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
[3522]	114
[3434]	115	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
	116
	117	END TYPE virt_mea
	118
	119	CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting
	120	CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type
	121	CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue
	122	CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names
	123	CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing
	124	CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations
	125	CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x
	126	CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y
	127	CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name
	128	CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable
	129	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
	130	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
	131	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
[3522]	132	!
	133	!-- MS: List requires careful revision!
	134	CHARACTER(LEN=10), DIMENSION(1:47), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM
[3471]	135	(/ 'hfls ', & ! surface latent heat flux (W/m2)
	136	'hfss ', & ! surface sensible heat flux (W/m2)
	137	'hur ', & ! relative humidity (-)
	138	'hus ', & ! specific humidity (g/kg)
	139	'haa ', & ! absolute atmospheric humidity (kg/m3)
	140	'mcpm1 ', & ! mass concentration of PM1 (kg/m3)
	141	'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3)
	142	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
	143	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
	144	'mcco ', & ! mass concentration of CO (kg/m3)
	145	'mcco2 ', & ! mass concentration of CO2 (kg/m3)
	146	'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3)
	147	'ncaa ', & ! number concentation of particles (1/m3)
[3522]	148	'mfco ', & ! mole fraction of CO (mol/mol)
[3471]	149	'mfco2 ', & ! mole fraction of CO2 (mol/mol)
	150	'mfch4 ', & ! mole fraction of methane (mol/mol)
	151	'mfnh3 ', & ! mole fraction of amonia (mol/mol)
	152	'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol)
	153	'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol)
	154	'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol)
	155	'mfh20 ', & ! mole fraction of water (mol/mol)
	156	'plev ', & ! ? air pressure - hydrostaic + perturbation?
	157	'rlds ', & ! surface downward longwave flux (W/m2)
	158	'rlus ', & ! surface upward longwave flux (W/m2)
	159	'rsds ', & ! surface downward shortwave flux (W/m2)
	160	'rsus ', & ! surface upward shortwave flux (W/m2)
	161	'ta ', & ! air temperature (degree C)
	162	't_va ', & ! virtual accoustic temperature (K)
	163	'theta ', & ! potential temperature (K)
	164	'tro3 ', & ! mole fraction of ozone air (mol/mol)
	165	'ts ', & ! scaling parameter of temperature (K)
	166	'wspeed ', & ! ? wind speed - horizontal?
	167	'wdir ', & ! wind direction
	168	'us ', & ! friction velocity
	169	'msoil ', & ! ? soil moisture - which depth?
	170	'tsoil ', & ! ? soil temperature - which depth?
	171	'u ', & ! u-component
	172	'ua ', & ! eastward wind (is there any difference to u?)
	173	'v ', & ! v-component
	174	'va ', & ! northward wind (is there any difference to v?)
	175	'w ', & ! w-component
	176	'rld ', & ! downward longwave radiative flux (W/m2)
	177	'rlu ', & ! upnward longwave radiative flux (W/m2)
	178	'rsd ', & ! downward shortwave radiative flux (W/m2)
	179	'rsu ', & ! upward shortwave radiative flux (W/m2)
	180	'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2)
	181	'rnds ' & ! surface net downward radiative flux (W/m2)
	182	/)
[3434]	183
	184	INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements
	185	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
	186	INTEGER(iwp) :: observation_coverage_xy = 0 !< horizontal distance from the measurement point where observations should be taken in the surrounding
	187	INTEGER(iwp) :: observation_coverage_z = 0 !< vertical distance from the measurement point where observations should be taken in the surrounding
	188
	189	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
	190	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
	191
	192	REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start
	193
	194
	195	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure
	196
	197
	198	INTERFACE vm_check_parameters
	199	MODULE PROCEDURE vm_check_parameters
	200	END INTERFACE vm_check_parameters
	201
	202	INTERFACE vm_init
	203	MODULE PROCEDURE vm_init
	204	END INTERFACE vm_init
	205
	206	INTERFACE vm_parin
	207	MODULE PROCEDURE vm_parin
	208	END INTERFACE vm_parin
	209
	210	INTERFACE vm_sampling
	211	MODULE PROCEDURE vm_sampling
	212	END INTERFACE vm_sampling
	213
	214	SAVE
	215
	216	PRIVATE
	217
	218	!
	219	!-- Public interfaces
	220	PUBLIC vm_check_parameters, vm_init, vm_parin, vm_sampling
	221
	222	!
	223	!-- Public variables
	224	PUBLIC vmea, vm_time_start
	225
	226	CONTAINS
	227
	228
	229	!------------------------------------------------------------------------------!
	230	! Description:
	231	! ------------
[3471]	232	!> Check parameters for virtual measurement module
[3434]	233	!------------------------------------------------------------------------------!
	234	SUBROUTINE vm_check_parameters
	235
	236	USE control_parameters, &
	237	ONLY: message_string, virtual_measurement
	238
	239	USE netcdf_data_input_mod, &
	240	ONLY: input_pids_static
	241
	242	IMPLICIT NONE
	243
	244	!
	245	!-- In case virtual measurements are taken, a static input file is required.
	246	!-- This is because UTM coordinates for the PALM domain origin are required
	247	!-- for correct mapping of the measurements.
	248	!-- ToDo: Revise this later and remove this requirement.
	249	IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN
	250	message_string = 'If virtual measurements are taken a static input ' // &
	251	'file is mandatory.'
	252	CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
	253	ENDIF
	254
	255	END SUBROUTINE vm_check_parameters
	256
	257	!------------------------------------------------------------------------------!
	258	! Description:
	259	! ------------
[3471]	260	!> Read namelist for the virtual measurement module
[3434]	261	!------------------------------------------------------------------------------!
	262	SUBROUTINE vm_parin
	263
	264	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	265
	266	NAMELIST /virtual_measurement_parameters/ use_virtual_measurement, &
	267	vm_time_start
	268
	269	line = ' '
	270
	271	!
	272	!-- Try to find stg package
	273	REWIND ( 11 )
	274	line = ' '
	275	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
	276	READ ( 11, '(A)', END=20 ) line
	277	ENDDO
	278	BACKSPACE ( 11 )
	279
	280	!
	281	!-- Read namelist
	282	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
	283
	284	!
[3471]	285	!-- Set flag that indicates that the virtual measurement module is switched on
[3434]	286	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
	287
	288	GOTO 20
	289
	290	10 BACKSPACE( 11 )
	291	READ( 11 , '(A)') line
	292	CALL parin_fail_message( 'virtual_measurement_parameters', line )
	293
	294	20 CONTINUE
	295
	296	END SUBROUTINE vm_parin
	297
	298
	299	!------------------------------------------------------------------------------!
	300	! Description:
	301	! ------------
	302	!> Initialize virtual measurements: read coordiante arrays and measured
	303	!> variables, set indicies indicating the measurement points, read further
	304	!> attributes, etc..
	305	!------------------------------------------------------------------------------!
	306	SUBROUTINE vm_init
	307
	308	USE arrays_3d, &
	309	ONLY: zu, zw
	310
	311	USE grid_variables, &
	312	ONLY: ddx, ddy, dx, dy
	313
	314	USE indices, &
	315	ONLY: nxl, nxr, nyn, nys
	316
	317	USE netcdf_data_input_mod, &
	318	ONLY: init_model, input_file_vm, &
	319	netcdf_data_input_get_dimension_length, &
	320	netcdf_data_input_att, netcdf_data_input_var
	321
	322	USE surface_mod, &
	323	ONLY: get_topography_top_index_ji
	324
	325	IMPLICIT NONE
	326
	327	CHARACTER(LEN=5) :: dum !< dummy string indicate station id
	328	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
[3522]	329	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
[3434]	330
	331	INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate
	332	INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate
	333	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
	334	INTEGER(iwp) :: dim_zag !< dimension size of height coordinate
	335	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
	336	INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement
	337	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
	338	INTEGER(iwp) :: j !< grid index of observation point in x-direction
	339	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
	340	INTEGER(iwp) :: k !< grid index of observation point in x-direction
[3522]	341	INTEGER(iwp) :: kl !< lower vertical index of surrounding grid points of an observation coordinate
[3434]	342	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
	343	INTEGER(iwp) :: ksurf !< topography top index
[3522]	344	INTEGER(iwp) :: ku !< upper vertical index of surrounding grid points of an observation coordinate
[3434]	345	INTEGER(iwp) :: l !< running index over all stations
	346	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
	347	INTEGER(iwp) :: ll !< running index over all measured variables in file
	348	INTEGER(iwp) :: lll !< running index over all allowed variables
	349	INTEGER(iwp) :: n !< running index over trajectory coordinates
	350	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
	351	INTEGER(iwp) :: t !< running index over number of trajectories
	352
[3522]	353	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: meas_flag !< mask array indicating measurement positions
	354
	355	LOGICAL :: chem_include !< flag indicating that chemical species is considered in modelled mechanism
	356	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
	357
[3434]	358	REAL(wp) :: fill_eutm !< _FillValue for coordinate array E_UTM
	359	REAL(wp) :: fill_nutm !< _FillValue for coordinate array N_UTM
	360	REAL(wp) :: fill_zag !< _FillValue for height coordinate
	361
[3437]	362	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: e_utm !< easting UTM coordinate, temporary variable
	363	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: n_utm !< northing UTM coordinate, temporary variable,
	364	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z_ag !< height coordinate relative to origin_z, temporary variable
[3434]	365	!
	366	!-- Obtain number of virtual measurement stations
	367	CALL netcdf_data_input_att( nvm, char_numstations, id_vm, input_file_vm, &
	368	global_attribute, 'open', '' )
[3522]	369
	370	! write(9,*) "num stationi", nvm
	371	! flush(9)
[3434]	372	!
	373	!-- ALLOCATE data structure
	374	ALLOCATE( vmea(1:nvm) )
	375	!
[3522]	376	!-- Allocate flag array
	377	ALLOCATE( meas_flag(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	378	meas_flag = 0
	379
	380	!
[3434]	381	!-- Read station coordinates and further attributes.
	382	!-- Note all coordinates are in UTM coordinates.
	383	DO l = 1, nvm
	384	!
	385	!-- Determine suffix which contains the ID
	386	IF( l < 10 ) THEN
	387	WRITE( dum, '(I1)') l
	388	ELSEIF( l < 100 ) THEN
	389	WRITE( dum, '(I2)') l
	390	ELSEIF( l < 1000 ) THEN
	391	WRITE( dum, '(I3)') l
	392	ELSEIF( l < 10000 ) THEN
	393	WRITE( dum, '(I4)') l
	394	ELSEIF( l < 100000 ) THEN
	395	WRITE( dum, '(I5)') l
	396	ENDIF
	397
	398	CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx &
	399	// TRIM( dum ), id_vm, '', global_attribute,&
	400	'', '' )
	401	CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy &
	402	// TRIM( dum ), id_vm, '', global_attribute,&
	403	'', '' )
	404	CALL netcdf_data_input_att( vmea(l)%site, char_site &
	405	// TRIM( dum ), id_vm, '', global_attribute,&
	406	'', '' )
	407	CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature &
	408	// TRIM( dum ), id_vm, '', global_attribute,&
	409	'', '' )
	410
	411	!
	412	!--- Set logicals depending on the type of the measurement
	413	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
	414	vmea(l)%timseries_profile = .TRUE.
	415	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
	416	vmea(l)%timseries = .TRUE.
	417	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
	418	vmea(l)%trajectory = .TRUE.
	419	ELSE
	420	!
	421	!-- Give error message
	422	message_string = 'Attribue featureType = ' // &
	423	TRIM( vmea(l)%feature_type ) // &
	424	' is not allowed.'
	425	CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
	426	ENDIF
	427	!
	428	!-- Read string with all measured variables at this station
	429	measured_variables_file = ''
	430	CALL netcdf_data_input_var( measured_variables_file, &
	431	char_mv // TRIM( dum ), id_vm )
	432	!
	433	!-- Count the number of measured variables which match with the variables
	434	!-- which are allowed to be measured in PALM. Please note, for some
	435	!-- NetCDF interal reasons characters end with a NULL, i.e. also empty
	436	!-- characters contain a NULL. Therefore, check the strings for a Null to
	437	!-- get the correct character length in order to compare them with the list
	438	!-- of allowed variables.
	439	vmea(l)%nvar = 0
	440	DO ll = 1, SIZE( measured_variables_file )
	441	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
	442	measured_variables_file(ll)(1:1) /= ' ') THEN
	443	!
	444	!-- Obtain character length of the character
	445	len_char = 1
	446	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
	447	.AND. measured_variables_file(ll)(len_char:len_char) /= ' ' )
	448	len_char = len_char + 1
	449	ENDDO
	450	len_char = len_char - 1
	451	!
	452	!-- Now, compare the measured variable with the list of allowed
	453	!-- variables.
	454	DO lll= 1, SIZE( list_allowed_variables )
	455	IF ( measured_variables_file(ll)(1:len_char) == &
	456	TRIM( list_allowed_variables(lll) ) ) THEN
	457	vmea(l)%nvar = vmea(l)%nvar + 1
	458	measured_variables(vmea(l)%nvar) = &
	459	measured_variables_file(ll)(1:len_char)
	460	ENDIF
	461	ENDDO
	462	ENDIF
	463	ENDDO
	464	!
	465	!-- Allocate array for the measured variables names for the station l.
	466	ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
	467
	468	DO ll = 1, vmea(l)%nvar
	469	vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
	470	ENDDO
	471	!
[3522]	472	!-- In case of chemistry, check if species is considered in the modelled
	473	!-- chemistry mechanism.
	474	IF ( air_chemistry ) THEN
	475	DO ll = 1, vmea(l)%nvar
	476	chem_include = .FALSE.
	477	DO n = 1, nspec
	478	IF ( TRIM( vmea(l)%measured_vars_name(ll) ) == &
	479	TRIM( chem_species(n)%name ) ) chem_include = .TRUE.
	480	ENDDO
	481	IF ( .NOT. chem_include ) THEN
	482	message_string = TRIM( vmea(l)%measured_vars_name(ll) ) // &
	483	' is not considered in the modelled ' // &
	484	'chemistry mechanism'
	485	CALL message( 'vm_init', 'PA0000', 0, 0, 0, 6, 0 )
	486	ENDIF
	487	ENDDO
	488	ENDIF
	489	!
[3434]	490	!-- For the actual measurement ID read the UTM coordinates. Based on these,
	491	!-- define the index space on each subdomain where measurements should be
	492	!-- taken. Note, the entire coordinate arrays will not be stored on data
	493	!-- type as this would exceed memory requirements, particularly for
	494	!-- trajectory measurements. If no variable will be virtually measured,
	495	!-- skip the reading.
	496	IF ( vmea(l)%nvar > 0 ) THEN
	497	!
	498	!-- For stationary measurements UTM coordinates are just one value and
	499	!-- its dimension is "station", while for mobile measurements UTM
	500	!-- coordinates are arrays. First, inquire dimension length for
	501	!-- UTM coordinates.
	502	IF ( vmea(l)%trajectory ) THEN
	503	!
	504	!-- For non-stationary measurements read the number of trajectories
	505	CALL netcdf_data_input_get_dimension_length( id_vm, &
	506	vmea(l)%ntraj, &
	507	"traj" // &
	508	TRIM( dum ) )
	509	CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime, &
	510	"ntime" // &
	511	TRIM( dum ) )
	512	!
	513	!-- For stationary measurements the dimension for UTM coordinates is 1
	514	ELSE
	515	vmea(l)%ntraj = 1
	516	dim_ntime = 1
	517	ENDIF
	518
	519	!
	520	!- Allocate array which defines individual time frame for each
	521	!-- trajectory or station
	522	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
	523	!
	524	!-- Allocate temporary arrays for UTM and height coordinates. Note,
	525	!-- on file UTM coordinates might be 1D or 2D variables
[3437]	526	ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	527	ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	528	ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) )
[3434]	529	!
	530	!-- Read _FillValue attributes
	531	CALL netcdf_data_input_att( fill_eutm, char_fillvalue, &
	532	id_vm, '', .NOT. global_attribute, '', &
	533	char_eutm // TRIM( dum ) )
	534	CALL netcdf_data_input_att( fill_nutm, char_fillvalue, &
	535	id_vm, '', .NOT. global_attribute, '', &
	536	char_nutm // TRIM( dum ) )
	537	CALL netcdf_data_input_att( fill_zag, char_fillvalue, &
	538	id_vm, '', .NOT. global_attribute, '', &
	539	char_zag // TRIM( dum ) )
	540	!
	541	!-- Read UTM and height coordinates coordinates for all trajectories and
	542	!-- times.
[3437]	543	IF ( vmea(l)%trajectory ) THEN
	544	CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm, &
	545	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	546	CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm, &
	547	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	548	CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm, &
	549	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	550	ELSE
[3471]	551	CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
[3437]	552	CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
	553	CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), id_vm )
	554	ENDIF
[3434]	555	!
[3522]	556	!-- For testing:
	557	e_utm = e_utm - e_utm(1,1)
	558	n_utm = n_utm - n_utm(1,1)
	559
	560	!
	561	!-- First, compute relative x- and y-coordinates with respect to the
	562	!-- lower-left origin of the model domain, which is the difference
	563	!-- betwen UTM coordinates.
	564	! e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t)) &
	565	! - init_model%origin_x
	566	! n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t)) &
	567	! - init_model%origin_y
	568
	569	!
[3434]	570	!-- Based on UTM coordinates, check if the measurement station or parts
	571	!-- of the trajectory is on subdomain. This case, setup grid index space
	572	!-- sample these quantities.
[3522]	573	meas_flag = 0
[3434]	574	DO t = 1, vmea(l)%ntraj
	575	!
	576	!-- Determine the individual time coordinate length for each station and
	577	!-- trajectory. This is required as several stations and trajectories
	578	!-- are merged into one file but they do not have the same number of
	579	!-- points in time, hence, missing values may occur and cannot be
	580	!-- processed further.
	581	vmea(l)%dim_t(t) = 0
	582	DO n = 1, dim_ntime
[3437]	583	IF ( e_utm(t,n) /= fill_eutm .AND. &
	584	n_utm(t,n) /= fill_nutm .AND. &
	585	z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n
[3434]	586	ENDDO
[3522]	587
[3434]	588	!
	589	!-- Compute grid indices relative to origin and check if these are
	590	!-- on the subdomain. Note, virtual measurements will be taken also
	591	!-- at grid points surrounding the station, hence, check also for
	592	!-- these grid points.
[3437]	593	DO n = 1, vmea(l)%dim_t(t)
	594	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
	595	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
[3434]	596	!
	597	!-- Is the observation point on subdomain?
	598	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
	599	js >= nys .AND. js <= nyn )
	600	!
[3522]	601	!-- Check if observation coordinate is on subdomain
[3434]	602	IF ( on_pe ) THEN
[3522]	603	!
	604	!-- Determine vertical index which correspond to the observation
	605	!-- height.
[3434]	606	ksurf = get_topography_top_index_ji( js, is, 's' )
[3437]	607	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
[3434]	608	!
[3522]	609	!-- Set mask array at the observation coordinates. Also, flag the
	610	!-- surrounding coordinate points, but first check whether the
	611	!-- surrounding coordinate points are on the subdomain.
	612	kl = MERGE( ks-1, ks, ks-1 >= nzb .AND. ks-1 > ksurf )
	613	ku = MERGE( ks+1, ks, ks+1 <= nzt+1 )
	614
	615	meas_flag(kl:ku,js-1:js+1,is-1:is+1) = &
	616	IBSET( meas_flag(kl:ku,js-1:js+1,is-1:is+1), 0 )
[3434]	617	ENDIF
	618	ENDDO
	619
	620	ENDDO
	621	!
[3522]	622	!-- Based on the flag array count the number of observation points.
	623	ns = 0
	624	DO is = nxl-1, nxr+1
	625	DO js = nys-1, nyn+1
	626	DO ks = nzb, nzt+1
	627	ns = ns + MERGE( 1, 0, BTEST( meas_flag(ks,js,is), 0 ) )
	628	ENDDO
	629	ENDDO
	630	ENDDO
	631	!
[3434]	632	!-- Store number of observation points on subdomain and allocate index
	633	!-- arrays.
	634	vmea(l)%ns = ns
[3522]	635	ns = 0
[3434]	636
	637	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
	638	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
	639	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
	640	!
[3522]	641	!-- Based on the flag array store the grid indices which correspond to
	642	!-- the observation coordinates.
	643	DO is = nxl-1, nxr+1
	644	DO js = nys-1, nyn+1
	645	DO ks = nzb, nzt+1
	646	IF ( BTEST( meas_flag(ks,js,is), 0 ) ) THEN
	647	ns = ns + 1
	648	vmea(l)%i(ns) = is
	649	vmea(l)%j(ns) = js
	650	vmea(l)%k(ns) = ks
	651	ENDIF
	652	ENDDO
[3434]	653	ENDDO
	654	ENDDO
[3522]	655
	656	! write(9,*) l, "NS: ", ns
	657	! IF ( ns > 0 ) THEN
	658	! write(9,*) "i ", vmea(l)%i
	659	! write(9,*) "j ", vmea(l)%j
	660	! write(9,*) "k ", vmea(l)%k
	661	! write(9,*)
	662	! ENDIF
[3434]	663	!
	664	!-- Allocate array to save the sampled values.
	665	!-- Todo: Is it better to allocate for all variables at a station
	666	!-- and store all the values before writing, or sample the variables
	667	!-- directly in the data output?
	668	ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
	669	!
	670	!-- Initialize with _FillValue
	671	vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
	672	!
	673	!-- Deallocate temporary coordinate arrays
	674	IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm )
	675	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
	676	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag )
	677	ENDIF
	678	ENDDO
[3522]	679	! flush(9)
[3434]	680
	681	!
	682	!-- Close input file for virtual measurements. Therefore, just call
	683	!-- the read attribute routine with the "close" option.
	684	CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '', &
	685	global_attribute, 'close', '' )
[3522]	686	!
	687	!-- Dellocate flag array
	688	DEALLOCATE( meas_flag )
	689
[3434]	690	END SUBROUTINE vm_init
	691
	692
	693	!------------------------------------------------------------------------------!
	694	! Description:
	695	! ------------
	696	!> Sampling of the actual quantities along the observation coordinates
	697	!------------------------------------------------------------------------------!
[3471]	698	SUBROUTINE vm_sampling
[3434]	699
[3522]	700	USE arrays_3d, &
	701	ONLY: exner, pt, q, u, v, w
[3471]	702
[3522]	703	USE basic_constants_and_equations_mod, &
	704	ONLY: pi
[3434]	705
[3522]	706	USE radiation_model_mod, &
	707	ONLY: radiation
	708
	709	USE surface_mod, &
	710	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	711
[3434]	712	IMPLICIT NONE
	713
	714	CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name
	715
[3522]	716	INTEGER(iwp) :: i !< grid index in x-direction
	717	INTEGER(iwp) :: j !< grid index in y-direction
	718	INTEGER(iwp) :: k !< grid index in z-direction
	719	INTEGER(iwp) :: l !< running index over the number of stations
	720	INTEGER(iwp) :: m !< running index over all virtual observation coordinates
	721	INTEGER(iwp) :: mm !< index of surface element which corresponds to the virtual observation coordinate
	722	INTEGER(iwp) :: n !< running index over all measured variables at a station
	723	INTEGER(iwp) :: nn !< running index over the number of chemcal species
[3434]	724	!
	725	!-- Loop over all stations. For each possible variable loop over all
	726	!-- observation points
	727	DO l = 1, nvm
	728	!
[3522]	729	!-- Loop over all measured variables at this station. Please note,
	730	!-- velocity components are interpolated onto scalar grid.
	731	DO n = 1, vmea(l)%nvar
	732
	733	SELECT CASE ( TRIM( vmea(l)%measured_vars_name(n) ) )
	734
	735	CASE ( 'theta' )
	736	IF ( .NOT. neutral ) THEN
	737	DO m = 1, vmea(l)%ns
	738	k = vmea(l)%k(m)
	739	j = vmea(l)%j(m)
	740	i = vmea(l)%i(m)
	741	vmea(l)%measured_vars(n,m) = pt(k,j,i)
	742	ENDDO
	743	ENDIF
	744
	745	CASE ( 'ta', 't_va' )
	746	IF ( .NOT. neutral ) THEN
	747	DO m = 1, vmea(l)%ns
	748	k = vmea(l)%k(m)
	749	j = vmea(l)%j(m)
	750	i = vmea(l)%i(m)
	751	vmea(l)%measured_vars(n,m) = pt(k,j,i) * exner( k )
	752	ENDDO
	753	ENDIF
	754
	755	CASE ( 'hus', 'haa' )
	756	IF ( humidity ) THEN
	757	DO m = 1, vmea(l)%ns
	758	k = vmea(l)%k(m)
	759	j = vmea(l)%j(m)
	760	i = vmea(l)%i(m)
	761	vmea(l)%measured_vars(n,m) = q(k,j,i)
	762	ENDDO
	763	ENDIF
	764
	765	CASE ( 'u', 'ua' )
	766	DO m = 1, vmea(l)%ns
	767	k = vmea(l)%k(m)
	768	j = vmea(l)%j(m)
	769	i = vmea(l)%i(m)
	770	vmea(l)%measured_vars(n,m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	771	ENDDO
	772
	773	CASE ( 'v', 'va' )
	774	DO m = 1, vmea(l)%ns
	775	k = vmea(l)%k(m)
	776	j = vmea(l)%j(m)
	777	i = vmea(l)%i(m)
	778	vmea(l)%measured_vars(n,m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	779	ENDDO
	780
	781	CASE ( 'w' )
	782	DO m = 1, vmea(l)%ns
	783	k = vmea(l)%k(m)
	784	j = vmea(l)%j(m)
	785	i = vmea(l)%i(m)
	786	vmea(l)%measured_vars(n,m) = w(k,j,i) !0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
	787	ENDDO
	788
	789	CASE ( 'wspeed' )
	790	DO m = 1, vmea(l)%ns
	791	k = vmea(l)%k(m)
	792	j = vmea(l)%j(m)
	793	i = vmea(l)%i(m)
	794	vmea(l)%measured_vars(n,m) = SQRT( &
	795	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 + &
	796	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 &
	797	)
	798	ENDDO
	799
	800	CASE ( 'wdir' )
	801	DO m = 1, vmea(l)%ns
	802	k = vmea(l)%k(m)
	803	j = vmea(l)%j(m)
	804	i = vmea(l)%i(m)
	805
	806	vmea(l)%measured_vars(n,m) = ATAN2( &
	807	- 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ), &
	808	- 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) &
	809	) * 180.0_wp / pi
	810	ENDDO
	811	!
	812	!-- MS: list of variables needs extension.
	813	CASE ( 'mcpm1', 'mcpm2p5', 'mcpm10', 'mcco', 'mcco2', 'mcbcda', &
	814	'ncaa', 'mfco2', 'mfco', 'mfch4', 'mfnh3', 'mfno', &
	815	'mfno2', 'mfso2', 'mfh20', 'tr03' )
	816	IF ( air_chemistry ) THEN
	817	DO nn = 1, nspec
	818	IF ( TRIM( vmea(l)%measured_vars_name(m) ) == &
	819	TRIM( chem_species(nn)%name ) ) THEN
	820	DO m = 1, vmea(l)%ns
	821	k = vmea(l)%k(m)
	822	j = vmea(l)%j(m)
	823	i = vmea(l)%i(m)
	824	vmea(l)%measured_vars(n,m) = &
	825	chem_species(nn)%conc(k,j,i)
	826	ENDDO
	827	ENDIF
	828	ENDDO
	829	ENDIF
	830
	831	CASE ( 'us' )
	832	DO m = 1, vmea(l)%ns
	833	!
	834	!-- Surface data is only available on inner subdomains, not
	835	!-- on ghost points. Hence, limit the indices.
	836	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	837	j = MERGE( j , nyn, j > nyn )
	838	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	839	i = MERGE( i , nxr, i > nxr )
	840
	841	DO mm = surf_def_h(0)%start_index(j,i), &
	842	surf_def_h(0)%end_index(j,i)
	843	vmea(l)%measured_vars(n,m) = surf_def_h(0)%us(mm)
	844	ENDDO
	845	DO mm = surf_lsm_h%start_index(j,i), &
	846	surf_lsm_h%end_index(j,i)
	847	vmea(l)%measured_vars(n,m) = surf_lsm_h%us(mm)
	848	ENDDO
	849	DO mm = surf_usm_h%start_index(j,i), &
	850	surf_usm_h%end_index(j,i)
	851	vmea(l)%measured_vars(n,m) = surf_usm_h%us(mm)
	852	ENDDO
	853	ENDDO
	854
	855	CASE ( 'ts' )
	856	DO m = 1, vmea(l)%ns
	857	!
	858	!-- Surface data is only available on inner subdomains, not
	859	!-- on ghost points. Hence, limit the indices.
	860	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	861	j = MERGE( j , nyn, j > nyn )
	862	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	863	i = MERGE( i , nxr, i > nxr )
	864
	865	DO mm = surf_def_h(0)%start_index(j,i), &
	866	surf_def_h(0)%end_index(j,i)
	867	vmea(l)%measured_vars(n,m) = surf_def_h(0)%ts(mm)
	868	ENDDO
	869	DO mm = surf_lsm_h%start_index(j,i), &
	870	surf_lsm_h%end_index(j,i)
	871	vmea(l)%measured_vars(n,m) = surf_lsm_h%ts(mm)
	872	ENDDO
	873	DO mm = surf_usm_h%start_index(j,i), &
	874	surf_usm_h%end_index(j,i)
	875	vmea(l)%measured_vars(n,m) = surf_usm_h%ts(mm)
	876	ENDDO
	877	ENDDO
	878
	879	CASE ( 'hfls' )
	880	DO m = 1, vmea(l)%ns
	881	!
	882	!-- Surface data is only available on inner subdomains, not
	883	!-- on ghost points. Hence, limit the indices.
	884	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	885	j = MERGE( j , nyn, j > nyn )
	886	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	887	i = MERGE( i , nxr, i > nxr )
	888
	889	DO mm = surf_def_h(0)%start_index(j,i), &
	890	surf_def_h(0)%end_index(j,i)
	891	vmea(l)%measured_vars(n,m) = surf_def_h(0)%qsws(mm)
	892	ENDDO
	893	DO mm = surf_lsm_h%start_index(j,i), &
	894	surf_lsm_h%end_index(j,i)
	895	vmea(l)%measured_vars(n,m) = surf_lsm_h%qsws(mm)
	896	ENDDO
	897	DO mm = surf_usm_h%start_index(j,i), &
	898	surf_usm_h%end_index(j,i)
	899	vmea(l)%measured_vars(n,m) = surf_usm_h%qsws(mm)
	900	ENDDO
	901	ENDDO
	902
	903	CASE ( 'hfss' )
	904	DO m = 1, vmea(l)%ns
	905	!
	906	!-- Surface data is only available on inner subdomains, not
	907	!-- on ghost points. Hence, limit the indices.
	908	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	909	j = MERGE( j , nyn, j > nyn )
	910	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	911	i = MERGE( i , nxr, i > nxr )
	912
	913	DO mm = surf_def_h(0)%start_index(j,i), &
	914	surf_def_h(0)%end_index(j,i)
	915	vmea(l)%measured_vars(n,m) = surf_def_h(0)%shf(mm)
	916	ENDDO
	917	DO mm = surf_lsm_h%start_index(j,i), &
	918	surf_lsm_h%end_index(j,i)
	919	vmea(l)%measured_vars(n,m) = surf_lsm_h%shf(mm)
	920	ENDDO
	921	DO mm = surf_usm_h%start_index(j,i), &
	922	surf_usm_h%end_index(j,i)
	923	vmea(l)%measured_vars(n,m) = surf_usm_h%shf(mm)
	924	ENDDO
	925	ENDDO
	926
	927	CASE ( 'rnds' )
	928	IF ( radiation ) THEN
	929	DO m = 1, vmea(l)%ns
	930	!
	931	!-- Surface data is only available on inner subdomains, not
	932	!-- on ghost points. Hence, limit the indices.
	933	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	934	j = MERGE( j , nyn, j > nyn )
	935	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	936	i = MERGE( i , nxr, i > nxr )
	937
	938	DO mm = surf_lsm_h%start_index(j,i), &
	939	surf_lsm_h%end_index(j,i)
	940	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_net(mm)
	941	ENDDO
	942	DO mm = surf_usm_h%start_index(j,i), &
	943	surf_usm_h%end_index(j,i)
	944	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_net(mm)
	945	ENDDO
	946	ENDDO
	947	ENDIF
	948
	949	CASE ( 'rsus', 'rsu' )
	950	IF ( radiation ) THEN
	951	DO m = 1, vmea(l)%ns
	952	!
	953	!-- Surface data is only available on inner subdomains, not
	954	!-- on ghost points. Hence, limit the indices.
	955	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	956	j = MERGE( j , nyn, j > nyn )
	957	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	958	i = MERGE( i , nxr, i > nxr )
	959
	960	DO mm = surf_lsm_h%start_index(j,i), &
	961	surf_lsm_h%end_index(j,i)
	962	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_out(mm)
	963	ENDDO
	964	DO mm = surf_usm_h%start_index(j,i), &
	965	surf_usm_h%end_index(j,i)
	966	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_out(mm)
	967	ENDDO
	968	ENDDO
	969	ENDIF
	970
	971	CASE ( 'rsds', 'rsd' )
	972	IF ( radiation ) THEN
	973	DO m = 1, vmea(l)%ns
	974	!
	975	!-- Surface data is only available on inner subdomains, not
	976	!-- on ghost points. Hence, limit the indices.
	977	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	978	j = MERGE( j , nyn, j > nyn )
	979	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	980	i = MERGE( i , nxr, i > nxr )
	981
	982	DO mm = surf_lsm_h%start_index(j,i), &
	983	surf_lsm_h%end_index(j,i)
	984	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_sw_in(mm)
	985	ENDDO
	986	DO mm = surf_usm_h%start_index(j,i), &
	987	surf_usm_h%end_index(j,i)
	988	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_sw_in(mm)
	989	ENDDO
	990	ENDDO
	991	ENDIF
	992
	993	CASE ( 'rlus', 'rlu' )
	994	IF ( radiation ) THEN
	995	DO m = 1, vmea(l)%ns
	996	!
	997	!-- Surface data is only available on inner subdomains, not
	998	!-- on ghost points. Hence, limit the indices.
	999	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1000	j = MERGE( j , nyn, j > nyn )
	1001	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1002	i = MERGE( i , nxr, i > nxr )
	1003
	1004	DO mm = surf_lsm_h%start_index(j,i), &
	1005	surf_lsm_h%end_index(j,i)
	1006	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_out(mm)
	1007	ENDDO
	1008	DO mm = surf_usm_h%start_index(j,i), &
	1009	surf_usm_h%end_index(j,i)
	1010	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_out(mm)
	1011	ENDDO
	1012	ENDDO
	1013	ENDIF
	1014
	1015	CASE ( 'rlds', 'rld' )
	1016	IF ( radiation ) THEN
	1017	DO m = 1, vmea(l)%ns
	1018	!
	1019	!-- Surface data is only available on inner subdomains, not
	1020	!-- on ghost points. Hence, limit the indices.
	1021	j = MERGE( vmea(l)%j(m), nys, vmea(l)%j(m) < nys )
	1022	j = MERGE( j , nyn, j > nyn )
	1023	i = MERGE( vmea(l)%i(m), nxl, vmea(l)%i(m) < nxl )
	1024	i = MERGE( i , nxr, i > nxr )
	1025
	1026	DO mm = surf_lsm_h%start_index(j,i), &
	1027	surf_lsm_h%end_index(j,i)
	1028	vmea(l)%measured_vars(n,m) = surf_lsm_h%rad_lw_in(mm)
	1029	ENDDO
	1030	DO mm = surf_usm_h%start_index(j,i), &
	1031	surf_usm_h%end_index(j,i)
	1032	vmea(l)%measured_vars(n,m) = surf_usm_h%rad_lw_in(mm)
	1033	ENDDO
	1034	ENDDO
	1035	ENDIF
	1036	!
	1037	!-- More will follow ...
	1038
	1039	END SELECT
	1040
[3494]	1041	ENDDO
[3434]	1042
	1043	ENDDO
	1044
[3471]	1045	END SUBROUTINE vm_sampling
[3434]	1046
	1047
[3471]	1048	END MODULE virtual_measurement_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |