Home

Context Navigation

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

Last change on this file since 3479 was 3473, checked in by suehring, 6 years ago
Bugfix for previous commit
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 35.9 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 3473 2018-10-30 20:50:15Z gronemeier $
[3473]	27	! Initial revision
	28	!
	29	! 3472 2018-10-30 20:43:50Z suehring
[3434]	30	!
	31	! Authors:
	32	! --------
	33	! @author Matthias Suehring
	34	!
	35	!
	36	!
	37	! Description:
	38	! ------------
[3471]	39	!> The module acts as an interface between 'real-world' observations and
	40	!> model simulations. Virtual measurements will be taken in the model at the
	41	!> coordinates representative for the 'real-world' measurement positions.
	42	!> More precisely, coordinates and measured quanties will be read from a
	43	!> NetCDF file which contains all required information. In the model,
	44	!> the same quantities (as long as all the required components are switched-on)
	45	!> will be sampled at the respective positions and output into an extra file,
	46	!> which allows for straight-forward comparison of model results with
	47	!> observations.
[3434]	48	!------------------------------------------------------------------------------!
[3471]	49	MODULE virtual_measurement_mod
[3434]	50
	51
	52	USE arrays_3d, &
	53	ONLY: q, pt, u, v, w, zu, zw
	54
	55	USE control_parameters, &
	56	ONLY: dz, message_string, virtual_measurement
	57
	58	USE cpulog, &
	59	ONLY: cpu_log, log_point
	60
	61	USE grid_variables, &
	62	ONLY: dx, dy
	63
	64	USE indices, &
	65	ONLY: nzb, nzt, nxl, nxr, nys, nyn
	66
	67	USE kinds
	68
	69
	70	IMPLICIT NONE
	71
	72	TYPE virt_mea
	73
	74	CHARACTER(LEN=100) :: feature_type !< type of the measurement
	75	CHARACTER(LEN=100) :: site !< name of the measurement site
	76
	77	CHARACTER(LEN=10), DIMENSION(:), ALLOCATABLE :: measured_vars_name !< name of the measured variables
	78
	79	INTEGER(iwp) :: ns !< total number of observation points for a site on subdomain, i.e. sum of all trajectories
	80	INTEGER(iwp) :: ntraj !< number of trajectories of a measurement
	81	INTEGER(iwp) :: nvar !< number of measured variables
	82
	83	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: dim_t !< number observations individual for each trajectory or station that are no _FillValues
	84	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ngp !< number of grid points where observations for a site took place,
	85	!<individual for each trajectory or station that are no _FillValues
	86
	87	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< grid index for measurement position in x-direction
	88	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< grid index for measurement position in y-direction
	89	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< grid index for measurement position in k-direction
	90
	91	LOGICAL :: trajectory = .FALSE. !< flag indicating that the observation is a mobile observation
	92	LOGICAL :: timseries = .FALSE. !< flag indicating that the observation is a stationary point measurement
	93	LOGICAL :: timseries_profile = .FALSE. !< flag indicating that the observation is a stationary profile measurement
	94
	95	REAL(wp) :: fill_eutm !< fill value for UTM coordinates in case of missing values
	96	REAL(wp) :: fill_nutm !< fill value for UTM coordinates in case of missing values
	97	REAL(wp) :: fill_zag !< fill value for heigth coordinates in case of missing values
	98	REAL(wp) :: fillout = -999.9 !< fill value for output in case a observation is taken from inside a building
	99	REAL(wp) :: origin_x_obs !< origin of the observation in UTM coordiates in x-direction
	100	REAL(wp) :: origin_y_obs !< origin of the observation in UTM coordiates in y-direction
	101
	102	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: xmea !< measurement x-position in absolute UTM coordinates
	103	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ymea !< measurement y-position in absolute UTM coordinates
	104	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmea !< measurement z-position in height above ground level
	105
	106	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: measured_vars !< measured variables
	107
	108	END TYPE virt_mea
	109
	110	CHARACTER(LEN=5) :: char_eutm = "E_UTM" !< dimension name for UTM coordinate easting
	111	CHARACTER(LEN=11) :: char_feature = "featureType" !< attribute name for feature type
	112	CHARACTER(LEN=10) :: char_fillvalue = "_FillValue" !< variable attribute name for _FillValue
	113	CHARACTER(LEN=18) :: char_mv = "measured_variables" !< variable name for the array with the measured variable names
	114	CHARACTER(LEN=5) :: char_nutm = "N_UTM" !< dimension name for UTM coordinate northing
	115	CHARACTER(LEN=18) :: char_numstations = "number_of_stations" !< attribute name for number of stations
	116	CHARACTER(LEN=8) :: char_origx = "origin_x" !< attribute name for station coordinate in x
	117	CHARACTER(LEN=8) :: char_origy = "origin_y" !< attribute name for station coordinate in y
	118	CHARACTER(LEN=4) :: char_site = "site" !< attribute name for site name
	119	CHARACTER(LEN=19) :: char_zag = "height_above_ground" !< attribute name for height above ground variable
	120	CHARACTER(LEN=10) :: type_ts = 'timeSeries' !< name of stationary point measurements
	121	CHARACTER(LEN=10) :: type_traj = 'trajectory' !< name of line measurements
	122	CHARACTER(LEN=17) :: type_tspr = 'timeSeriesProfile' !< name of stationary profile measurements
	123
[3471]	124	CHARACTER(LEN=10), DIMENSION(1:53), PARAMETER :: list_allowed_variables = & !< variables that can be sampled in PALM
	125	(/ 'hfls ', & ! surface latent heat flux (W/m2)
	126	'hfss ', & ! surface sensible heat flux (W/m2)
	127	'hur ', & ! relative humidity (-)
	128	'hus ', & ! specific humidity (g/kg)
	129	'haa ', & ! absolute atmospheric humidity (kg/m3)
	130	'mcpm1 ', & ! mass concentration of PM1 (kg/m3)
	131	'mcpm2p5 ', & ! mass concentration of PM2.5 (kg/m3)
	132	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
	133	'mcpm10 ', & ! mass concentration of PM10 (kg/m3)
	134	'mcco ', & ! mass concentration of CO (kg/m3)
	135	'mcco2 ', & ! mass concentration of CO2 (kg/m3)
	136	'mcbcda ', & ! mass concentration of black carbon paritcles (kg/m3)
	137	'ncaa ', & ! number concentation of particles (1/m3)
	138	'mfco2 ', & ! mole fraction of CO (mol/mol)
	139	'mfco2 ', & ! mole fraction of CO2 (mol/mol)
	140	'mfch4 ', & ! mole fraction of methane (mol/mol)
	141	'mfnh3 ', & ! mole fraction of amonia (mol/mol)
	142	'mfno ', & ! mole fraction of nitrogen monoxide (mol/mol)
	143	'mfno2 ', & ! mole fraction of nitrogen dioxide (mol/mol)
	144	'mfso2 ', & ! mole fraction of sulfur dioxide (mol/mol)
	145	'mfh20 ', & ! mole fraction of water (mol/mol)
	146	'plev ', & ! ? air pressure - hydrostaic + perturbation?
	147	'rlds ', & ! surface downward longwave flux (W/m2)
	148	'rlus ', & ! surface upward longwave flux (W/m2)
	149	'rsds ', & ! surface downward shortwave flux (W/m2)
	150	'rsus ', & ! surface upward shortwave flux (W/m2)
	151	'ta ', & ! air temperature (degree C)
	152	't_va ', & ! virtual accoustic temperature (K)
	153	'theta ', & ! potential temperature (K)
	154	'tro3 ', & ! mole fraction of ozone air (mol/mol)
	155	'ts ', & ! scaling parameter of temperature (K)
	156	'wspeed ', & ! ? wind speed - horizontal?
	157	'wdir ', & ! wind direction
	158	'us ', & ! friction velocity
	159	'msoil ', & ! ? soil moisture - which depth?
	160	'tsoil ', & ! ? soil temperature - which depth?
	161	'u ', & ! u-component
	162	'ua ', & ! eastward wind (is there any difference to u?)
	163	'uw ', & ! ? vertical momentum flux - total ?
	164	'utheta ', & ! ? horizontal heat flux - total ?
	165	'uv ', & ! upward-northward horizontal momentum flux
	166	'v ', & ! v-component
	167	'va ', & ! northward wind (is there any difference to v?)
	168	'vw ', & ! ? vertical momentum flux - total ?
	169	'vtheta ', & ! ? horizontal heat flux - total ?
	170	'w ', & ! w-component
	171	'wtheta ', & ! ? vertical heat flux - total ?
	172	'rld ', & ! downward longwave radiative flux (W/m2)
	173	'rlu ', & ! upnward longwave radiative flux (W/m2)
	174	'rsd ', & ! downward shortwave radiative flux (W/m2)
	175	'rsu ', & ! upward shortwave radiative flux (W/m2)
	176	'rsddif ', & ! downward shortwave diffuse radiative flux (W/m2)
	177	'rnds ' & ! surface net downward radiative flux (W/m2)
	178	/)
[3434]	179
	180	INTEGER(iwp) :: id_vm !< NetCDF file id for virtual measurements
	181	INTEGER(iwp) :: nvm = 0 !< number of virtual measurements
	182	INTEGER(iwp) :: observation_coverage_xy = 0 !< horizontal distance from the measurement point where observations should be taken in the surrounding
	183	INTEGER(iwp) :: observation_coverage_z = 0 !< vertical distance from the measurement point where observations should be taken in the surrounding
	184
	185	LOGICAL :: use_virtual_measurement = .FALSE. !< Namelist parameter
	186	LOGICAL :: global_attribute = .TRUE. !< flag indicating a global attribute
	187
	188	REAL(wp) :: vm_time_start = 0.0 !< time after virtual measurements should start
	189
	190
	191	TYPE( virt_mea ), DIMENSION(:), ALLOCATABLE :: vmea !< virtual measurement data structure
	192
	193
	194	INTERFACE vm_check_parameters
	195	MODULE PROCEDURE vm_check_parameters
	196	END INTERFACE vm_check_parameters
	197
	198	INTERFACE vm_init
	199	MODULE PROCEDURE vm_init
	200	END INTERFACE vm_init
	201
	202	INTERFACE vm_parin
	203	MODULE PROCEDURE vm_parin
	204	END INTERFACE vm_parin
	205
	206	INTERFACE vm_sampling
	207	MODULE PROCEDURE vm_sampling
	208	END INTERFACE vm_sampling
	209
	210	SAVE
	211
	212	PRIVATE
	213
	214	!
	215	!-- Public interfaces
	216	PUBLIC vm_check_parameters, vm_init, vm_parin, vm_sampling
	217
	218	!
	219	!-- Public variables
	220	PUBLIC vmea, vm_time_start
	221
	222	CONTAINS
	223
	224
	225	!------------------------------------------------------------------------------!
	226	! Description:
	227	! ------------
[3471]	228	!> Check parameters for virtual measurement module
[3434]	229	!------------------------------------------------------------------------------!
	230	SUBROUTINE vm_check_parameters
	231
	232	USE control_parameters, &
	233	ONLY: message_string, virtual_measurement
	234
	235	USE netcdf_data_input_mod, &
	236	ONLY: input_pids_static
	237
	238	IMPLICIT NONE
	239
	240	!
	241	!-- In case virtual measurements are taken, a static input file is required.
	242	!-- This is because UTM coordinates for the PALM domain origin are required
	243	!-- for correct mapping of the measurements.
	244	!-- ToDo: Revise this later and remove this requirement.
	245	IF ( virtual_measurement .AND. .NOT. input_pids_static ) THEN
	246	message_string = 'If virtual measurements are taken a static input ' // &
	247	'file is mandatory.'
	248	CALL message( 'vm_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
	249	ENDIF
	250
	251	END SUBROUTINE vm_check_parameters
	252
	253	!------------------------------------------------------------------------------!
	254	! Description:
	255	! ------------
[3471]	256	!> Read namelist for the virtual measurement module
[3434]	257	!------------------------------------------------------------------------------!
	258	SUBROUTINE vm_parin
	259
	260	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	261
	262	NAMELIST /virtual_measurement_parameters/ use_virtual_measurement, &
	263	vm_time_start
	264
	265	line = ' '
	266
	267	!
	268	!-- Try to find stg package
	269	REWIND ( 11 )
	270	line = ' '
	271	DO WHILE ( INDEX( line, '&virtual_measurement_parameters' ) == 0 )
	272	READ ( 11, '(A)', END=20 ) line
	273	ENDDO
	274	BACKSPACE ( 11 )
	275
	276	!
	277	!-- Read namelist
	278	READ ( 11, virtual_measurement_parameters, ERR = 10, END = 20 )
	279
	280	!
[3471]	281	!-- Set flag that indicates that the virtual measurement module is switched on
[3434]	282	IF ( use_virtual_measurement ) virtual_measurement = .TRUE.
	283
	284	GOTO 20
	285
	286	10 BACKSPACE( 11 )
	287	READ( 11 , '(A)') line
	288	CALL parin_fail_message( 'virtual_measurement_parameters', line )
	289
	290	20 CONTINUE
	291
	292	END SUBROUTINE vm_parin
	293
	294
	295	!------------------------------------------------------------------------------!
	296	! Description:
	297	! ------------
	298	!> Initialize virtual measurements: read coordiante arrays and measured
	299	!> variables, set indicies indicating the measurement points, read further
	300	!> attributes, etc..
	301	!------------------------------------------------------------------------------!
	302	SUBROUTINE vm_init
	303
	304	USE arrays_3d, &
	305	ONLY: zu, zw
	306
	307	USE control_parameters, &
	308	ONLY: message_string
	309
	310	USE grid_variables, &
	311	ONLY: ddx, ddy, dx, dy
	312
	313	USE indices, &
	314	ONLY: nxl, nxr, nyn, nys
	315
	316	USE netcdf_data_input_mod, &
	317	ONLY: init_model, input_file_vm, &
	318	netcdf_data_input_get_dimension_length, &
	319	netcdf_data_input_att, netcdf_data_input_var
	320
	321	USE surface_mod, &
	322	ONLY: get_topography_top_index_ji
	323
	324	IMPLICIT NONE
	325
	326	CHARACTER(LEN=5) :: dum !< dummy string indicate station id
	327	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables_file = '' !< array with all measured variables read from NetCDF
	328	CHARACTER(LEN=10), DIMENSION(50) :: measured_variables = '' !< dummy array with all measured variables that are allowed
	329
	330	LOGICAL :: on_pe !< flag indicating that the respective measurement coordinate is on subdomain
	331
	332	INTEGER(iwp) :: dim_eutm !< dimension size of UTM easting coordinate
	333	INTEGER(iwp) :: dim_nutm !< dimension size of UTM northing coordinate
	334	INTEGER(iwp) :: dim_ntime !< dimension size of time coordinate
	335	INTEGER(iwp) :: dim_zag !< dimension size of height coordinate
	336	INTEGER(iwp) :: i !< grid index of virtual observation point in x-direction
	337	INTEGER(iwp) :: icov !< index range where observations should be taken in x-direction
	338	INTEGER(iwp) :: ii !< running index over all coordinate points of a measurement
[3471]	339	INTEGER(iwp) :: i_prev !< grid index along x for UTM coordinate at previous observation time step
[3434]	340	INTEGER(iwp) :: is !< grid index of real observation point of the respective station in x-direction
	341	INTEGER(iwp) :: j !< grid index of observation point in x-direction
	342	INTEGER(iwp) :: jcov !< index range where observations should be taken in y-direction
[3471]	343	INTEGER(iwp) :: j_prev !< grid index along y for UTM coordinate at previous observation time step
[3434]	344	INTEGER(iwp) :: js !< grid index of real observation point of the respective station in y-direction
	345	INTEGER(iwp) :: k !< grid index of observation point in x-direction
	346	INTEGER(iwp) :: kcov !< index range where observations should be taken in z-direction
	347	INTEGER(iwp) :: ks !< grid index of real observation point of the respective station in z-direction
[3471]	348	INTEGER(iwp) :: k_prev !< grid index along z for UTM coordinate at previous observation time step
[3434]	349	INTEGER(iwp) :: ksurf !< topography top index
	350	INTEGER(iwp) :: l !< running index over all stations
	351	INTEGER(iwp) :: len_char !< character length of single measured variables without Null character
	352	INTEGER(iwp) :: ll !< running index over all measured variables in file
	353	INTEGER(iwp) :: lll !< running index over all allowed variables
	354	INTEGER(iwp) :: n !< running index over trajectory coordinates
	355	INTEGER(iwp) :: ns !< counter variable for number of observation points on subdomain
	356	INTEGER(iwp) :: t !< running index over number of trajectories
	357
	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', '' )
	369	!
	370	!-- ALLOCATE data structure
	371	ALLOCATE( vmea(1:nvm) )
	372
[3437]	373	! print*, "nvm", nvm
[3434]	374	!
	375	!-- Read station coordinates and further attributes.
	376	!-- Note all coordinates are in UTM coordinates.
	377	DO l = 1, nvm
	378	!
	379	!-- Determine suffix which contains the ID
	380	IF( l < 10 ) THEN
	381	WRITE( dum, '(I1)') l
	382	ELSEIF( l < 100 ) THEN
	383	WRITE( dum, '(I2)') l
	384	ELSEIF( l < 1000 ) THEN
	385	WRITE( dum, '(I3)') l
	386	ELSEIF( l < 10000 ) THEN
	387	WRITE( dum, '(I4)') l
	388	ELSEIF( l < 100000 ) THEN
	389	WRITE( dum, '(I5)') l
	390	ENDIF
	391
	392	CALL netcdf_data_input_att( vmea(l)%origin_x_obs, char_origx &
	393	// TRIM( dum ), id_vm, '', global_attribute,&
	394	'', '' )
	395	CALL netcdf_data_input_att( vmea(l)%origin_y_obs, char_origy &
	396	// TRIM( dum ), id_vm, '', global_attribute,&
	397	'', '' )
	398	CALL netcdf_data_input_att( vmea(l)%site, char_site &
	399	// TRIM( dum ), id_vm, '', global_attribute,&
	400	'', '' )
	401	CALL netcdf_data_input_att( vmea(l)%feature_type, char_feature &
	402	// TRIM( dum ), id_vm, '', global_attribute,&
	403	'', '' )
	404
	405	!
	406	!--- Set logicals depending on the type of the measurement
	407	IF ( INDEX( vmea(l)%feature_type, type_tspr ) /= 0 ) THEN
	408	vmea(l)%timseries_profile = .TRUE.
	409	ELSEIF ( INDEX( vmea(l)%feature_type, type_ts ) /= 0 ) THEN
	410	vmea(l)%timseries = .TRUE.
	411	ELSEIF ( INDEX( vmea(l)%feature_type, type_traj ) /= 0 ) THEN
	412	vmea(l)%trajectory = .TRUE.
	413	ELSE
	414	!
	415	!-- Give error message
	416	message_string = 'Attribue featureType = ' // &
	417	TRIM( vmea(l)%feature_type ) // &
	418	' is not allowed.'
	419	CALL message( 'vm_init', 'PA0000', 1, 2, 0, 6, 0 )
	420	ENDIF
	421	!
	422	!-- Read string with all measured variables at this station
	423	measured_variables_file = ''
	424	CALL netcdf_data_input_var( measured_variables_file, &
	425	char_mv // TRIM( dum ), id_vm )
	426	!
	427	!-- Count the number of measured variables which match with the variables
	428	!-- which are allowed to be measured in PALM. Please note, for some
	429	!-- NetCDF interal reasons characters end with a NULL, i.e. also empty
	430	!-- characters contain a NULL. Therefore, check the strings for a Null to
	431	!-- get the correct character length in order to compare them with the list
	432	!-- of allowed variables.
	433	vmea(l)%nvar = 0
	434	DO ll = 1, SIZE( measured_variables_file )
	435	IF ( measured_variables_file(ll)(1:1) /= CHAR(0) .AND. &
	436	measured_variables_file(ll)(1:1) /= ' ') THEN
	437	!
	438	!-- Obtain character length of the character
	439	len_char = 1
	440	DO WHILE ( measured_variables_file(ll)(len_char:len_char) /= CHAR(0)&
	441	.AND. measured_variables_file(ll)(len_char:len_char) /= ' ' )
	442	len_char = len_char + 1
	443	ENDDO
	444	len_char = len_char - 1
	445	!
	446	!-- Now, compare the measured variable with the list of allowed
	447	!-- variables.
	448	DO lll= 1, SIZE( list_allowed_variables )
	449	IF ( measured_variables_file(ll)(1:len_char) == &
	450	TRIM( list_allowed_variables(lll) ) ) THEN
	451	vmea(l)%nvar = vmea(l)%nvar + 1
	452	measured_variables(vmea(l)%nvar) = &
	453	measured_variables_file(ll)(1:len_char)
	454	ENDIF
	455	ENDDO
	456	ENDIF
	457	ENDDO
	458	!
	459	!-- Allocate array for the measured variables names for the station l.
	460	ALLOCATE( vmea(l)%measured_vars_name(1:vmea(l)%nvar) )
	461
	462	DO ll = 1, vmea(l)%nvar
	463	vmea(l)%measured_vars_name(ll) = TRIM( measured_variables(ll) )
	464	ENDDO
	465
[3437]	466	! print*, "numvars", vmea(l)%nvar, vmea(l)%measured_vars_name(1:vmea(l)%nvar)
[3434]	467	!
	468	!-- For the actual measurement ID read the UTM coordinates. Based on these,
	469	!-- define the index space on each subdomain where measurements should be
	470	!-- taken. Note, the entire coordinate arrays will not be stored on data
	471	!-- type as this would exceed memory requirements, particularly for
	472	!-- trajectory measurements. If no variable will be virtually measured,
	473	!-- skip the reading.
	474	IF ( vmea(l)%nvar > 0 ) THEN
	475	!
	476	!-- For stationary measurements UTM coordinates are just one value and
	477	!-- its dimension is "station", while for mobile measurements UTM
	478	!-- coordinates are arrays. First, inquire dimension length for
	479	!-- UTM coordinates.
	480	IF ( vmea(l)%trajectory ) THEN
	481	!
	482	!-- For non-stationary measurements read the number of trajectories
	483	CALL netcdf_data_input_get_dimension_length( id_vm, &
	484	vmea(l)%ntraj, &
	485	"traj" // &
	486	TRIM( dum ) )
	487	CALL netcdf_data_input_get_dimension_length( id_vm, dim_ntime, &
	488	"ntime" // &
	489	TRIM( dum ) )
	490	!
	491	!-- For stationary measurements the dimension for UTM coordinates is 1
	492	ELSE
	493	vmea(l)%ntraj = 1
	494	dim_ntime = 1
	495	ENDIF
	496
	497	!
	498	!- Allocate array which defines individual time frame for each
	499	!-- trajectory or station
	500	ALLOCATE( vmea(l)%dim_t(1:vmea(l)%ntraj) )
[3437]	501	ALLOCATE( vmea(l)%ngp(1:vmea(l)%ntraj) )
[3434]	502	!
	503	!-- Allocate temporary arrays for UTM and height coordinates. Note,
	504	!-- on file UTM coordinates might be 1D or 2D variables
[3437]	505	ALLOCATE( e_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	506	ALLOCATE( n_utm(1:vmea(l)%ntraj,1:dim_ntime) )
	507	ALLOCATE( z_ag(1:vmea(l)%ntraj,1:dim_ntime) )
[3434]	508	!
	509	!-- Read _FillValue attributes
	510	CALL netcdf_data_input_att( fill_eutm, char_fillvalue, &
	511	id_vm, '', .NOT. global_attribute, '', &
	512	char_eutm // TRIM( dum ) )
	513	CALL netcdf_data_input_att( fill_nutm, char_fillvalue, &
	514	id_vm, '', .NOT. global_attribute, '', &
	515	char_nutm // TRIM( dum ) )
	516	CALL netcdf_data_input_att( fill_zag, char_fillvalue, &
	517	id_vm, '', .NOT. global_attribute, '', &
	518	char_zag // TRIM( dum ) )
	519	!
	520	!-- Read UTM and height coordinates coordinates for all trajectories and
	521	!-- times.
[3437]	522	IF ( vmea(l)%trajectory ) THEN
	523	CALL netcdf_data_input_var( e_utm, char_eutm // TRIM( dum ), id_vm, &
	524	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	525	CALL netcdf_data_input_var( n_utm, char_nutm // TRIM( dum ), id_vm, &
	526	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	527	CALL netcdf_data_input_var( z_ag, char_zag // TRIM( dum ), id_vm, &
	528	0, dim_ntime-1, 0, vmea(l)%ntraj-1 )
	529	ELSE
[3471]	530	CALL netcdf_data_input_var( e_utm(1,:), char_eutm // TRIM( dum ), id_vm )
[3437]	531	CALL netcdf_data_input_var( n_utm(1,:), char_nutm // TRIM( dum ), id_vm )
	532	CALL netcdf_data_input_var( z_ag(1,:), char_zag // TRIM( dum ), id_vm )
	533	ENDIF
[3434]	534	!
	535	!-- Based on UTM coordinates, check if the measurement station or parts
	536	!-- of the trajectory is on subdomain. This case, setup grid index space
	537	!-- sample these quantities.
	538	ns = 0
	539	DO t = 1, vmea(l)%ntraj
	540	!
	541	!-- Determine the individual time coordinate length for each station and
	542	!-- trajectory. This is required as several stations and trajectories
	543	!-- are merged into one file but they do not have the same number of
	544	!-- points in time, hence, missing values may occur and cannot be
	545	!-- processed further.
	546	vmea(l)%dim_t(t) = 0
	547	DO n = 1, dim_ntime
[3437]	548	IF ( e_utm(t,n) /= fill_eutm .AND. &
	549	n_utm(t,n) /= fill_nutm .AND. &
	550	z_ag(t,n) /= fill_zag ) vmea(l)%dim_t(t) = n
[3434]	551	ENDDO
	552	!
	553	!-- First, compute relative x- and y-coordinates with respect to the
	554	!-- lower-left origin of the model domain, which is the difference
	555	!-- betwen UTM coordinates.
[3471]	556	! write(9,*) l, "before Eutm", e_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_x
	557	! write(9,*) l, "before Nutm", n_utm(t,1:vmea(l)%dim_t(t)), "orig", init_model%origin_y
	558
[3437]	559	e_utm(t,1:vmea(l)%dim_t(t)) = e_utm(t,1:vmea(l)%dim_t(t)) &
	560	- init_model%origin_x
	561	n_utm(t,1:vmea(l)%dim_t(t)) = n_utm(t,1:vmea(l)%dim_t(t)) &
	562	- init_model%origin_y
[3471]	563	! write(9,*) l, "Eutm", e_utm(t,1:vmea(l)%dim_t(t))
	564	! write(9,*) l, "Nutm", n_utm(t,1:vmea(l)%dim_t(t))
	565	! write(9,*)
[3434]	566	!
	567	!-- Compute grid indices relative to origin and check if these are
	568	!-- on the subdomain. Note, virtual measurements will be taken also
	569	!-- at grid points surrounding the station, hence, check also for
	570	!-- these grid points.
	571	vmea(l)%ngp(t) = 0
	572	k_prev = -999
	573	j_prev = -999
	574	i_prev = -999
[3437]	575	DO n = 1, vmea(l)%dim_t(t)
	576	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
	577	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
[3434]	578	!
	579	!-- Is the observation point on subdomain?
	580	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
	581	js >= nys .AND. js <= nyn )
	582	!
	583	!-- If the measurement is on subdomain, determine vertical index
	584	!-- which refers to the observation height above ground level.
	585	ks = k_prev
	586	IF ( on_pe ) THEN
[3471]	587	! write(9,*) " I am onpe"
[3434]	588	ksurf = get_topography_top_index_ji( js, is, 's' )
[3437]	589	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
[3434]	590	ENDIF
	591	!
	592	!-- Count the number of observation points in index space on
	593	!-- subdomain. Only increment if grid indices are different from
	594	!-- the previous one.
	595	IF ( on_pe .AND. is /= i_prev .AND. js /= j_prev .AND. &
	596	ks /= k_prev ) THEN
	597	ns = ns + 1
	598	vmea(l)%ngp(t) = vmea(l)%ngp(t) + 1
	599	ENDIF
[3437]	600
[3434]	601	!-- Store arrays for next iteration - avoid double counting
	602	i_prev = is
	603	j_prev = js
	604	k_prev = ks
	605	ENDDO
	606
	607	ENDDO
	608
	609	!
	610	!-- Store number of observation points on subdomain and allocate index
	611	!-- arrays.
	612	vmea(l)%ns = ns
	613
	614	ALLOCATE( vmea(l)%i(1:vmea(l)%ns) )
	615	ALLOCATE( vmea(l)%j(1:vmea(l)%ns) )
	616	ALLOCATE( vmea(l)%k(1:vmea(l)%ns) )
	617
[3437]	618	! print*, "Num ns: ", vmea(l)%ns, "per traj", vmea(l)%ngp(:)
[3434]	619	!
	620	!-- Repeat the prior loop and save the grid indices relevant for
	621	!-- sampling.
	622	ns = 0
	623	DO t = 1, vmea(l)%ntraj
	624	!
	625	!-- Compute grid indices relative to origin and check if these are
	626	!-- on the subdomain. Note, virtual measurements will be taken also
	627	!-- at grid points surrounding the station, hence, check also for
	628	!-- these grid points.
	629	k_prev = -999
	630	j_prev = -999
	631	i_prev = -999
[3437]	632	DO n = 1, vmea(l)%dim_t(t)
	633	is = INT( ( e_utm(t,n) + 0.5_wp * dx ) * ddx, KIND = iwp )
	634	js = INT( ( n_utm(t,n) + 0.5_wp * dy ) * ddy, KIND = iwp )
[3434]	635	!
	636	!-- Is the observation point on subdomain?
	637	on_pe = ( is >= nxl .AND. is <= nxr .AND. &
	638	js >= nys .AND. js <= nyn )
	639	!
	640	!-- If the measurement is on subdomain, determine vertical index
	641	!-- which refers to the observation height above ground level.
	642	ks = k_prev
	643	IF ( on_pe ) THEN
	644	ksurf = get_topography_top_index_ji( js, is, 's' )
[3437]	645	ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
[3434]	646	ENDIF
	647	!
	648	!-- Count the number of observation points in index space on
	649	!-- subdomain. Only increment if grid indices are different from
	650	!-- the previous one.
	651	IF ( on_pe .AND. is /= i_prev .AND. js /= j_prev .AND. &
	652	ks /= k_prev ) THEN
	653	ns = ns + 1
	654	vmea(l)%i(ns) = is
	655	vmea(l)%j(ns) = js
	656	vmea(l)%k(ns) = ks
[3471]	657	! write(9,*) TRIM( vmea(l)%feature_type), l, "ns", ns, "ijk", vmea(l)%i(ns), vmea(l)%j(ns), vmea(l)%k(ns)
[3434]	658	ENDIF
	659	!
	660	!-- Store arrays for next iteration - avoid double counting
	661	i_prev = is
	662	j_prev = js
	663	k_prev = ks
	664	ENDDO
	665
	666	ENDDO
	667	!
	668	!-- Allocate array to save the sampled values.
	669	!-- Todo: Is it better to allocate for all variables at a station
	670	!-- and store all the values before writing, or sample the variables
	671	!-- directly in the data output?
	672	ALLOCATE( vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) )
	673	!
	674	!-- Initialize with _FillValue
	675	vmea(l)%measured_vars(1:vmea(l)%nvar,1:vmea(l)%ns) = vmea(l)%fillout
	676	!
	677	!-- Deallocate temporary coordinate arrays
	678	IF ( ALLOCATED( e_utm ) ) DEALLOCATE( e_utm )
	679	IF ( ALLOCATED( n_utm ) ) DEALLOCATE( n_utm )
	680	IF ( ALLOCATED( z_ag ) ) DEALLOCATE( z_ag )
	681	ENDIF
	682	ENDDO
[3471]	683	flush(9)
[3434]	684
	685	!
	686	!-- Close input file for virtual measurements. Therefore, just call
	687	!-- the read attribute routine with the "close" option.
	688	CALL netcdf_data_input_att( nvm, char_numstations, id_vm, '', &
	689	global_attribute, 'close', '' )
	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
	700	USE arrays_3d !, &
	701	! ONLY: pt
[3471]	702
	703	USE surface_mod
[3434]	704
	705	IMPLICIT NONE
	706
	707	CHARACTER(LEN=10) :: trimvar !< dummy for the measured variable name
	708
	709	INTEGER(iwp) :: l !<
	710	INTEGER(iwp) :: m !<
	711	INTEGER(iwp) :: var !<
[3471]	712
	713	INTEGER(iwp) :: mm, j, i
[3434]	714
[3471]	715	! write(9,*) "sampling"
	716	! flush(9)
[3434]	717	!
	718	!-- Loop over all stations. For each possible variable loop over all
	719	!-- observation points
	720	DO l = 1, nvm
	721	!
	722	!-- Loop over all measured variables. Please note, for the moment
	723	!-- the same indices for scalar and velocity components are used.
	724	!-- ToDo: Revise this later.
[3471]	725	! DO m = 1, vmea(l)%ns
	726	! j = vmea(l)%j(m)
	727	! i = vmea(l)%i(m)
	728	!
	729	! IF ( i >= nxl .AND. i <= nxr .AND. &
	730	! j >= nys .AND. j <= nyn ) THEN
	731	! IF ( surf_def_h(0)%start_index(j,i) <= &
	732	! surf_def_h(0)%end_index(j,i) ) THEN
	733	!
	734	! write(9,*) "sampled"
	735	! flush(9)
	736	! mm = surf_def_h(0)%end_index(j,i)
	737	!
	738	! surf_def_h(0)%pt_surface(mm) = -99.0
	739	! ENDIF
	740	! ENDIF
	741	! ENDDO
	742	! DO var = 1, vmea(l)%nvar
	743	! trimvar = TRIM( vmea(l)%measured_vars_name(var) )
	744	!
	745	! IF ( TRIM( trimvar ) == 'theta' ) THEN
	746	! DO m = 1, vmea(l)%ns
	747	! vmea(l)%measured_vars(var,m) = pt(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
	748	! ENDDO
	749	! ENDIF
	750	! IF ( TRIM( trimvar ) == 'w' ) THEN
	751	! DO m = 1, vmea(l)%ns
	752	! vmea(l)%measured_vars(var,m) = w(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
	753	! ENDDO
	754	! ENDIF
	755	! IF ( TRIM( trimvar ) == 'v' ) THEN
	756	! DO m = 1, vmea(l)%ns
	757	! vmea(l)%measured_vars(var,m) = v(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
	758	! ENDDO
	759	! ENDIF
	760	! IF ( TRIM( trimvar ) == 'u' ) THEN
	761	! DO m = 1, vmea(l)%ns
	762	! vmea(l)%measured_vars(var,m) = u(vmea(l)%k(m),vmea(l)%j(m),vmea(l)%i(m))
	763	! ENDDO
	764	! ENDIF
	765	! ENDDO
[3434]	766
	767	ENDDO
	768
[3471]	769	END SUBROUTINE vm_sampling
[3434]	770
	771
[3471]	772	END MODULE virtual_measurement_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |