Home

Context Navigation

source: palm/trunk/UTIL/agent_preprocessing/agent_preprocessing.f90 @ 4558

Last change on this file since 4558 was 4370, checked in by raasch, 5 years ago
bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added
Property svn:keywords set to `Id`
File size: 117.1 KB

Rev	Line
[3208]	1	!> @agent_preprocessing.f90
[3159]	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	!
[4370]	17	! Copyright 1997-2020 Leibniz Universitaet Hannover
[3159]	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! ------------------
	22	!
	23	!
	24	! Former revisions:
	25	! -----------------!
[3168]	26	! $Id: agent_preprocessing.f90 4370 2020-01-10 14:00:44Z moh.hefny $
[3400]	27	! Changed getcwd to GET_ENVIRONMENT_VARIABLE and changed X to 1X in format
	28	! statements to have code conform to fortran 2008 standard
	29	!
	30	! 3259 2018-09-18 09:53:18Z sward
[3259]	31	! Removed unused variables and fixed real to real comparisons
	32	!
	33	! 3216 2018-08-29 10:22:12Z sward
[3216]	34	! Bubfix for gfortran: reordering of type definitions
	35	!
	36	! 3210 2018-08-28 07:31:13Z sward
[3210]	37	! Bugfix: changed intrinsic SIZEOF to STORAGE_SIZE
	38	!
	39	! 3208 2018-08-27 13:10:50Z sward
[3208]	40	! Renamed nav_mesh to agent_preprocessing, adapted terminal output
	41	!
	42	! 3198 2018-08-15 09:23:10Z sward
[3198]	43	! Reduced tolerance_dp to 3 entries, fixed its initialization
	44	!
	45	! 3168 2018-07-25 06:40:29Z sward
[3168]	46	! Updated NetCDF ororgraphy and building input
	47	!
[3159]	48	! Initial revision
	49	!
	50	!
	51	! Description:
	52	! ------------
	53	!> Reads topography and building data and converts this information into a
	54	!> navigation mesh (NavMesh, a visibility graph). This mesh is necessary for
	55	!> the use of the Multi Agent System in PALM for the agents to navigate in
	56	!> complex (urban) terrain.
	57	!------------------------------------------------------------------------------!
	58
	59	MODULE kinds
	60
	61	IMPLICIT NONE
	62
	63	!
	64	!-- Floating point kinds
	65	INTEGER, PARAMETER :: sp = 4 !< single precision (32 bit)
	66	INTEGER, PARAMETER :: dp = 8 !< double precision (64 bit)
	67
	68	!
	69	!-- Integer kinds
	70	INTEGER, PARAMETER :: isp = SELECTED_INT_KIND( 9 ) !< single precision (32 bit)
	71	INTEGER, PARAMETER :: idp = SELECTED_INT_KIND( 14 ) !< double precision (64 bit)
	72
	73	!
	74	!-- Set kinds to be used as defaults
	75	INTEGER, PARAMETER :: wp = dp !< default real kind
	76	INTEGER, PARAMETER :: iwp = isp !< default integer kind
	77
	78	SAVE
	79
	80	END MODULE kinds
	81
	82	MODULE variables
	83
	84	USE kinds
	85
	86	CHARACTER(LEN=3) :: char_lod = 'lod' !< name of level-of-detail attribute in NetCDF file
	87	CHARACTER(LEN=10) :: char_fill = '_FillValue' !< name of fill value attribute in NetCDF file
	88	CHARACTER(LEN=128) :: runname !< Run name
[3168]	89
[3159]	90	LOGICAL :: internal_buildings = .FALSE. !< Flag that indicates whether buildings within closed courtyards should be deleted
[3168]	91	LOGICAL :: flag_2d = .FALSE. !< Flag that indicates that 2d buildings will be used in all cases
[3159]	92
	93	INTEGER(iwp) :: i !< Index along x
	94	INTEGER(iwp) :: j !< Index along y
	95	INTEGER(iwp) :: nx = 99999 !< Number of grid points in x-direction
	96	INTEGER(iwp) :: ny = 99999 !< Number of grid points in x-direction
	97	INTEGER(iwp) :: nov !< Number of vertices
	98	INTEGER(iwp) :: polygon_counter !< Iterator for the number of building polygons
	99	INTEGER(iwp) :: number_of_connections = 0 !< Counter for number of connections in mesh
	100	INTEGER(iwp) :: i_cn !< Min number of corners left in polygons after Douglas Poiker algorithm
	101	INTEGER(iwp) :: i_sc !< Cycle number for Douglas-Peucker algorithm
	102	INTEGER(iwp) :: nc_stat !< return value of nf90 function call
	103	INTEGER(iwp) :: vertex_counter !< Counter: total number of vertices
	104
	105	INTEGER, DIMENSION(:,:), ALLOCATABLE :: wall_flags_0 !< Bit-array containing surface information
	106	INTEGER, DIMENSION(:,:), ALLOCATABLE :: polygon_id !< Identifies each grid point as part of exactly one building
	107
	108	REAL(wp) :: ddx !< inverse of dx
	109	REAL(wp) :: ddy !< inverse of dy
	110	REAL(wp) :: dx = 99999.9_wp !< grid spacing in x-direction
	111	REAL(wp) :: dy = 99999.9_wp !< grid spacing in x-direction
	112	REAL(wp) :: dz = 99999.9_wp !< grid spacing in x-direction
	113	REAL(wp) :: finish !< variable for CPU time measurement
	114	REAL(wp) :: start !< variable for CPU time measurement
	115
[3198]	116	REAL(wp), DIMENSION(0:2) :: tolerance_dp = 999999.0_wp !< tolerance in Douglas-Peucker algorithm
[3159]	117
	118	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: obstacle_height !< height of obstacles
	119
	120	!
	121	!-- Define data structure where the dimension and type of the input depends
	122	!-- on the given level of detail.
	123	!-- For buildings, the input is either 2D float, or 3d byte.
	124	TYPE build_in
	125	INTEGER(iwp) :: lod = 1 !< level of detail
	126	INTEGER(KIND=1) :: fill2 = -127 !< fill value for lod = 2
	127	INTEGER(iwp) :: nz !< number of vertical layers in file
	128	INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: var_3d !< 3d variable (lod = 2)
	129	REAL(wp), DIMENSION(:), ALLOCATABLE :: z !< vertical coordinate for 3D building, used for consistency check
	130	LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
	131	REAL(wp) :: fill1 = -9999.9_wp !< fill values for lod = 1
	132	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2d !< 2d variable (lod = 1)
	133	END TYPE build_in
	134
	135	!
	136	!-- Topography grid point
	137	TYPE grid_point
	138	LOGICAL :: checked !< Flag to indicate whether this grid point has been evaluated already
	139	INTEGER(iwp) :: i !< x-index
	140	INTEGER(iwp) :: j !< y-index
	141	INTEGER(iwp) :: polygon_id !< ID of the polygon this grid point belongs to
	142	END TYPE grid_point
	143
	144	!
	145	!-- Node in the visibility graph navigation mesh
	146	TYPE mesh_point
	147	INTEGER(iwp) :: polygon_id !< Polygon the point belongs to
	148	INTEGER(iwp) :: vertex_id !< Vertex in the polygon
	149	INTEGER(iwp) :: noc !< number of connections
	150	INTEGER(iwp) :: origin_id !< ID of previous mesh point on path (A*)
	151	REAL(wp) :: cost_so_far !< Cost to reach this mesh point (A*)
	152	REAL(wp) :: x !< x-coordinate
	153	REAL(wp) :: y !< y-coordinate
	154	REAL(wp) :: x_s !< corner shifted outward from building by 1m (x)
	155	REAL(wp) :: y_s !< corner shifted outward from building by 1m (y)
	156	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: connected_vertices !< Index of connected vertices
	157	REAL(wp), DIMENSION(:), ALLOCATABLE :: distance_to_vertex !< Distance to each vertex
	158	END TYPE mesh_point
	159
	160	!
[3216]	161	!-- Vertex of a polygon
	162	TYPE vertex_type
	163	LOGICAL :: delete !< Flag to mark vertex for deletion
	164	REAL(wp) :: x !< x-coordinate
	165	REAL(wp) :: y !< y-coordinate
	166	END TYPE vertex_type
	167
	168	!
[3159]	169	!-- Polygon containing a number of vertices
	170	TYPE polygon_type
	171	INTEGER(iwp) :: nov !< Number of vertices in this polygon
	172	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: vertices !< Array of vertices
	173	END TYPE polygon_type
	174
	175	!
	176	!-- Define data type to read 2D real variables
	177	TYPE real_2d
	178	LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
	179	REAL(wp) :: fill = -9999.9_wp !< fill value
	180	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var !< respective variable
	181	END TYPE real_2d
	182
	183	!
	184	!-- Define data type to read 2D real variables
	185	TYPE real_3d
	186	LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
	187	INTEGER(iwp) :: nz !< number of grid points along vertical dimension
	188	REAL(wp) :: fill = -9999.9_wp !< fill value
	189	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: var !< respective variable
	190	END TYPE real_3d
	191
	192	TYPE(grid_point), DIMENSION(:,:), ALLOCATABLE :: grid !< 2d Topography grid
	193
	194	TYPE(mesh_point), DIMENSION(:), ALLOCATABLE, TARGET :: mesh !< Navigation mesh
	195
	196	TYPE(real_2d) :: terrain_height_f !< input variable for terrain height
	197
	198	TYPE(vertex_type) :: dummy_vertex !< placeholder vertex used for data copying
	199	TYPE(vertex_type) :: null_vertex !< placeholder vertex used for initialisation
	200
	201	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: dummy_v_list !< Dummy for reallocation of polygon array
	202
	203	TYPE(polygon_type), POINTER :: polygon !< Current polygon
	204
	205	TYPE(polygon_type), DIMENSION(:), ALLOCATABLE, TARGET :: polygons !< Building polygons
	206
	207	TYPE(build_in) :: buildings_f !< input variable for buildings
	208
	209
	210	END MODULE variables
	211
	212	MODULE data_input
	213
	214	USE kinds
	215	USE variables
	216
	217	#if defined ( __netcdf )
	218	USE NETCDF
	219	#endif
	220
	221	INTERFACE get_variable
	222	MODULE PROCEDURE get_variable_1d_int
	223	MODULE PROCEDURE get_variable_1d_real
	224	MODULE PROCEDURE get_variable_2d_int8
	225	MODULE PROCEDURE get_variable_2d_int32
	226	MODULE PROCEDURE get_variable_2d_real
	227	MODULE PROCEDURE get_variable_3d_int8
	228	MODULE PROCEDURE get_variable_3d_real
	229	MODULE PROCEDURE get_variable_4d_real
	230	END INTERFACE get_variable
	231
	232	INTERFACE get_attribute
	233	MODULE PROCEDURE get_attribute_real
	234	MODULE PROCEDURE get_attribute_int8
	235	MODULE PROCEDURE get_attribute_int32
	236	MODULE PROCEDURE get_attribute_string
	237	END INTERFACE get_attribute
	238
	239	CONTAINS
	240
	241	!------------------------------------------------------------------------------!
	242	! Description:
	243	! ------------
	244	!> Reads orography and building information.
	245	!------------------------------------------------------------------------------!
	246	SUBROUTINE netcdf_data_input_topo ( input_trunk )
	247
	248	IMPLICIT NONE
	249
	250	CHARACTER(LEN=*) :: input_trunk !< run path
	251	CHARACTER(LEN=200) :: input_filename !< filename
	252
	253	CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE :: var_names !< variable names in static input file
	254
	255
	256	INTEGER(iwp) :: i !< running index along x-direction
	257	INTEGER(iwp) :: k_head !< minimum k index for agents to walk underneath overhanging buildings
	258	INTEGER(iwp) :: id_topo !< NetCDF id of topograhy input file
	259	INTEGER(iwp) :: j !< running index along y-direction
	260	INTEGER(iwp) :: num_vars !< number of variables in netcdf input file
	261
	262	LOGICAL :: netcdf_flag = .FALSE. !< indicates whether netcdf file is used for input
	263	LOGICAL :: lod_flag = .FALSE. !< true if 3d building data is used
	264	LOGICAL :: topo_file_flag = .FALSE. !< true if 3d building data is used
	265
[3400]	266	WRITE(*,'((1X,A,/))') 'Looking for topography/building information'
[3159]	267	INQUIRE( FILE = TRIM( input_trunk )//'_static', EXIST = netcdf_flag )
	268
	269	IF ( netcdf_flag ) THEN
	270	input_filename = TRIM( input_trunk )//'_static'
	271	WRITE(*,'(2(3X,A,/))') 'Topography/building data will be used from', &
	272	TRIM( input_trunk )//'_static'
	273	ELSE
	274	WRITE(*,'(2(3X,A,/))') 'No static driver was found.', &
	275	'Trying to read building data from _topo file.'
	276	input_filename = TRIM( input_trunk )//'_topo'
	277	INQUIRE( FILE = TRIM( input_filename ), EXIST = topo_file_flag )
	278	IF ( .NOT. topo_file_flag ) THEN
	279	WRITE(*,'(6(3X,A,/))') &
	280	'No ASCII topography file was found in INPUT directory.', &
	281	'Make sure you provided building data in the form of either', &
	282	' (A) a static driver (<runname>_static) or', &
	283	' (B) an ASCII topography file (<runname>_topo).', &
	284	NEW_LINE('A')//'Aborting nav_mesh program...'
	285	STOP
	286	ENDIF
	287	WRITE(*,'(2(3X,A,/))') 'Topography/building data will be used from', &
	288	TRIM( input_trunk )//'_topo'
	289	ENDIF
	290
	291	!
	292	!-- Input via palm-input data standard
	293	IF ( netcdf_flag ) THEN
	294	#if defined ( __netcdf )
	295	!
	296	!-- Open file in read-only mode
	297	CALL open_read_file( TRIM(input_filename) , id_topo )
	298
	299	!
	300	!-- At first, inquire all variable names.
	301	!-- This will be used to check whether an input variable exist or not.
	302	nc_stat = NF90_INQUIRE( id_topo, NVARIABLES = num_vars )
	303	CALL handle_error( 'inquire_num_variables', 534 )
	304	!
	305	!-- Allocate memory to store variable names and inquire them.
	306	ALLOCATE( var_names(1:num_vars) )
	307	CALL inquire_variable_names( id_topo, var_names )
	308	!
	309	!-- Terrain height. First, get variable-related _FillValue attribute
[3168]	310	IF ( check_existence( var_names, 'zt' ) ) THEN
[3159]	311	terrain_height_f%from_file = .TRUE.
	312	CALL get_attribute( id_topo, char_fill, &
	313	terrain_height_f%fill, &
[3168]	314	.FALSE., 'zt' )
[3159]	315	!
	316	!-- PE-wise reading of 2D terrain height.
	317	ALLOCATE ( terrain_height_f%var(0:ny,0:nx) )
	318	DO i = 0, nx
[3168]	319	CALL get_variable( id_topo, 'zt', &
[3159]	320	i, terrain_height_f%var(:,i) )
	321	ENDDO
	322	ELSE
	323	terrain_height_f%from_file = .FALSE.
	324	ENDIF
	325
	326	!
	327	!-- Read building height. First, read its _FillValue attribute,
	328	!-- as well as lod attribute
[3168]	329	buildings_f%from_file = .FALSE.
	330	IF ( check_existence( var_names, 'buildings_2d' ) ) THEN
[3159]	331	buildings_f%from_file = .TRUE.
	332	CALL get_attribute( id_topo, char_lod, buildings_f%lod, &
[3168]	333	.FALSE., 'buildings_2d' )
[3159]	334
	335	CALL get_attribute( id_topo, char_fill, &
	336	buildings_f%fill1, &
[3168]	337	.FALSE., 'buildings_2d' )
[3159]	338
	339	!
	340	!-- Read 2D topography
	341	IF ( buildings_f%lod == 1 ) THEN
	342	ALLOCATE ( buildings_f%var_2d(0:ny,0:nx) )
	343	DO i = 0, nx
[3168]	344	CALL get_variable( id_topo, 'buildings_2d', &
[3159]	345	i, buildings_f%var_2d(:,i) )
	346	ENDDO
	347	ELSE
	348	WRITE(*,'(A)') 'NetCDF attribute lod ' // &
	349	'(level of detail) is not set properly.'
	350	ENDIF
	351	ENDIF
	352	!
	353	!-- If available, also read 3D building information. If both are
[3168]	354	!-- available, use 3D information. Do this only if the flag that indicates
	355	!-- that 2d buildings shall be used no matter what is false.
	356	IF ( check_existence( var_names, 'buildings_3d' ) &
	357	.AND. .NOT. flag_2d ) &
	358	THEN
[3159]	359	lod_flag = .TRUE.
	360	buildings_f%from_file = .TRUE.
	361	CALL get_attribute( id_topo, char_lod, buildings_f%lod, &
[3168]	362	.FALSE., 'buildings_3d' )
[3159]	363
	364	CALL get_attribute( id_topo, char_fill, &
	365	buildings_f%fill2, &
[3168]	366	.FALSE., 'buildings_3d' )
[3159]	367
	368	CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
	369
	370	IF ( buildings_f%lod == 2 ) THEN
	371	ALLOCATE( buildings_f%z(0:buildings_f%nz-1) )
	372	CALL get_variable( id_topo, 'z', buildings_f%z )
	373
	374	ALLOCATE( buildings_f%var_3d(0:buildings_f%nz-1, &
	375	0:ny,0:nx) )
	376	buildings_f%var_3d = 0
	377	!
	378	!-- Read data PE-wise. Read yz-slices.
	379	DO i = 0, nx
	380	DO j = 0, ny
[3168]	381	CALL get_variable( id_topo, 'buildings_3d', &
[3159]	382	i, j, &
	383	buildings_f%var_3d(:,j,i) )
	384	ENDDO
	385	ENDDO
	386	ELSE
	387	WRITE(*,'(A)') 'NetCDF attribute lod ' // &
	388	'(level of detail) is not set properly.'
	389	ENDIF
	390	ENDIF
	391
	392	!
	393	!-- Close topography input file
	394	CALL close_input_file( id_topo )
	395	#endif
	396	!
	397	!-- ASCII input
	398	ELSE
	399
	400	OPEN( 90, FILE= input_filename, &
	401	STATUS='OLD', FORM='FORMATTED' )
	402	!
	403	!-- Read data from nyn to nys and nxl to nxr. Therefore, skip
	404	!-- column until nxl-1 is reached
	405	ALLOCATE ( buildings_f%var_2d(0:ny,0:nx) )
	406	DO j = ny, 0, -1
	407	READ( 90, *, ERR=11, END=11 ) &
	408	( buildings_f%var_2d(j,i), i = 0, nx )
	409	ENDDO
	410
	411	GOTO 12
	412
	413	11 WRITE(*,'(2A)') 'errors in file ',input_filename
	414
	415	12 CLOSE( 90 )
	416	buildings_f%from_file = .TRUE.
	417
	418	ENDIF
	419	!
	420	!-- In case no terrain height is provided by static input file, allocate
	421	!-- array nevertheless and set terrain height to 0, which simplifies
	422	!-- topography initialization.
	423	IF ( .NOT. terrain_height_f%from_file ) THEN
	424	ALLOCATE ( terrain_height_f%var(0:ny,0:nx) )
	425	terrain_height_f%var = 0.0_wp
	426	ENDIF
	427
	428	!
	429	!-- Transfer read data to uniform format: For agents the only relevant
	430	!-- information is whether they can walk or not at ground level.
	431	k_head = CEILING(2./dz)
	432	IF ( buildings_f%from_file ) THEN
	433	IF ( lod_flag ) THEN
	434	obstacle_height(0:nx,0:ny) = 1.
	435	DO j = 0, ny
	436	DO i = 0, nx
	437	!
	438	!-- For this purpose, an overhanging structure that an angent
	439	!-- can walk beneath (e.g. a doorway) is not considered an
	440	!-- obstacle.
	441	IF ( ALL( buildings_f%var_3d(0:k_head,j,i) == 0 ) ) THEN
	442	obstacle_height(i,j) = 0.
	443	ENDIF
	444	ENDDO
	445	ENDDO
	446	ELSE
	447	DO j = 0, ny
	448	DO i = 0, nx
	449	obstacle_height(i,j) = buildings_f%var_2d(j,i)
	450	ENDDO
	451	ENDDO
	452	ENDIF
	453	ELSE
	454	WRITE(,) 'No building data was read from file. There will be no' //&
	455	'navigation data available to agents.'
	456	ENDIF
	457
	458	END SUBROUTINE netcdf_data_input_topo
	459
	460	!------------------------------------------------------------------------------!
	461	! Description:
	462	! ------------
	463	!> Checks if a given variables is on file
	464	!------------------------------------------------------------------------------!
	465	FUNCTION check_existence( vars_in_file, var_name )
	466
	467	IMPLICIT NONE
	468
	469	CHARACTER(LEN=*) :: var_name !< variable to be checked
	470	CHARACTER(LEN=*), DIMENSION(:) :: vars_in_file !< list of variables in file
	471
	472	INTEGER(iwp) :: i !< loop variable
	473
	474	LOGICAL :: check_existence !< flag indicating whether a variable exist or not - actual return value
	475
	476	i = 1
	477	check_existence = .FALSE.
	478	DO WHILE ( i <= SIZE( vars_in_file ) )
	479	check_existence = TRIM( vars_in_file(i) ) == TRIM( var_name ) .OR. &
	480	check_existence
	481	i = i + 1
	482	ENDDO
	483
	484	RETURN
	485
	486	END FUNCTION check_existence
	487
	488
	489	!------------------------------------------------------------------------------!
	490	! Description:
	491	! ------------
	492	!> Closes an existing netCDF file.
	493	!------------------------------------------------------------------------------!
	494	SUBROUTINE close_input_file( id )
	495	#if defined( __netcdf )
	496
	497	IMPLICIT NONE
	498
	499	INTEGER(iwp), INTENT(INOUT) :: id !< file id
	500
	501	nc_stat = NF90_CLOSE( id )
	502	CALL handle_error( 'close', 537 )
	503	#endif
	504	END SUBROUTINE close_input_file
	505
	506	!------------------------------------------------------------------------------!
	507	! Description:
	508	! ------------
	509	!> Opens an existing netCDF file for reading only and returns its id.
	510	!------------------------------------------------------------------------------!
	511	SUBROUTINE open_read_file( filename, id )
	512	#if defined( __netcdf )
	513
	514	IMPLICIT NONE
	515
	516	CHARACTER (LEN=*), INTENT(IN) :: filename !< filename
	517	INTEGER(iwp), INTENT(INOUT) :: id !< file id
	518
	519	nc_stat = NF90_OPEN( filename, NF90_NOWRITE, id )
	520
	521	CALL handle_error( 'open_read_file', 536 )
	522
	523	#endif
	524	END SUBROUTINE open_read_file
	525
	526
	527
	528	!------------------------------------------------------------------------------!
	529	! Description:
	530	! ------------
	531	!> Get dimension array for a given dimension
	532	!------------------------------------------------------------------------------!
	533	SUBROUTINE get_dimension_length( id, dim_len, variable_name )
	534	#if defined( __netcdf )
	535
	536	IMPLICIT NONE
	537
	538	CHARACTER(LEN=*) :: variable_name !< dimension name
	539	CHARACTER(LEN=100) :: dum !< dummy variable to receive return character
	540
	541	INTEGER(iwp) :: dim_len !< dimension size
	542	INTEGER(iwp), INTENT(IN) :: id !< file id
	543	INTEGER(iwp) :: id_dim !< dimension id
	544
	545	!
	546	!-- First, inquire dimension ID
	547	nc_stat = NF90_INQ_DIMID( id, TRIM( variable_name ), id_dim )
	548	CALL handle_error( 'get_dimension_length', 526 )
	549	!
	550	!-- Inquire dimension length
	551	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, dum, LEN = dim_len )
	552	CALL handle_error( 'get_dimension_length', 526 )
	553
	554	#endif
	555	END SUBROUTINE get_dimension_length
	556
	557	!------------------------------------------------------------------------------!
	558	! Description:
	559	! ------------
	560	!> Reads a 1D integer variable from file.
	561	!------------------------------------------------------------------------------!
	562	SUBROUTINE get_variable_1d_int( id, variable_name, var )
	563
	564	IMPLICIT NONE
	565
	566	CHARACTER(LEN=*) :: variable_name !< variable name
	567
	568	INTEGER(iwp), INTENT(IN) :: id !< file id
	569	INTEGER(iwp) :: id_var !< dimension id
	570
	571	INTEGER(iwp), DIMENSION(:), INTENT(INOUT) :: var !< variable to be read
	572	#if defined( __netcdf )
	573
	574	!
	575	!-- First, inquire variable ID
	576	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	577	CALL handle_error( 'get_variable_1d_int', 527 )
	578	!
	579	!-- Inquire dimension length
	580	nc_stat = NF90_GET_VAR( id, id_var, var )
	581	CALL handle_error( 'get_variable_1d_int', 527 )
	582
	583	#endif
	584	END SUBROUTINE get_variable_1d_int
	585
	586	!------------------------------------------------------------------------------!
	587	! Description:
	588	! ------------
	589	!> Reads a 1D float variable from file.
	590	!------------------------------------------------------------------------------!
	591	SUBROUTINE get_variable_1d_real( id, variable_name, var )
	592
	593	IMPLICIT NONE
	594
	595	CHARACTER(LEN=*) :: variable_name !< variable name
	596
	597	INTEGER(iwp), INTENT(IN) :: id !< file id
	598	INTEGER(iwp) :: id_var !< dimension id
	599
	600	REAL(wp), DIMENSION(:), INTENT(INOUT) :: var !< variable to be read
	601	#if defined( __netcdf )
	602
	603	!
	604	!-- First, inquire variable ID
	605	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	606	CALL handle_error( 'get_variable_1d_real', 527 )
	607	!
	608	!-- Inquire dimension length
	609	nc_stat = NF90_GET_VAR( id, id_var, var )
	610	CALL handle_error( 'get_variable_1d_real', 527 )
	611
	612	#endif
	613	END SUBROUTINE get_variable_1d_real
	614
	615	!------------------------------------------------------------------------------!
	616	! Description:
	617	! ------------
	618	!> Reads a 2D REAL variable from a file. Reading is done processor-wise,
	619	!> i.e. each core reads its own domain in slices along x.
	620	!------------------------------------------------------------------------------!
	621	SUBROUTINE get_variable_2d_real( id, variable_name, i, var )
	622
	623	IMPLICIT NONE
	624
	625	CHARACTER(LEN=*) :: variable_name !< variable name
	626
	627	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	628	INTEGER(iwp), INTENT(IN) :: id !< file id
	629	INTEGER(iwp) :: id_var !< variable id
	630
	631	REAL(wp), DIMENSION(0:ny), INTENT(INOUT) :: var !< variable to be read
	632	#if defined( __netcdf )
	633	!
	634	!-- Inquire variable id
	635	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	636	!
	637	!-- Get variable
	638	nc_stat = NF90_GET_VAR( id, id_var, var(0:ny), &
	639	start = (/ i+1, 1 /), &
	640	count = (/ 1, ny + 1 /) )
	641
	642	CALL handle_error( 'get_variable_2d_real', 528 )
	643	#endif
	644	END SUBROUTINE get_variable_2d_real
	645
	646	!------------------------------------------------------------------------------!
	647	! Description:
	648	! ------------
	649	!> Reads a 2D 32-bit INTEGER variable from file. Reading is done processor-wise,
	650	!> i.e. each core reads its own domain in slices along x.
	651	!------------------------------------------------------------------------------!
	652	SUBROUTINE get_variable_2d_int32( id, variable_name, i, var )
	653
	654	IMPLICIT NONE
	655
	656	CHARACTER(LEN=*) :: variable_name !< variable name
	657
	658	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	659	INTEGER(iwp), INTENT(IN) :: id !< file id
	660	INTEGER(iwp) :: id_var !< variable id
	661	INTEGER(iwp), DIMENSION(0:ny), INTENT(INOUT) :: var !< variable to be read
	662	#if defined( __netcdf )
	663	!
	664	!-- Inquire variable id
	665	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	666	!
	667	!-- Get variable
	668	nc_stat = NF90_GET_VAR( id, id_var, var(0:ny), &
	669	start = (/ i+1, 1 /), &
	670	count = (/ 1, ny + 1 /) )
	671
	672	CALL handle_error( 'get_variable_2d_int32', 529 )
	673	#endif
	674	END SUBROUTINE get_variable_2d_int32
	675
	676	!------------------------------------------------------------------------------!
	677	! Description:
	678	! ------------
	679	!> Reads a 2D 8-bit INTEGER variable from file. Reading is done processor-wise,
	680	!> i.e. each core reads its own domain in slices along x.
	681	!------------------------------------------------------------------------------!
	682	SUBROUTINE get_variable_2d_int8( id, variable_name, i, var )
	683
	684	IMPLICIT NONE
	685
	686	CHARACTER(LEN=*) :: variable_name !< variable name
	687
	688	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	689	INTEGER(iwp), INTENT(IN) :: id !< file id
	690	INTEGER(iwp) :: id_var !< variable id
	691	INTEGER(KIND=1), DIMENSION(0:ny), INTENT(INOUT) :: var !< variable to be read
	692	#if defined( __netcdf )
	693	!
	694	!-- Inquire variable id
	695	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	696	!
	697	!-- Get variable
	698	nc_stat = NF90_GET_VAR( id, id_var, var(0:ny), &
	699	start = (/ i+1, 1 /), &
	700	count = (/ 1, ny + 1 /) )
	701
	702	CALL handle_error( 'get_variable_2d_int8', 530 )
	703	#endif
	704	END SUBROUTINE get_variable_2d_int8
	705
	706	!------------------------------------------------------------------------------!
	707	! Description:
	708	! ------------
	709	!> Reads a 3D 8-bit INTEGER variable from file.
	710	!------------------------------------------------------------------------------!
	711	SUBROUTINE get_variable_3d_int8( id, variable_name, i, j, var )
	712
	713	IMPLICIT NONE
	714
	715	CHARACTER(LEN=*) :: variable_name !< variable name
	716
	717	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	718	INTEGER(iwp), INTENT(IN) :: id !< file id
	719	INTEGER(iwp) :: id_var !< variable id
	720	INTEGER(iwp), INTENT(IN) :: j !< index along y direction
	721	INTEGER(iwp) :: n_file !< number of data-points along 3rd dimension
	722
	723	INTEGER(iwp), DIMENSION(1:3) :: id_dim
	724
	725	INTEGER( KIND = 1 ), DIMENSION(0:buildings_f%nz-1), INTENT(INOUT) :: var !< variable to be read
	726	#if defined( __netcdf )
	727
	728	!
	729	!-- Inquire variable id
	730	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	731	!
	732	!-- Get length of first dimension, required for the count parameter.
	733	!-- Therefore, first inquired dimension ids
	734	nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
	735	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n_file )
	736	!
	737	!-- Get variable
	738	nc_stat = NF90_GET_VAR( id, id_var, var, &
	739	start = (/ i+1, j+1, 1 /), &
	740	count = (/ 1, 1, n_file /) )
	741
	742	CALL handle_error( 'get_variable_3d_int8', 531 )
	743	#endif
	744	END SUBROUTINE get_variable_3d_int8
	745
	746
	747	!------------------------------------------------------------------------------!
	748	! Description:
	749	! ------------
	750	!> Reads a 3D float variable from file.
	751	!------------------------------------------------------------------------------!
	752	SUBROUTINE get_variable_3d_real( id, variable_name, i, j, var )
	753
	754	IMPLICIT NONE
	755
	756	CHARACTER(LEN=*) :: variable_name !< variable name
	757
	758	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	759	INTEGER(iwp), INTENT(IN) :: id !< file id
	760	INTEGER(iwp) :: id_var !< variable id
	761	INTEGER(iwp), INTENT(IN) :: j !< index along y direction
	762	INTEGER(iwp) :: n3 !< number of data-points along 3rd dimension
	763
	764	INTEGER(iwp), DIMENSION(3) :: id_dim
	765
	766	REAL(wp), DIMENSION(:), INTENT(INOUT) :: var !< variable to be read
	767	#if defined( __netcdf )
	768
	769	!
	770	!-- Inquire variable id
	771	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	772	!
	773	!-- Get length of first dimension, required for the count parameter.
	774	!-- Therefore, first inquired dimension ids
	775	nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
	776	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
	777	!
	778	!-- Get variable
	779	nc_stat = NF90_GET_VAR( id, id_var, var, &
	780	start = (/ i+1, j+1, 1 /), &
	781	count = (/ 1, 1, n3 /) )
	782
	783	CALL handle_error( 'get_variable_3d_real', 532 )
	784	#endif
	785	END SUBROUTINE get_variable_3d_real
	786
	787
	788	!------------------------------------------------------------------------------!
	789	! Description:
	790	! ------------
	791	!> Reads a 4D float variable from file. Note, in constrast to 3D versions,
	792	!> dimensions are already inquired and passed so that they are known here.
	793	!------------------------------------------------------------------------------!
	794	SUBROUTINE get_variable_4d_real( id, variable_name, i, j, var, n3, n4 )
	795
	796	IMPLICIT NONE
	797
	798	CHARACTER(LEN=*) :: variable_name !< variable name
	799
	800	INTEGER(iwp), INTENT(IN) :: i !< index along x direction
	801	INTEGER(iwp), INTENT(IN) :: id !< file id
	802	INTEGER(iwp) :: id_var !< variable id
	803	INTEGER(iwp), INTENT(IN) :: j !< index along y direction
	804	INTEGER(iwp), INTENT(IN) :: n3 !< number of data-points along 3rd dimension
	805	INTEGER(iwp), INTENT(IN) :: n4 !< number of data-points along 4th dimension
	806
	807	REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: var !< variable to be read
	808	#if defined( __netcdf )
	809
	810	!
	811	!-- Inquire variable id
	812	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	813	!
	814	!-- Get variable
	815	nc_stat = NF90_GET_VAR( id, id_var, var, &
	816	start = (/ i+1, j+1, 1, 1 /), &
	817	count = (/ 1, 1, n3, n4 /) )
	818
	819	CALL handle_error( 'get_variable_4d_real', 533 )
	820	#endif
	821	END SUBROUTINE get_variable_4d_real
	822
	823	!------------------------------------------------------------------------------!
	824	! Description:
	825	! ------------
	826	!> Prints out a text message corresponding to the current status.
	827	!------------------------------------------------------------------------------!
	828	SUBROUTINE handle_error( routine_name, errno )
	829
	830	IMPLICIT NONE
	831
	832	CHARACTER(LEN=6) :: message_identifier
	833	CHARACTER(LEN=*) :: routine_name
	834	CHARACTER(LEN=100) :: message_string
	835
	836	INTEGER(iwp) :: errno
	837	#if defined( __netcdf )
	838
	839	IF ( nc_stat /= NF90_NOERR ) THEN
	840
	841	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
	842
	843	message_string = TRIM( NF90_STRERROR( nc_stat ) )
	844
	845	WRITE(,) routine_name,' ', message_identifier,' ', TRIM(message_string)
	846	WRITE(,) 'Aborting NavMesh-tool'
	847
	848	ENDIF
	849
	850	#endif
	851	END SUBROUTINE handle_error
	852
	853
	854	!------------------------------------------------------------------------------!
	855	! Description:
	856	! ------------
	857	!> Inquires the variable names belonging to a file.
	858	!------------------------------------------------------------------------------!
	859	SUBROUTINE inquire_variable_names( id, var_names )
	860
	861	IMPLICIT NONE
	862
	863	CHARACTER(LEN=*), DIMENSION(:), INTENT(INOUT) :: var_names !< return variable - variable names
	864	INTEGER(iwp) :: i !< loop variable
	865	INTEGER(iwp), INTENT(IN) :: id !< file id
	866	INTEGER(iwp) :: num_vars !< number of variables (unused return parameter)
	867	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: varids !< dummy array to strore variable ids temporarily
	868	#if defined( __netcdf )
	869
	870	ALLOCATE( varids(1:SIZE(var_names)) )
	871	nc_stat = NF90_INQ_VARIDS( id, NVARS = num_vars, VARIDS = varids )
	872	CALL handle_error( 'inquire_variable_names', 535 )
	873
	874	DO i = 1, SIZE(var_names)
	875	nc_stat = NF90_INQUIRE_VARIABLE( id, varids(i), NAME = var_names(i) )
	876	CALL handle_error( 'inquire_variable_names', 535 )
	877	ENDDO
	878
	879	DEALLOCATE( varids )
	880	#endif
	881	END SUBROUTINE inquire_variable_names
	882
	883	!------------------------------------------------------------------------------!
	884	! Description:
	885	! ------------
	886	!> Reads global or variable-related attributes of type INTEGER (32-bit)
	887	!------------------------------------------------------------------------------!
	888	SUBROUTINE get_attribute_int32( id, attribute_name, value, global, &
	889	variable_name )
	890
	891	IMPLICIT NONE
	892
	893	CHARACTER(LEN=*) :: attribute_name !< attribute name
	894	CHARACTER(LEN=*), OPTIONAL :: variable_name !< variable name
	895
	896	INTEGER(iwp), INTENT(IN) :: id !< file id
	897	INTEGER(iwp) :: id_var !< variable id
	898	INTEGER(iwp), INTENT(INOUT) :: value !< read value
	899
	900	LOGICAL, INTENT(IN) :: global !< flag indicating global attribute
	901	#if defined( __netcdf )
	902
	903	!
	904	!-- Read global attribute
	905	IF ( global ) THEN
	906	nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
	907	CALL handle_error( 'get_attribute_int32 global', 522 )
	908	!
	909	!-- Read attributes referring to a single variable. Therefore, first inquire
	910	!-- variable id
	911	ELSE
	912	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	913	CALL handle_error( 'get_attribute_int32', 522 )
	914	nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
	915	CALL handle_error( 'get_attribute_int32', 522 )
	916	ENDIF
	917	#endif
	918	END SUBROUTINE get_attribute_int32
	919
	920	!------------------------------------------------------------------------------!
	921	! Description:
	922	! ------------
	923	!> Reads global or variable-related attributes of type INTEGER (8-bit)
	924	!------------------------------------------------------------------------------!
	925	SUBROUTINE get_attribute_int8( id, attribute_name, value, global, &
	926	variable_name )
	927
	928	IMPLICIT NONE
	929
	930	CHARACTER(LEN=*) :: attribute_name !< attribute name
	931	CHARACTER(LEN=*), OPTIONAL :: variable_name !< variable name
	932
	933	INTEGER(iwp), INTENT(IN) :: id !< file id
	934	INTEGER(iwp) :: id_var !< variable id
	935	INTEGER(KIND=1), INTENT(INOUT) :: value !< read value
	936
	937	LOGICAL, INTENT(IN) :: global !< flag indicating global attribute
	938	#if defined( __netcdf )
	939
	940	!
	941	!-- Read global attribute
	942	IF ( global ) THEN
	943	nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
	944	CALL handle_error( 'get_attribute_int8 global', 523 )
	945	!
	946	!-- Read attributes referring to a single variable. Therefore, first inquire
	947	!-- variable id
	948	ELSE
	949	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	950	CALL handle_error( 'get_attribute_int8', 523 )
	951	nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
	952	CALL handle_error( 'get_attribute_int8', 523 )
	953	ENDIF
	954	#endif
	955	END SUBROUTINE get_attribute_int8
	956
	957	!------------------------------------------------------------------------------!
	958	! Description:
	959	! ------------
	960	!> Reads global or variable-related attributes of type REAL
	961	!------------------------------------------------------------------------------!
	962	SUBROUTINE get_attribute_real( id, attribute_name, value, global, &
	963	variable_name )
	964
	965	IMPLICIT NONE
	966
	967	CHARACTER(LEN=*) :: attribute_name !< attribute name
	968	CHARACTER(LEN=*), OPTIONAL :: variable_name !< variable name
	969
	970	INTEGER(iwp), INTENT(IN) :: id !< file id
	971	INTEGER(iwp) :: id_var !< variable id
	972
	973	LOGICAL, INTENT(IN) :: global !< flag indicating global attribute
	974
	975	REAL(wp), INTENT(INOUT) :: value !< read value
	976	#if defined( __netcdf )
	977
	978
	979	!
	980	!-- Read global attribute
	981	IF ( global ) THEN
	982	nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
	983	CALL handle_error( 'get_attribute_real global', 524 )
	984	!
	985	!-- Read attributes referring to a single variable. Therefore, first inquire
	986	!-- variable id
	987	ELSE
	988	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	989	CALL handle_error( 'get_attribute_real', 524 )
	990	nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
	991	CALL handle_error( 'get_attribute_real', 524 )
	992	ENDIF
	993	#endif
	994	END SUBROUTINE get_attribute_real
	995
	996	!------------------------------------------------------------------------------!
	997	! Description:
	998	! ------------
	999	!> Reads global or variable-related attributes of type CHARACTER
	1000	!> Remark: reading attributes of type NF_STRING return an error code 56 -
	1001	!> Attempt to convert between text & numbers.
	1002	!------------------------------------------------------------------------------!
	1003	SUBROUTINE get_attribute_string( id, attribute_name, value, global, &
	1004	variable_name )
	1005
	1006	IMPLICIT NONE
	1007
	1008	CHARACTER(LEN=*) :: attribute_name !< attribute name
	1009	CHARACTER(LEN=*), OPTIONAL :: variable_name !< variable name
	1010	CHARACTER(LEN=*), INTENT(INOUT) :: value !< read value
	1011
	1012	INTEGER(iwp), INTENT(IN) :: id !< file id
	1013	INTEGER(iwp) :: id_var !< variable id
	1014
	1015	LOGICAL, INTENT(IN) :: global !< flag indicating global attribute
	1016	#if defined( __netcdf )
	1017
	1018	!
	1019	!-- Read global attribute
	1020	IF ( global ) THEN
	1021	nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
	1022	CALL handle_error( 'get_attribute_string global', 525 )
	1023	!
	1024	!-- Read attributes referring to a single variable. Therefore, first inquire
	1025	!-- variable id
	1026	ELSE
	1027	nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
	1028	CALL handle_error( 'get_attribute_string', 525 )
	1029
	1030	nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
	1031	CALL handle_error( 'get_attribute_string',525 )
	1032	ENDIF
	1033	#endif
	1034	END SUBROUTINE get_attribute_string
	1035
	1036	END MODULE
	1037
	1038	MODULE mod_functions
	1039
	1040	USE kinds
	1041
	1042	PRIVATE
	1043	PUBLIC dist_point_to_edge, intersect, is_left, is_right
	1044
	1045	CONTAINS
	1046
	1047	!
	1048	!-- Calculates distance of point P to edge (A,B). If A = B, calculates
	1049	!-- point-to-point distance from A/B to P
	1050	FUNCTION dist_point_to_edge ( a_x, a_y, b_x, b_y, p_x, p_y )
	1051
	1052	IMPLICIT NONE
	1053
	1054	REAL(wp) :: ab_x !< x-coordinate of vector from A to B
	1055	REAL(wp) :: ab_y !< y-coordinate of vector from A to B
	1056	REAL(wp) :: ab_d !< inverse length of vector from A to B
	1057	REAL(wp) :: ab_u_x !< x-coordinate of vector with direction of ab and length 1
	1058	REAL(wp) :: ab_u_y !< y-coordinate of vector with direction of ab and length 1
	1059	REAL(wp) :: ba_x !< x-coordinate of vector from B to A
	1060	REAL(wp) :: ba_y !< y-coordinate of vector from B to A
	1061	REAL(wp) :: ap_x !< x-coordinate of vector from A to P
	1062	REAL(wp) :: ap_y !< y-coordinate of vector from A to P
	1063	REAL(wp) :: bp_x !< x-coordinate of vector from B to P
	1064	REAL(wp) :: bp_y !< y-coordinate of vector from B to P
	1065	REAL(wp) :: a_x !< x-coordinate of point A of edge
	1066	REAL(wp) :: a_y !< y-coordinate of point A of edge
	1067	REAL(wp) :: b_x !< x-coordinate of point B of edge
	1068	REAL(wp) :: b_y !< y-coordinate of point B of edge
	1069	REAL(wp) :: p_x !< x-coordinate of point P
	1070	REAL(wp) :: p_y !< y-coordinate of point P
	1071	REAL(wp) :: dist_x !< x-coordinate of point P
	1072	REAL(wp) :: dist_y !< y-coordinate of point P
	1073	REAL(wp) :: dist_point_to_edge !< y-coordinate of point P
	1074
	1075	ab_x = - a_x + b_x
	1076	ab_y = - a_y + b_y
	1077	ba_x = - b_x + a_x
	1078	ba_y = - b_y + a_y
	1079	ap_x = - a_x + p_x
	1080	ap_y = - a_y + p_y
	1081	bp_x = - b_x + p_x
	1082	bp_y = - b_y + p_y
	1083
	1084	IF ( ab_x * ap_x + ab_y * ap_y <= 0. ) THEN
	1085	dist_point_to_edge = SQRT((a_x - p_x)2 + (a_y - p_y)2)
	1086	ELSEIF ( ba_x * bp_x + ba_y * bp_y <= 0. ) THEN
	1087	dist_point_to_edge = SQRT((b_x - p_x)2 + (b_y - p_y)2)
	1088	ELSE
	1089	ab_d = 1./SQRT((ab_x)2+(ab_y)2)
	1090	ab_u_x = ab_x*ab_d
	1091	ab_u_y = ab_y*ab_d
	1092	dist_x = ap_x - (ap_xab_u_x+ap_yab_u_y)*ab_u_x
	1093	dist_y = ap_y - (ap_xab_u_x+ap_yab_u_y)*ab_u_y
	1094	dist_point_to_edge = SQRT( dist_x2 + dist_y2 )
	1095	ENDIF
	1096
	1097	RETURN
	1098
	1099	END FUNCTION dist_point_to_edge
	1100
	1101	!
	1102	!-- Returns true if the line segments AB and PQ share an intersection
	1103	FUNCTION intersect ( ax, ay, bx, by, px, py, qx, qy )
	1104
	1105	IMPLICIT NONE
	1106
	1107	LOGICAL :: intersect !< return value; TRUE if intersection was found
	1108	LOGICAL :: la !< T if a is left of PQ
	1109	LOGICAL :: lb !< T if b is left of PQ
	1110	LOGICAL :: lp !< T if p is left of AB
	1111	LOGICAL :: lq !< T if q is left of AB
	1112	LOGICAL :: poss !< flag that indicates if an intersection is still possible
	1113	LOGICAL :: ra !< T if a is right of PQ
	1114	LOGICAL :: rb !< T if b is right of PQ
	1115	LOGICAL :: rp !< T if p is right of AB
	1116	LOGICAL :: rq !< T if q is right of AB
	1117
	1118	REAL(wp) :: ax !< x-coordinate of point A
	1119	REAL(wp) :: ay !< y-coordinate of point A
	1120	REAL(wp) :: bx !< x-coordinate of point B
	1121	REAL(wp) :: by !< y-coordinate of point B
	1122	REAL(wp) :: px !< x-coordinate of point P
	1123	REAL(wp) :: py !< y-coordinate of point P
	1124	REAL(wp) :: qx !< x-coordinate of point Q
	1125	REAL(wp) :: qy !< y-coordinate of point Q
	1126
	1127	intersect = .FALSE.
	1128	poss = .FALSE.
	1129	!
	1130	!-- Intersection is possible only if P and Q are on opposing sides of AB
	1131	lp = is_left(ax,ay,bx,by,px,py)
	1132	rq = is_right(ax,ay,bx,by,qx,qy)
	1133	IF ( lp .AND. rq ) poss = .TRUE.
	1134	IF ( .NOT. poss ) THEN
	1135	lq = is_left(ax,ay,bx,by,qx,qy)
	1136	rp = is_right(ax,ay,bx,by,px,py)
	1137	IF ( lq .AND. rp ) poss = .TRUE.
	1138	ENDIF
	1139	!
	1140	!-- Intersection occurs only if above test (poss) was true AND
	1141	!-- A and B are on opposing sides of PQ
	1142	IF ( poss ) THEN
	1143	la = is_left(px,py,qx,qy,ax,ay)
	1144	rb = is_right(px,py,qx,qy,bx,by)
	1145	IF ( la .AND. rb ) intersect = .TRUE.
	1146	IF ( .NOT. intersect ) THEN
	1147	lb = is_left(px,py,qx,qy,bx,by)
	1148	ra = is_right(px,py,qx,qy,ax,ay)
	1149	IF ( lb .AND. ra ) intersect = .TRUE.
	1150	ENDIF
	1151	ENDIF
	1152
	1153	RETURN
	1154
	1155	END FUNCTION intersect
	1156
	1157	!
	1158	!-- Calculates if point P is left of the infinite
	1159	!-- line that contains A and B (direction: A to B)
	1160	!-- Concept: 2D rotation of two vectors
	1161	FUNCTION is_left ( ax, ay, bx, by, px, py )
	1162
	1163	IMPLICIT NONE
	1164
	1165	LOGICAL :: is_left !< return value; TRUE if P is left of AB
	1166
	1167	REAL(wp) :: ax !< x-coordinate of point A
	1168	REAL(wp) :: ay !< y-coordinate of point A
	1169	REAL(wp) :: bx !< x-coordinate of point B
	1170	REAL(wp) :: by !< y-coordinate of point B
	1171	REAL(wp) :: px !< x-coordinate of point P
	1172	REAL(wp) :: py !< y-coordinate of point P
	1173	!
	1174	!-- 2D-rotation
	1175	is_left = (bx-ax)(py-ay)-(px-ax)(by-ay) > 0
	1176	!
	1177	!-- False if the point is on the line (or very close)
	1178	IF ( (ABS(ax-px) < .001 .AND. ABS(ay-py) < .001) .OR. &
	1179	(ABS(bx-px) < .001 .AND. ABS(by-py) < .001) ) &
	1180	THEN
	1181	is_left = .FALSE.
	1182	ENDIF
	1183
	1184	RETURN
	1185
	1186	END FUNCTION is_left
	1187
	1188	!
	1189	!-- Calculates if point P is right of the infinite
	1190	!-- line that contains A and B (direction: A to B)
	1191	!-- Concept: 2D rotation of two vectors
	1192	FUNCTION is_right ( ax, ay, bx, by, px, py )
	1193
	1194	IMPLICIT NONE
	1195
	1196	LOGICAL :: is_right !< return value; TRUE if P is right of AB
	1197
	1198	REAL(wp), INTENT(IN) :: ax !< x-coordinate of point A
	1199	REAL(wp), INTENT(IN) :: ay !< y-coordinate of point A
	1200	REAL(wp), INTENT(IN) :: bx !< x-coordinate of point B
	1201	REAL(wp), INTENT(IN) :: by !< y-coordinate of point B
	1202	REAL(wp), INTENT(IN) :: px !< x-coordinate of point P
	1203	REAL(wp), INTENT(IN) :: py !< y-coordinate of point P
	1204
	1205	!
	1206	!-- 2D-rotation
	1207	is_right = (bx-ax)(py-ay)-(px-ax)(by-ay) < 0
	1208	!
	1209	!-- False if the point is on the line (or very close)
	1210	IF ( (ABS(ax-px) < .001 .AND. ABS(ay-py) < .001) .OR. &
	1211	(ABS(bx-px) < .001 .AND. ABS(by-py) < .001) ) &
	1212	THEN
	1213	is_right = .FALSE.
	1214	ENDIF
	1215
	1216	RETURN
	1217
	1218	END FUNCTION is_right
	1219
	1220	END MODULE mod_functions
	1221
	1222	MODULE polygon_creation
	1223
	1224	USE kinds
	1225
	1226	USE mod_functions
	1227
	1228	USE variables
	1229
	1230	USE data_input
	1231
	1232	CONTAINS
	1233
	1234	!------------------------------------------------------------------------------!
	1235	! Description:
	1236	! ------------
	1237	!> Initialisation, allocation, and reading of some input
	1238	!------------------------------------------------------------------------------!
	1239	SUBROUTINE init
	1240
	1241	IMPLICIT NONE
	1242
[3400]	1243	CHARACTER(LEN=255) :: dirname !< dummy to read current working directory
	1244	CHARACTER(LEN=255) :: rundir !< base run directory
	1245	CHARACTER(LEN=255) :: input_trunk !< base filename for run including path
	1246	CHARACTER(LEN=80) :: line !< string to identify namelist
	1247	CHARACTER(LEN=80) :: line_dum !< line dummy for error output
[3159]	1248
[3400]	1249	CHARACTER(LEN=2),DIMENSION(1:5) :: run_pars !< parameters from other namelist
[3159]	1250
[3400]	1251	INTEGER(iwp) :: ie !< end index (string manipulation)
	1252	INTEGER(iwp) :: is !< start index (string manipulation)
	1253	INTEGER(iwp) :: line_counter !< line on which reading error occured
	1254
[3159]	1255	LOGICAL :: p3d_flag = .FALSE. !< indicates whether p3d file was found
	1256
[3168]	1257	NAMELIST /prepro_par/ flag_2d, internal_buildings, tolerance_dp
[3159]	1258
[3400]	1259	WRITE(*,'(1X,A)') &
[3159]	1260	"o----------------------------------------------o", &
	1261	"\| o------------------------------------------o \|", &
[3208]	1262	"\| \| __ ____ ____ ____ \| \|", &
	1263	"\| \| / _\ ( _ \(_ _) ___ ( _ \ \| \|", &
	1264	"\| \| / \ ) __/ )( (___) ) __/ \| \|", &
	1265	"\| \| \_/\_/(__) (__) (__) \| \|", &
[3159]	1266	"\| \| \| \|", &
[3208]	1267	"\| \| Agent Preprocessing Tool for PALM \| \|", &
[3159]	1268	"\| o------------------------------------------o \|", &
	1269	"o----------------------------------------------o"
[3208]	1270
[3159]	1271	!
	1272	!-- Identify run name and Input files
[3400]	1273	CALL GET_ENVIRONMENT_VARIABLE('PWD', dirname)
[3159]	1274	ie = INDEX(dirname, '/', BACK=.TRUE.)
	1275	is = INDEX(dirname(1:ie-1), '/', BACK=.TRUE.)
[3208]	1276	IF ( TRIM(ADJUSTL(dirname(ie+1:))) /= 'INPUT' ) THEN
	1277	WRITE(*,'(3X,A)') 'NavMesh was called from', &
	1278	' ', TRIM(ADJUSTL(dirname)), ' ', &
	1279	'and is now aborting. Please call this tool', &
	1280	'from the INPUT-folder of your job directory.'
	1281	STOP
	1282	ENDIF
[3159]	1283	runname = TRIM(ADJUSTL(dirname(is+1:ie-1)))
	1284	rundir = TRIM(ADJUSTL(dirname(1:ie)))
	1285	input_trunk = TRIM(rundir)//'INPUT/'//TRIM(runname)
[3208]	1286
[3159]	1287	!
	1288	!-- Check for parameter file
	1289	INQUIRE( FILE = TRIM( input_trunk )//'_p3d', EXIST = p3d_flag )
	1290	IF ( .NOT. p3d_flag ) THEN
	1291	WRITE(*,'(3(3X,A,/))') 'No _p3d file was found. Aborting.', &
	1292	'I was looking for the file', &
	1293	TRIM( input_trunk )//'_p3d'
	1294	STOP
	1295	ELSE
	1296	WRITE(*,'(2(3X,A,/))') 'The following input file will be used:', &
	1297	TRIM( input_trunk )//'_p3d'
	1298	ENDIF
	1299
	1300	!
	1301	!-- Read run parameters from run parameter file (_p3d), though not from
	1302	!-- namelist.
	1303	run_pars = (/'dx','dy','dz','nx','ny'/)
	1304	OPEN ( 11, FILE=TRIM(input_trunk)//'_p3d', FORM='FORMATTED', &
	1305	STATUS='OLD' )
	1306	DO i = 1, SIZE(run_pars)
	1307	REWIND ( 11 )
	1308	line = ' '
	1309	DO WHILE ( INDEX( line, run_pars(i) ) == 0 )
	1310	READ ( 11, '(A)', END=10 ) line
	1311	IF ( INDEX(line, '!') /= 0 ) THEN
	1312	IF ( INDEX(line, run_pars(i)) > INDEX(line, '!' ) ) THEN
	1313	line = ' '
	1314	CYCLE
	1315	ENDIF
	1316	ENDIF
	1317	ENDDO
	1318	line = TRIM(ADJUSTL(line(INDEX(line,'=')+1:INDEX(line,',')-1)))
	1319	SELECT CASE (i)
	1320	CASE(1)
	1321	READ(line,*) dx
	1322	CASE(2)
	1323	READ(line,*) dy
	1324	CASE(3)
	1325	READ(line,*) dz
	1326	CASE(4)
	1327	READ(line,*) nx
	1328	CASE(5)
	1329	READ(line,*) ny
	1330	CASE DEFAULT
	1331	END SELECT
	1332	ENDDO
	1333	10 CONTINUE
	1334
	1335	!
	1336	!-- Try to find prepro package
	1337	REWIND ( 11 )
	1338	line = ' '
	1339	DO WHILE ( INDEX( line, '&prepro_par' ) == 0 )
[3400]	1340	READ ( 11, '(A)', END=40 ) line
[3159]	1341	ENDDO
	1342	BACKSPACE ( 11 )
	1343
	1344	!
	1345	!-- Read user-defined namelist
[3400]	1346	READ ( 11, prepro_par, ERR = 20, END = 40 )
	1347	GOTO 40
[3159]	1348
[3400]	1349	20 BACKSPACE( 11 )
	1350	READ( 11 , '(A)') line
	1351
	1352	line_dum = ' '
	1353	line_counter = 0
	1354
	1355	REWIND( 11 )
	1356	DO WHILE ( INDEX( line_dum, TRIM(line) ) == 0 )
	1357	READ ( 11, '(A)', END=30 ) line_dum
	1358	line_counter = line_counter + 1
	1359	ENDDO
	1360
	1361	30 WRITE( *, '(A,/,A,I3,A,/,A)' ) 'Error(s) in NAMELIST prepro_par.', &
	1362	'Reading fails on line ', line_counter, &
	1363	' at ', TRIM(ADJUSTL(line))
	1364	STOP
	1365
	1366	40 CONTINUE
[3159]	1367	CLOSE( 11 )
	1368
	1369	!
	1370	!-- If tolerance_dp was not set, put in default values
[3198]	1371	DO i = 0, 2
[3259]	1372	IF ( tolerance_dp(i) > 999998.0_wp ) THEN
[3198]	1373	tolerance_dp(i) = SQRT(dxdy)1.41/(2**i)
	1374	ELSE
	1375	tolerance_dp(i) = tolerance_dp(i)SQRT(dxdy)
	1376	ENDIF
	1377	ENDDO
	1378
[3159]	1379	!
	1380	!-- Allocate arrays
	1381	ALLOCATE(obstacle_height(-3:nx+3,-3:ny+3), polygon_id(-3:nx+3,-3:ny+3), &
	1382	wall_flags_0(-3:nx+3,-3:ny+3), grid(-3:nx+3,-3:ny+3))
	1383	!
	1384	!-- Set null_vertex
[3259]	1385	CALL set_vertex(null_vertex,0.0_wp,0.0_wp)
[3159]	1386	!
	1387	!-- Some initializations
	1388	ddx = 1./dx
	1389	ddy = 1./dy
	1390
	1391	polygon_id = 0
	1392	obstacle_height = 0.
	1393
	1394	grid%checked = .FALSE.
	1395	grid(-3:-1,:)%checked = .TRUE.
	1396	grid(nx+1:nx+3,:)%checked = .TRUE.
	1397	grid(:,-3:-1)%checked = .TRUE.
	1398	grid(:,ny+1:ny+3)%checked = .TRUE.
	1399	grid%polygon_id = 0
	1400	!
	1401	!-- Open files and topography/building data
	1402	CALL netcdf_data_input_topo ( TRIM(input_trunk) )
	1403
	1404	END SUBROUTINE init
	1405
	1406	!------------------------------------------------------------------------------!
	1407	! Description:
	1408	! ------------
	1409	!> Identifies all grid boxes that belong to a building and assigns a building
	1410	!> number to each grid box.
	1411	!> Method: Scans each grid point. If a grid point was not previously checked
	1412	!> and contains a building, it is added to a new polygon and marked as
	1413	!> checked. Then, all its neighbors that also contain a building are
	1414	!> added to the same polygon and marked as checked. All neighbors of
	1415	!> neighbors are subsequently found, added and checked until none are
	1416	!> left. Then, the initial scan continues, skipping already checked
	1417	!> grid points. Once a grid point with a new building is found, the
	1418	!> polygon_id increases and the grid point and all others that belong
	1419	!> to the same building are added as described above.
	1420	!> NOTE: This procedure will identify grid points that are only connected
	1421	!> diagonally (share only one point with each other) as connected and
	1422	!> have them belonging to the same building. This is necessary, as an
	1423	!> agent will not be able to traverse between these grid points and the
	1424	!> navigation mesh will therefore have to make him circumvent this
	1425	!> point.
	1426	!------------------------------------------------------------------------------!
	1427	SUBROUTINE identify_polygons
	1428
	1429	IMPLICIT NONE
	1430
	1431	INTEGER(iwp) :: ii !< local counter
	1432	INTEGER(iwp) :: il !< local counter
	1433	INTEGER(iwp) :: jj !< local counter
	1434	INTEGER(iwp) :: jl !< local counter
	1435	INTEGER(iwp) :: gpil !< number of grid points in list
	1436	INTEGER(iwp) :: gpta !< number of grid points to add to grid_list
	1437
	1438	TYPE(grid_point), DIMENSION(1:7) :: add_to_grid_list !< grid points to be added to the list
	1439
	1440	TYPE(grid_point), DIMENSION(:), ALLOCATABLE :: dummy_grid_list !< dummy for reallocation of grid_list
	1441	TYPE(grid_point), DIMENSION(:), ALLOCATABLE :: grid_list !< list of grid points that belong to the current building but whose neighbors have not been checked yet
	1442
	1443	!
	1444	!-- Initialize wall_flags array: 1 where no buildings, 0 where buildings
	1445	wall_flags_0 = 1
	1446	DO i = 0, nx
	1447	DO j = 0, ny
	1448	IF ( obstacle_height(i,j) > 0 ) THEN
	1449	wall_flags_0(i,j) = 0
	1450	ENDIF
	1451	ENDDO
	1452	ENDDO
	1453	DEALLOCATE(obstacle_height)
	1454	polygon_counter = 0
	1455	gpil = 0
	1456	gpta = 0
	1457	ALLOCATE(grid_list(1:100))
	1458	!
	1459	!-- Declare all grid points that contain no buildings as already checked.
	1460	!-- This way, these points will be skipped in the following calculation and
	1461	!-- will have polygon_id = 0
	1462	DO i = 0, nx
	1463	DO j = 0, ny
	1464	IF ( BTEST( wall_flags_0(i,j), 0 ) ) THEN
	1465	grid(i,j)%checked = .TRUE.
	1466	ENDIF
	1467	ENDDO
	1468	ENDDO
	1469	!
	1470	!-- Check all grid points and process them
	1471	DO i = 0, nx
	1472	DO j = 0, ny
	1473	!
	1474	!-- If the current grid point has not been checked, mark it as checked.
	1475	!-- As it is the first point belonging to a new building, increase the
	1476	!-- polygon_id counter and associate the grid point with that id.
	1477	IF ( .NOT. grid(i,j)%checked ) THEN
	1478	polygon_counter = polygon_counter + 1
	1479	grid(i,j)%polygon_id = polygon_counter
	1480	grid(i,j)%checked = .TRUE.
	1481	!
	1482	!-- Check if any neighbors of the found grid point are part of a
	1483	!-- building too. If so, add them to the list of grid points
	1484	!-- that have to be checked and associate them with the same polygon
	1485	gpta = 0
	1486	DO ii = i-1, i+1
	1487	DO jj = j-1, j+1
	1488	IF ( ii == i .AND. jj == j ) CYCLE
	1489	IF ( .NOT. grid(ii,jj)%checked ) THEN
	1490	gpta = gpta + 1
	1491	add_to_grid_list(gpta)%i = ii
	1492	add_to_grid_list(gpta)%j = jj
	1493	ENDIF
	1494	ENDDO
	1495	ENDDO
	1496
	1497	!
	1498	!-- Change size of grid_list if it becomes too small
	1499	IF ( gpil + gpta > SIZE(grid_list) ) THEN
	1500	ALLOCATE(dummy_grid_list(1:gpil))
	1501	dummy_grid_list = grid_list(1:gpil)
	1502	DEALLOCATE(grid_list)
	1503	ALLOCATE(grid_list(1:2*(gpil+gpta)))
	1504	grid_list(1:gpil) = dummy_grid_list(1:gpil)
	1505	DEALLOCATE(dummy_grid_list)
	1506	ENDIF
	1507	!
	1508	!-- If there are grid points to add to grid_list, add them
	1509	IF ( gpta > 0 ) THEN
	1510	grid_list(gpil+1:gpil+gpta) = add_to_grid_list(1:gpta)
	1511	gpil = gpil + gpta
	1512	ENDIF
	1513	!
	1514	!-- Handle all grid points in grid_list until there are none left
	1515	DO WHILE (gpil>0)
	1516	il = grid_list(gpil)%i
	1517	jl = grid_list(gpil)%j
	1518	grid(il,jl)%polygon_id = polygon_counter
	1519	grid(il,jl)%checked = .TRUE.
	1520	!
	1521	!-- this grid point in the list is processed, so the number of
	1522	!-- grid points in the list can be reduced by one
	1523	gpil = gpil - 1
	1524	gpta = 0
	1525	!
	1526	!-- For the current grid point, check if any unchecked
	1527	!-- neighboring grid points also contain a building. All such
	1528	!-- grid points are added to the list of grid points to be
	1529	!-- handled in this loop
	1530	DO ii = il-1, il+1
	1531	DO jj = jl-1, jl+1
	1532	IF ( jj == jl .AND. ii == il ) CYCLE
	1533	IF ( .NOT. grid(ii,jj)%checked ) &
	1534	THEN
	1535	gpta = gpta + 1
	1536	add_to_grid_list(gpta)%i = ii
	1537	add_to_grid_list(gpta)%j = jj
	1538	ENDIF
	1539	ENDDO
	1540	ENDDO
	1541	!
	1542	!-- Change size of grid list if it becomes too small
	1543	IF ( gpil + gpta > SIZE(grid_list) ) THEN
	1544	ALLOCATE(dummy_grid_list(1:gpil))
	1545	dummy_grid_list = grid_list(1:gpil)
	1546	DEALLOCATE(grid_list)
	1547	ALLOCATE(grid_list(1:2*(gpil+gpta)))
	1548	grid_list(1:gpil) = dummy_grid_list(1:gpil)
	1549	DEALLOCATE(dummy_grid_list)
	1550	ENDIF
	1551	!
	1552	!-- If there are grid points to add to list, add them
	1553	IF ( gpta > 0 ) THEN
	1554	grid_list(gpil+1:gpil+gpta) = add_to_grid_list(1:gpta)
	1555	gpil = gpil + gpta
	1556	ENDIF
	1557	ENDDO
	1558	ENDIF
	1559	ENDDO
	1560	ENDDO
	1561	DEALLOCATE(grid_list)
	1562	!
	1563	!-- Set size of polygon array and initialize
	1564	ALLOCATE(polygons(1:polygon_counter))
	1565	polygons%nov = 0
	1566
	1567	END SUBROUTINE identify_polygons
	1568
	1569	!------------------------------------------------------------------------------!
	1570	! Description:
	1571	! ------------
	1572	!> Identifies the corners of the PALM building topography and adds them to
	1573	!> a specific polygon for each building as vertices. This converts the gridded
	1574	!> building data into one polygon per building that contains the coordinates of
	1575	!> each inner and outer corner of that building as vertices.
	1576	!> A grid point contains an outer corner if it's part of a building and exactly
	1577	!> one of its horizontal and one of its vertical neighbors is also part of a
	1578	!> building (4 cases).
	1579	!> A grid point contains an inner corner if it's not part of a building and
	1580	!> exactly one of its horizontal, one of its diagonal and one of its vertical
	1581	!> neighbors are each part of a building and in turn neighbors
	1582	!> to each other (4 cases).
	1583	!------------------------------------------------------------------------------!
	1584	SUBROUTINE identify_corners
	1585
	1586	IMPLICIT NONE
	1587
	1588	INTEGER(iwp) :: p_id !< current polygon_id
	1589	!
	1590	!-- For all grid points, check whether it contains one or more corners
	1591	DO i = 0, nx
	1592	DO j = 0, ny
	1593	!
	1594	!-- First, check if grid contains topography and has a corner.
	1595	IF ( .NOT. BTEST( wall_flags_0(i,j), 0 ) ) THEN
	1596	!
	1597	!-- Corner at south left edge of grid cell
	1598	IF ( BTEST( wall_flags_0(i-1,j), 0 ) .AND. &
	1599	BTEST( wall_flags_0(i,j-1), 0 )) &
	1600	THEN
	1601	p_id = grid(i,j)%polygon_id
	1602	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1603	nov = polygons(p_id)%nov
[3259]	1604	CALL set_vertex(dummy_vertex, idx, jdy)
[3159]	1605	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1606	ENDIF
	1607	!
	1608	!-- Corner at north left edge of grid cell
	1609	IF ( BTEST( wall_flags_0(i-1,j), 0 ) .AND. &
	1610	BTEST( wall_flags_0(i,j+1), 0 )) &
	1611	THEN
	1612	p_id = grid(i,j)%polygon_id
	1613	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1614	nov = polygons(p_id)%nov
[3259]	1615	CALL set_vertex(dummy_vertex, idx, (j+1)dy)
[3159]	1616	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1617	ENDIF
	1618	!
	1619	!-- Corner at north right edge of grid cell
	1620	IF ( BTEST( wall_flags_0(i+1,j), 0 ) .AND. &
	1621	BTEST( wall_flags_0(i,j+1), 0 )) &
	1622	THEN
	1623	p_id = grid(i,j)%polygon_id
	1624	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1625	nov = polygons(p_id)%nov
[3259]	1626	CALL set_vertex(dummy_vertex, (i+1)dx, (j+1)dy)
[3159]	1627	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1628	ENDIF
	1629	!
	1630	!-- Corner at south right edge of grid cell
	1631	IF ( BTEST( wall_flags_0(i+1,j), 0 ) .AND. &
	1632	BTEST( wall_flags_0(i,j-1), 0 )) &
	1633	THEN
	1634	p_id = grid(i,j)%polygon_id
	1635	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1636	nov = polygons(p_id)%nov
[3259]	1637	CALL set_vertex(dummy_vertex, (i+1)dx, jdy)
[3159]	1638	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1639	ENDIF
	1640	!
	1641	!-- Second, check if grid contains no topography and has a corner.
	1642	ELSE
	1643	!
	1644	!-- Corner at south left edge of grid cell
	1645	IF ( .NOT. BTEST( wall_flags_0(i-1,j), 0 ) .AND. &
	1646	.NOT. BTEST( wall_flags_0(i,j-1), 0 ) .AND. &
	1647	.NOT. BTEST( wall_flags_0(i-1,j-1), 0 ) ) &
	1648	THEN
	1649	p_id = grid(i-1,j-1)%polygon_id
	1650	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1651	nov = polygons(p_id)%nov
[3259]	1652	CALL set_vertex(dummy_vertex, idx, jdy)
[3159]	1653	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1654	ENDIF
	1655	!
	1656	!-- Corner at north left edge of grid cell
	1657	IF ( .NOT. BTEST( wall_flags_0(i-1,j), 0 ) .AND. &
	1658	.NOT. BTEST( wall_flags_0(i,j+1), 0 ) .AND. &
	1659	.NOT. BTEST( wall_flags_0(i-1,j+1), 0 ) ) &
	1660	THEN
	1661	p_id = grid(i-1,j+1)%polygon_id
	1662	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1663	nov = polygons(p_id)%nov
[3259]	1664	CALL set_vertex(dummy_vertex, idx, (j+1)dy)
[3159]	1665	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1666	ENDIF
	1667	!
	1668	!-- Corner at north right edge of grid cell
	1669	IF ( .NOT. BTEST( wall_flags_0(i+1,j), 0 ) .AND. &
	1670	.NOT. BTEST( wall_flags_0(i,j+1), 0 ) .AND. &
	1671	.NOT. BTEST( wall_flags_0(i+1,j+1), 0 ) ) &
	1672	THEN
	1673	p_id = grid(i+1,j+1)%polygon_id
	1674	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1675	nov = polygons(p_id)%nov
[3259]	1676	CALL set_vertex(dummy_vertex, (i+1)dx, (j+1)dy)
[3159]	1677	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1678	ENDIF
	1679	!
	1680	!-- Corner at south right edge of grid cell
	1681	IF ( .NOT. BTEST( wall_flags_0(i+1,j), 0 ) .AND. &
	1682	.NOT. BTEST( wall_flags_0(i,j-1), 0 ) .AND. &
	1683	.NOT. BTEST( wall_flags_0(i+1,j-1), 0 ) ) &
	1684	THEN
	1685	p_id = grid(i+1,j-1)%polygon_id
	1686	polygons(p_id)%nov = polygons(p_id)%nov + 1
	1687	nov = polygons(p_id)%nov
[3259]	1688	CALL set_vertex(dummy_vertex, (i+1)dx, jdy)
[3159]	1689	CALL add_vertex_to_polygon(dummy_vertex, p_id, nov)
	1690	ENDIF
	1691	ENDIF
	1692	ENDDO
	1693	ENDDO
	1694
	1695	END SUBROUTINE identify_corners
	1696
	1697	!------------------------------------------------------------------------------!
	1698	! Description:
	1699	! ------------
	1700	!> Initializes a vertex
	1701	!------------------------------------------------------------------------------!
[3259]	1702	SUBROUTINE set_vertex (in_vertex, x, y)
[3159]	1703
	1704	IMPLICIT NONE
	1705
	1706	REAL(wp) :: x !< x-coordinate of vertex position
	1707	REAL(wp) :: y !< y-coordinate of vertex position
	1708
	1709	TYPE(vertex_type) :: in_vertex !< vertex to be set
	1710
	1711	in_vertex%delete = .FALSE.
	1712	in_vertex%x = x
	1713	in_vertex%y = y
	1714
	1715	END SUBROUTINE set_vertex
	1716
	1717	!------------------------------------------------------------------------------!
	1718	! Description:
	1719	! ------------
	1720	!> Adds an existing vertex to the polygon with ID p_id at position in_nov
	1721	!------------------------------------------------------------------------------!
	1722	SUBROUTINE add_vertex_to_polygon ( in_vertex, p_id, in_nov)
	1723
	1724	IMPLICIT NONE
	1725
	1726	INTEGER(iwp) :: in_nov !< counter of vertex being added to polygon
	1727	INTEGER(iwp) :: p_id !< polygon ID
	1728	INTEGER(iwp) :: sop !< size of vertices array
	1729
	1730	TYPE(vertex_type) :: in_vertex !< vertex to be added
	1731
	1732	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: dummy_v_list !< for reallocation
	1733
	1734	polygon => polygons(p_id)
	1735	!
	1736	!-- Allocate and initialize the vertex array of the polygon, if necessary
	1737	IF ( .NOT. ALLOCATED(polygon%vertices) ) THEN
	1738	ALLOCATE(polygon%vertices(1:100))
	1739	polygon%vertices = null_vertex
	1740	ENDIF
	1741	!
	1742	!-- Adjust size of polygon, if necessary
	1743	sop = SIZE(polygon%vertices)
	1744	IF ( in_nov > sop ) THEN
	1745	ALLOCATE(dummy_v_list(1:sop))
	1746	dummy_v_list(1:sop) = polygon%vertices(1:sop)
	1747	DEALLOCATE(polygon%vertices)
	1748	ALLOCATE(polygon%vertices(1:2*sop))
	1749	polygon%vertices = null_vertex
	1750	polygon%vertices(1:sop) = dummy_v_list(1:sop)
	1751	DEALLOCATE(dummy_v_list)
	1752	ENDIF
	1753	polygon%vertices(in_nov) = in_vertex
	1754	END SUBROUTINE add_vertex_to_polygon
	1755
	1756	!------------------------------------------------------------------------------!
	1757	! Description:
	1758	! ------------
	1759	!> Sorts the vertices of a polygon in a counter-clockwise fashion. During this
	1760	!> process, all vertices that are not part of the hull of the building
	1761	!> (inner courtyards) are deleted.
	1762	!------------------------------------------------------------------------------!
	1763	SUBROUTINE sort_polygon(i_p)
	1764
	1765	IMPLICIT NONE
	1766
	1767	LOGICAL :: starting_vertex_found
	1768
	1769	INTEGER(iwp) :: counter !< counter for potential starting vertices
	1770	INTEGER(iwp) :: id_neighbor !< final ID of neighboring vertex
	1771	INTEGER(iwp) :: id_neighbor1 !< ID of first potential neighbor
	1772	INTEGER(iwp) :: id_neighbor2 !< ID of second potential neighbor
	1773	INTEGER(iwp) :: il !< local counter
	1774	INTEGER(iwp) :: i_p !< index of the current polygon
	1775	INTEGER(iwp) :: noc !< number of candidates
	1776	INTEGER(iwp) :: nosv !< number of sorted vertices
	1777	INTEGER(iwp) :: xe !< x-end-index for building search
	1778	INTEGER(iwp) :: xs !< x-start-index for building search
	1779	INTEGER(iwp) :: ye !< y-end-index for building search
	1780	INTEGER(iwp) :: ys !< y-start-index for building search
	1781
	1782	INTEGER, DIMENSION(:), ALLOCATABLE :: candidate_id !< ID of the potential neighbors stored in 'candidates'
	1783	INTEGER, DIMENSION(:), ALLOCATABLE :: dummy_id_arr !< used for resizing
	1784
	1785	REAL(wp) :: dist !< distance of one vertex to its neighbor
	1786	REAL(wp) :: m_x !< min/max x-value of polygon used for starting vertex
	1787	REAL(wp) :: m_y !< min/max y-value of polygon used for starting vertex
	1788
	1789	TYPE(vertex_type) :: current_v !< current vertex
	1790	TYPE(vertex_type) :: dummy_vertex !< dummy vertex for reordering
	1791
	1792	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: candidates !< potential neighbors of the current vertex
	1793	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: dummy_vertex_arr !< used for resizing
	1794	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: sorted_p !< vertices that have been sorted
	1795
	1796	starting_vertex_found = .FALSE.
	1797	ALLOCATE(sorted_p(1:nov))
	1798	sorted_p(1:nov) = polygon%vertices(1:nov)
	1799	!
	1800	!-- Identify a vertex that is certainly a part of the outer hull of the
	1801	!-- current polygon: Get rightmost border of polygon (or if that
	1802	!-- coincides with model border, leftmost) border. Then of those points,
	1803	!-- get northmost (or if that coincides with model domain border, southmost).
	1804	!-- This identifies exactly one point that is then set to the first index.
	1805	counter = 0
	1806	IF ( MAXVAL(sorted_p%x) < nx*dx ) THEN
	1807	m_x = MAXVAL(sorted_p%x)
	1808	ELSE
	1809	m_x = MINVAL(sorted_p%x)
	1810	ENDIF
	1811	DO il = 1, nov
[3259]	1812	IF ( ABS( sorted_p(il)%x - m_x ) < .01 * dx ) THEN
[3159]	1813	counter = counter + 1
	1814	dummy_vertex = sorted_p(il)
	1815	sorted_p(il) = sorted_p(counter)
	1816	sorted_p(counter) = dummy_vertex
	1817	ENDIF
	1818	ENDDO
	1819	IF ( MAXVAL(sorted_p(1:counter)%y) < ny*dy ) THEN
	1820	m_y = MAXVAL(sorted_p(1:counter)%y)
	1821	ELSE
	1822	m_y = MINVAL(sorted_p(1:counter)%y)
	1823	ENDIF
	1824	DO il = 1, counter
[3259]	1825	IF ( ABS(sorted_p(il)%y - m_y) < .01 * dy ) THEN
[3159]	1826	dummy_vertex = sorted_p(il)
	1827	sorted_p(il) = sorted_p(1)
	1828	sorted_p(1) = dummy_vertex
	1829	starting_vertex_found = .TRUE.
	1830	EXIT
	1831	ENDIF
	1832	ENDDO
	1833	!
	1834	!-- If no starting vertex was found for the current polygon, it will be
	1835	!-- deleted and an error message thrown
	1836	IF ( .NOT. starting_vertex_found ) THEN
[3400]	1837	WRITE(*,'(A,/,A,1X,I6,/,A)') &
[3159]	1838	'An error occured during polygon sorting:', &
	1839	'no starting vertex could be found for polygon', &
	1840	i_p, 'This polygon contains the following vertices (x/y)'
	1841	DO il = 1, nov
[3400]	1842	WRITE(*,'(4X,F8.1,1X,F8.1)') &
[3159]	1843	polygon%vertices(il)%x, polygon%vertices(il)%x
	1844	ENDDO
	1845	WRITE(*,'(A,/,A)') &
	1846	'This polygon will be skipped during sorting and deleted',&
	1847	'For details on the procedure, see SUBROUTINE sort_polygon.'
	1848	polygon%vertices%delete = .TRUE.
	1849	polygons(i_p)%nov = 0
	1850	CALL delete_empty_polygons
	1851	!
	1852	!-- Find the unique neighbor of the current vertex. For this, first
	1853	!-- determine all possible candidates. Of those, keep only the ones that
	1854	!-- are connected to the current vertex along a building edge (the polygon
	1855	!-- is sorted counter-clockwise. Therefore, the building is always on the
	1856	!-- left-hand side of the connecting line from the current vertex to its
	1857	!-- potential neighbor). This leaves a maximum of two possible neighbors.
	1858	!-- This is only the case if the current vertex is the point that diagonally
	1859	!-- connects two parts of the same building. In that case, the vertex that
	1860	!-- lies to the right of the connecting line between the previous and
	1861	!-- current vertex is the neighbor.
	1862	ELSE
	1863	DO nosv = 1, nov
	1864	current_v = sorted_p(nosv)
	1865	noc = 0
	1866	ALLOCATE(candidates(1:100), candidate_id(1:100))
	1867	!
	1868	!-- Find all candidates that could be neighbors of current vertex:
	1869	!-- these are those vertices that share the same x- or y-coordinate
	1870	!-- with the current vertex, as the vertices are all inner and outer
	1871	!-- corners of the gridded building data
	1872	IF ( nosv < nov ) THEN
	1873	DO il = nosv+1, nov
[3259]	1874	IF ( ABS( current_v%x - sorted_p(il)%x ) < .01 * dx .OR. &
	1875	ABS( current_v%y - sorted_p(il)%y ) < .01 * dy ) &
[3159]	1876	THEN
	1877	!
	1878	!-- If necessary, resize arrays for candidates
	1879	IF ( noc >= SIZE(candidates) ) THEN
	1880	ALLOCATE(dummy_vertex_arr(1:noc), dummy_id_arr(1:noc))
	1881	dummy_vertex_arr(1:noc) = candidates(1:noc)
	1882	dummy_id_arr(1:noc) = candidate_id(1:noc)
	1883	DEALLOCATE(candidates, candidate_id)
	1884	ALLOCATE(candidates(1:2noc), candidate_id(1:2noc))
	1885	candidates(1:noc) = dummy_vertex_arr(1:noc)
	1886	candidate_id(1:noc) = dummy_id_arr(1:noc)
	1887	DEALLOCATE(dummy_vertex_arr, dummy_id_arr)
	1888	ENDIF
	1889	noc = noc +1
	1890	candidates(noc) = sorted_p(il)
	1891	candidate_id(noc) = il
	1892	ENDIF
	1893	ENDDO
	1894	ENDIF
	1895	!
	1896	!-- Check which one of the candidates is the neighbor of the current
	1897	!-- vertex. This is done by several tests that would exclude the
	1898	!-- candidate from being the neighbor. Each successful test will
	1899	!-- therefore result in a cycle to the next candidate. Only if all
	1900	!-- all tests fail, is the candidate one of a maximum of two possible
	1901	!-- neighbors.
[3259]	1902	id_neighbor = -999
[3159]	1903	id_neighbor1 = -999
	1904	id_neighbor2 = -999
	1905	DO il = 1, noc
	1906	!
	1907	!-- Exclude the possibility of a vertex with the same coordinates
	1908	!-- being chosen as the neighbor. (dist < .9*dx)
	1909	!-- NOTE: this could happen, if part of a building is only connected
	1910	!-- to the rest of the building diagonally. In that case, the
	1911	!-- same point is added to the polygon twice. This is necessary
	1912	!-- and not redundant! Two such points can never be neighbors.
	1913	!-- Example: the north right corner of grid point i,j
	1914	!-- AND the south left corner of grid point i+1,j+1.
	1915	!-- See SUBROUTINE identify_corners for the identification
	1916	!-- method.
	1917	!-- Also, exclude a connection back to the coordinates of the
	1918	!-- previous vertex.
	1919	dist = SQRT( (candidates(il)%x - current_v%x)**2 + &
	1920	(candidates(il)%y - current_v%y)**2 )
	1921	IF ( nosv > 1 ) THEN
[3259]	1922	IF ( dist < .9 * dx .OR. &
	1923	( ( ABS( sorted_p(nosv-1)%x &
	1924	- candidates(il)%x ) < .01 * dx ) .AND. &
	1925	( ABS( sorted_p(nosv-1)%y &
	1926	- candidates(il)%y ) < .01 * dy ) ) ) &
[3159]	1927	THEN
	1928	CYCLE
	1929	ENDIF
	1930	ENDIF
	1931	!
	1932	!-- Check if there is a building all along only the left-hand side
	1933	!-- of the connecting line from current vertex to potential neighbor
	1934	!-- (4 cases)
	1935	!-- First: for vertical connection
[3259]	1936	IF ( ABS( candidates(il)%x - current_v%x ) < .01 * dx ) THEN
[3159]	1937	xs = NINT(current_v%x*ddx)-1
	1938	xe = xs + 1
	1939	!
	1940	!-- Case 1: ys < ye, edge from south to north, building must be
	1941	!-- exclusively in all grid cells left of the edge
	1942	IF ( current_v%y < candidates(il)%y ) THEN
	1943	ys = NINT(current_v%y*ddy)
	1944	ye = NINT(candidates(il)%y*ddy)-1
	1945	IF ( .NOT.( ALL( .NOT. BTEST( wall_flags_0(xs,ys:ye), 0))&
	1946	.AND.( ALL( BTEST( wall_flags_0(xe,ys:ye), 0 ) ) )))&
	1947	THEN
	1948	CYCLE
	1949	ENDIF
	1950	!
	1951	!-- Case 2: ys > ye, edge from north to south, building must be
	1952	!-- exclusively in all grid cells right of the edge
	1953	ELSEIF ( current_v%y > candidates(il)%y ) THEN
	1954	ys = NINT(current_v%y*ddy)-1
	1955	ye = NINT(candidates(il)%y*ddy)
	1956	IF ( .NOT.( ALL( .NOT. BTEST( wall_flags_0(xe,ye:ys), 0))&
	1957	.AND.( ALL( BTEST( wall_flags_0(xs,ye:ys), 0 ) ) )))&
	1958	THEN
	1959	CYCLE
	1960	ENDIF
	1961	ENDIF
	1962	!
	1963	!-- Horizontal connection
[3259]	1964	ELSEIF ( ABS( candidates(il)%y - current_v%y ) < .01 * dy ) THEN
[3159]	1965
	1966	ys = NINT(current_v%y*ddy)-1
	1967	ye = ys + 1
	1968	!
	1969	!-- Case 3: xs > xe, edge from right to left, building must be
	1970	!-- exclusively in all grid cells south of the edge
	1971	IF ( current_v%x > candidates(il)%x ) THEN
	1972	xs = NINT(current_v%x*ddx)-1
	1973	xe = NINT(candidates(il)%x*ddx)
	1974	IF ( .NOT.( ALL( .NOT. BTEST( wall_flags_0(xe:xs,ys), 0))&
	1975	.AND.( ALL( BTEST( wall_flags_0(xe:xs,ye), 0 ) ) )))&
	1976	THEN
	1977	CYCLE
	1978	ENDIF
	1979	!
	1980	!-- Case 4: xs < xe, edge from left to right, building must be
	1981	!-- exclusively in all grid cells north of the edge
	1982	ELSEIF ( current_v%x < candidates(il)%x ) THEN
	1983	xs = NINT(current_v%x*ddx)
	1984	xe = NINT(candidates(il)%x*ddx)-1
	1985	IF ( .NOT.( ALL( .NOT. BTEST( wall_flags_0(xs:xe,ye), 0))&
	1986	.AND.( ALL( BTEST( wall_flags_0(xs:xe,ys), 0 ) ) )))&
	1987	THEN
	1988	CYCLE
	1989	ENDIF
	1990	ENDIF
	1991	ENDIF
	1992	!
	1993	!-- After the tests, only two potential neighbors are possible. The
	1994	!-- one found first will get id_neighbor1, the possible 2nd one will
	1995	!-- get id_neighbor2
	1996	IF ( id_neighbor1 == -999 ) THEN
	1997	id_neighbor1 = candidate_id(il)
	1998	ELSEIF ( id_neighbor1 /= -999 .AND. &
[3259]	1999	( ( ABS( sorted_p(id_neighbor1)%x - candidates(il)%x ) &
	2000	> .01 * dx ) .OR. &
	2001	( ABS( sorted_p(id_neighbor1)%y - candidates(il)%y ) &
	2002	> .01 * dy ) ) ) &
[3159]	2003	THEN
	2004	id_neighbor2 = candidate_id(il)
	2005	ENDIF
	2006	ENDDO
	2007	!
	2008	!-- If two potential neighbors were found, determine the one that is on
	2009	!-- the right hand side of the line connecting the current and previous
	2010	!-- vertex. It is the real neighbor.
	2011	IF ( id_neighbor2 /= -999 .AND. nosv > 1 ) THEN
	2012	IF ( is_right(sorted_p(nosv-1)%x,sorted_p(nosv-1)%y, &
	2013	current_v%x,current_v%y, &
	2014	sorted_p(id_neighbor1)%x,sorted_p(id_neighbor1)%y) )&
	2015	THEN
	2016	id_neighbor = id_neighbor1
	2017	ELSEIF ( is_right(sorted_p(nosv-1)%x,sorted_p(nosv-1)%y, &
	2018	current_v%x,current_v%y, &
	2019	sorted_p(id_neighbor2)%x,sorted_p(id_neighbor2)%y) )&
	2020	THEN
	2021	id_neighbor = id_neighbor2
	2022	ENDIF
	2023	ELSE
	2024	id_neighbor = id_neighbor1
	2025	ENDIF
	2026	!
	2027	!-- Put the found neighbor at next index in sorted array and move the
	2028	!-- unsorted vertices back one index. This way, only yet unsorted
	2029	!-- vertices are eligible to be candidates during the next iteration.
	2030	IF (id_neighbor /= nosv + 1 .AND. id_neighbor /= -999) THEN
	2031	dummy_vertex = sorted_p(id_neighbor)
	2032	sorted_p(nosv+2:id_neighbor) = sorted_p(nosv+1:id_neighbor-1)
	2033	sorted_p(nosv+1) = dummy_vertex
	2034	!
	2035	!-- If no neighbor was found, sorting is done for this polygon
	2036	ELSEIF ( id_neighbor == -999 ) THEN
	2037	DEALLOCATE(candidates,candidate_id)
	2038	EXIT
	2039	ENDIF
	2040	DEALLOCATE(candidates,candidate_id)
	2041	ENDDO
	2042	!
	2043	!-- Sorting is done. Reduce size (which means get rid of vertices
	2044	!-- that are not part of the outer hull of the building: holes)
	2045	!-- of sorted polygon and put it back in polygon%vertices.
	2046	!-- Also add first vertex to the end of polygon and last vertex
	2047	!-- before the beginning of polygon.
	2048	DEALLOCATE(polygon%vertices)
	2049	ALLOCATE(polygon%vertices(0:nosv+1))
	2050	polygon%vertices(1:nosv) = sorted_p(1:nosv)
	2051	polygon%vertices(0) = sorted_p(nosv)
	2052	polygon%vertices(nosv+1) = sorted_p(1)
	2053	polygons(i_p)%nov = nosv
	2054	nov = polygons(i_p)%nov
	2055	DEALLOCATE(sorted_p)
	2056	ENDIF
	2057
	2058	END SUBROUTINE sort_polygon
	2059
	2060	!------------------------------------------------------------------------------!
	2061	! Description:
	2062	! ------------
	2063	!> Reduces the number of vertices in a polygon using the
	2064	!> Douglas-Poiker-Algorithm (1973)
	2065	!------------------------------------------------------------------------------!
	2066	RECURSIVE SUBROUTINE simplify_polygon( id_s, id_e, tol )
	2067
	2068	IMPLICIT NONE
	2069
	2070	INTEGER(iwp) :: max_dist_ind !< Index of vertex with maximum distance
	2071	INTEGER(iwp) :: il !< counter
	2072	INTEGER(iwp) :: id_s !< End index in polygon
	2073	INTEGER(iwp) :: id_e !< End index in polygon
	2074
	2075	REAL(wp) :: max_dist !< Maximum distance from line
	2076	REAL(wp) :: dum_dist !< Distance from line: dummy
	2077	REAL(wp) :: tol !< factor that determines how far a vertex can be from the polygon approximation so that the approximation is still accepted
	2078
	2079	max_dist = 0.
	2080	max_dist_ind = -999999
	2081	!
	2082	!-- Find vertex with max distance to id_s and id_e
	2083	DO il = id_s + 1, id_e -1
	2084	dum_dist = dist_point_to_edge(polygon%vertices(id_s)%x, &
	2085	polygon%vertices(id_s)%y, &
	2086	polygon%vertices(id_e)%x, &
	2087	polygon%vertices(id_e)%y, &
	2088	polygon%vertices(il)%x, &
	2089	polygon%vertices(il)%y)
	2090	IF ( dum_dist > max_dist ) THEN
	2091	max_dist = dum_dist
	2092	max_dist_ind = il
	2093	ENDIF
	2094	ENDDO
	2095
	2096	IF ( max_dist > tol ) THEN
	2097	CALL simplify_polygon( id_s, max_dist_ind, tol )
	2098	CALL simplify_polygon( max_dist_ind, id_e, tol )
	2099	ELSE
	2100	polygon%vertices(id_s+1:id_e-1)%delete = .TRUE.
	2101	ENDIF
	2102
	2103	END SUBROUTINE simplify_polygon
	2104
	2105	!------------------------------------------------------------------------------!
	2106	! Description:
	2107	! ------------
	2108	!> Checks if a vertex of a polygon is inside another polygon and if so, deletes
	2109	!> it. The check is done using the crossing number algorithm. If a straight
	2110	!> ray starting at a point crosses the borders of one polygon an odd
	2111	!> number of times, the point is inside that polygon.
	2112	!> This algorithm detects buildings that are completely surrounded by
	2113	!> another building. They can be deleted since they can never be navigated.
	2114	!> TODO: Maybe add a flag to turn this off and on as it might not be needed.
	2115	!> also, if the domain has buildings at all boundary points, there would
	2116	!> only be one giant building and everything in it deleted. So nothing
	2117	!> could be navigated. relevant?!
	2118	!------------------------------------------------------------------------------!
	2119	SUBROUTINE inside_other_polygon( i_p )
	2120
	2121	IMPLICIT NONE
	2122
	2123	LOGICAL :: exit_flag !< flag to exit loops if an odd crossing number was found for any of a polygons vertices
	2124
	2125	INTEGER(iwp) :: cn !< number of crossings
	2126	INTEGER(iwp) :: i_p !< index of current polygon
	2127	INTEGER(iwp) :: il !< index of tested polygon
	2128	INTEGER(iwp) :: nov_test !< no. of vertices of test-polygon
	2129	INTEGER(iwp) :: ref_vert !< vertex currently being tested if it is inside another polygon
	2130	INTEGER(iwp) :: test_edge !< index of edge being tested
	2131
	2132	REAL(wp) :: px !< x-coord of the point at the crossing of the ray and the vertex
	2133	REAL(wp) :: xe !< x-coordinate of end point of edge
	2134	REAL(wp) :: xr !< x-coordinate of reference point
	2135	REAL(wp) :: xs !< x-coordinate of start point of edge
	2136	REAL(wp) :: ye !< y-coordinate of end point of edge
	2137	REAL(wp) :: yr !< y-coordinate of reference point
	2138	REAL(wp) :: ys !< y-coordinate of start point of edge
	2139
	2140	TYPE(polygon_type), POINTER :: test_pol !< Polygon to be tested
	2141
	2142	exit_flag = .FALSE.
	2143	!
	2144	!-- Loop over all polygons other than the one being tested
	2145	DO il = 1, polygon_counter
	2146	IF ( il == i_p ) CYCLE
	2147	test_pol => polygons(il)
	2148	nov_test = polygons(il)%nov
	2149	!
	2150	!-- Inclusion test is done for every vertex of the polygon
	2151	DO ref_vert = 1, nov
	2152	cn = 0
	2153	xr = polygon%vertices(ref_vert)%x
	2154	yr = polygon%vertices(ref_vert)%y
	2155	!
	2156	!-- All edges of the every polygon il is tested for ray crossing
	2157	DO test_edge = 1, nov_test
	2158
	2159	!-- It is tested wether the current edge crosses a ray that extends
	2160	!-- from the current point to the right indefinitely.
	2161	!-- Check if start point of edge is lower than end point. If they
	2162	!-- are the same, ignore, since horizontal edges are excluded
	2163	IF ( test_pol%vertices(test_edge)%y < &
	2164	test_pol%vertices(test_edge+1)%y ) &
	2165	THEN
	2166	xs = test_pol%vertices(test_edge)%x
	2167	xe = test_pol%vertices(test_edge+1)%x
	2168	ys = test_pol%vertices(test_edge)%y
	2169	ye = test_pol%vertices(test_edge+1)%y
	2170	ELSEIF ( test_pol%vertices(test_edge)%y > &
	2171	test_pol%vertices(test_edge+1)%y ) &
	2172	THEN
	2173	xs = test_pol%vertices(test_edge+1)%x
	2174	xe = test_pol%vertices(test_edge)%x
	2175	ys = test_pol%vertices(test_edge+1)%y
	2176	ye = test_pol%vertices(test_edge)%y
	2177	ELSE
	2178	CYCLE
	2179	ENDIF
	2180	!
	2181	!-- Only such edges where the starting point of the edge is south of
	2182	!-- (or equal to) the reference point and the end point is north of
	2183	!-- it are relevant. Note: an edge includes its southern endpoint
	2184	!-- and excludes its northern endpoint.
	2185	IF ( .NOT. (ys <= yr .AND. ye > yr )) CYCLE
	2186	!
	2187	!-- Only edges that are crossed on the right side of the reference
	2188	!-- point are relevant, those on the left are ignored
	2189	IF ( xs <= xr .AND. xe <= xr ) CYCLE
	2190	IF ( ( xs <= xr .AND. xe >= xr ) .OR. &
	2191	( xs >= xr .AND. xe <= xr ) ) &
	2192	THEN
	2193	px = xe - (xe-xs)*(ye-yr)/(ye-ys)
	2194	IF ( px <= xr ) CYCLE
	2195	ENDIF
	2196	!
	2197	!-- If none of the previous if clauses were true, a crossing with
	2198	!-- an eligible edge was found and the count increases
	2199	cn = cn + 1
	2200	ENDDO
	2201	!
	2202	!-- If the number of crossings is odd, the point is inside another
	2203	!-- polyon. The polygon associated with the point will be deleted
	2204	IF ( MOD(cn, 2) /= 0 ) THEN
	2205	exit_flag = .TRUE.
	2206	EXIT
	2207	ENDIF
	2208	ENDDO
	2209	IF ( exit_flag ) EXIT
	2210	ENDDO
	2211	IF ( exit_flag ) polygon%vertices%delete = .TRUE.
	2212
	2213	END SUBROUTINE inside_other_polygon
	2214
	2215	!------------------------------------------------------------------------------!
	2216	! Description:
	2217	! ------------
	2218	!> Deletes thoses vertices that are marked for deletion (%delete flag) and
	2219	!> resizes the polygon
	2220	!------------------------------------------------------------------------------!
	2221	SUBROUTINE delete_extra_vertices (i_p)
	2222
	2223	IMPLICIT NONE
	2224
	2225	INTEGER(iwp) :: il !< local counter
	2226	INTEGER(iwp) :: vcounter !< vertex counter
	2227	INTEGER(iwp) :: i_p !< polygon ID
	2228
	2229	TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: dummy_pol !< Temporarily stores non-deleted vertices
	2230
	2231	ALLOCATE(dummy_pol(1:nov))
	2232	vcounter = 0
	2233	!
	2234	!-- Check all vertices and only keep those not marked for deletion
	2235	DO il = 1, nov
	2236	IF ( .NOT. polygon%vertices(il)%delete ) THEN
	2237	vcounter = vcounter + 1
	2238	dummy_pol(vcounter) = polygon%vertices(il)
	2239	ENDIF
	2240	ENDDO
	2241	!
	2242	!-- Set new number of vertices in the polygon
	2243	nov = vcounter
	2244	polygons(i_p)%nov = nov
	2245	!
	2246	!-- Resize
	2247	DEALLOCATE(polygon%vertices)
	2248	ALLOCATE(polygon%vertices(0:nov+1))
	2249	polygon%vertices(1:nov) = dummy_pol(1:nov)
	2250	polygon%vertices(0) = polygon%vertices(nov)
	2251	polygon%vertices(nov+1) = polygon%vertices(1)
	2252	DEALLOCATE(dummy_pol)
	2253
	2254	END SUBROUTINE delete_extra_vertices
	2255
	2256	!------------------------------------------------------------------------------!
	2257	! Description:
	2258	! ------------
	2259	!> Deletes polygons that contain no vertices (happens for those polygons that
	2260	!> were entirely encompassed by another polygon)
	2261	!------------------------------------------------------------------------------!
	2262	SUBROUTINE delete_empty_polygons
	2263
	2264	IMPLICIT NONE
	2265
	2266	INTEGER(iwp) :: il !< local counter
	2267	INTEGER(iwp) :: pc !< number of nonempty polygons
	2268	INTEGER(iwp) :: sv !< size of vertex array
	2269
	2270	TYPE(polygon_type), DIMENSION(:), ALLOCATABLE :: dummy_polygons !< temporarily stores non-deletd polygons
	2271
	2272	pc = 0
	2273	sv = 0
	2274	ALLOCATE( dummy_polygons(1:polygon_counter) )
	2275	!
	2276	!-- Keep only those polygons that contain any vertices, skip the rest
	2277	DO il = 1, polygon_counter
	2278	IF ( polygons(il)%nov > 0 ) THEN
	2279	pc = pc + 1
	2280	sv = SIZE(polygons(il)%vertices)
	2281	ALLOCATE(dummy_polygons(pc)%vertices(0:sv-1))
	2282	dummy_polygons(pc) = polygons(il)
	2283	ENDIF
	2284	ENDDO
	2285	polygon_counter = pc
	2286	!
	2287	!-- Resize polygon array
	2288	DEALLOCATE(polygons)
	2289	ALLOCATE(polygons(1:polygon_counter))
	2290	DO il = 1, polygon_counter
	2291	!
	2292	!-- give each %vertices array the correct size and information
	2293	sv = SIZE(dummy_polygons(il)%vertices)
	2294	polygons(il)%nov = sv - 2
	2295	ALLOCATE(polygons(il)%vertices(0:sv-1))
	2296	polygons(il) = dummy_polygons(il)
	2297	ENDDO
	2298	DEALLOCATE(dummy_polygons)
	2299
	2300	END SUBROUTINE delete_empty_polygons
	2301
	2302	END MODULE polygon_creation
	2303
	2304	MODULE mesh_creation
	2305
	2306	USE kinds
	2307
	2308	USE mod_functions
	2309
	2310	USE variables
	2311
	2312	CONTAINS
	2313
	2314	!------------------------------------------------------------------------------!
	2315	! Description:
	2316	! ------------
	2317	!> Creates the navigation mesh:
	2318	!> 1) Finds eligible vertices (those that are locally convex)
	2319	!> 2) Adds them to the mesh
	2320	!> 3) Adds connections between mesh points if they are in line of sight
	2321	!> of each other and the connecting line does not point into either of
	2322	!> the originating polygons (this is known as a visibility graph)
	2323	!------------------------------------------------------------------------------!
	2324	SUBROUTINE create_nav_mesh
	2325
	2326	IMPLICIT NONE
	2327
	2328	LOGICAL :: add !< flag for second cycle of add loop
	2329	LOGICAL :: intersection_found !< flag to indicate a found intersection
	2330
	2331	INTEGER(iwp) :: cmp !< counter: current mesh point
	2332	INTEGER(iwp) :: il !< local counter
	2333	INTEGER(iwp) :: jl !< local counter
	2334	INTEGER(iwp) :: pid !< polygon id of current mesh point
	2335	INTEGER(iwp) :: pid_t !< polygon id of tested mesh point
	2336	INTEGER(iwp) :: pl !< polygon counter
	2337	INTEGER(iwp) :: vid !< vertex id of current mesh point
	2338	INTEGER(iwp) :: vid_t !< vertex id of tested mesh point
	2339	INTEGER(iwp) :: vl !< vertex counter
[3259]	2340
[3159]	2341	REAL(wp) :: v1x !< x-coordinate of test vertex 1 for intersection test
	2342	REAL(wp) :: v1y !< y-coordinate of test vertex 1 for intersection test
	2343	REAL(wp) :: v2x !< x-coordinate of test vertex 2 for intersection test
	2344	REAL(wp) :: v2y !< y-coordinate of test vertex 2 for intersection test
	2345	REAL(wp) :: x !< x-coordinate of current mesh point
	2346	REAL(wp) :: x_t !< x-coordinate of tested mesh point
	2347	REAL(wp) :: y !< y-coordinate of current mesh point
	2348	REAL(wp) :: y_t !< y-coordinate of tested mesh point
	2349	REAL(wp) :: corner_x !< x-coordinate of shifted corner
	2350	REAL(wp) :: corner_x_e !< x-coordinate of end of corner gate
	2351	REAL(wp) :: corner_y !< y-coordinate of shifted corner
	2352	REAL(wp) :: corner_y_e !< y-coordinate of end of corner gate
	2353	REAL(wp) :: t_start !< CPU measure: start
	2354	REAL(wp) :: t_inter !< CPU measure: output test time
	2355	REAL(wp) :: t_inter1 !< CPU measure: output test time
	2356	REAL(wp) :: t_end !< CPU measure: end
	2357	REAL(wp) :: t_left !< CPU measure: estimate for time left
	2358	REAL(wp) :: t_done !< CPU measure: elapsed time
	2359	REAL(wp) :: percent_done !< CPU measure: proportion of mesh points checked
	2360
	2361	!
	2362	!-- Add all convex vertices to the mesh.
	2363	!-- DO loop will be executed twice. Once to count the mesh points to be
	2364	!-- added and allocate the mesh point array, the second time (add == .TRUE.)
	2365	!-- to fill the mesh point array.
[3400]	2366	WRITE(*,'(1X,A)') 'Adding polygon vertices to mesh ...'
[3159]	2367	add = .FALSE.
	2368	DO
	2369	cmp = 0
	2370	DO il = 1, polygon_counter
	2371	polygon => polygons(il)
	2372	nov = polygons(il)%nov
	2373	DO jl = 1, nov
	2374	!
	2375	!-- In a polygon that is sorted counter-clockwise, if the next vertex
	2376	!-- is left of the line connecting the previous and the current vertex,
	2377	!-- the current vertex is locally convex.
	2378	IF ( is_left(polygon%vertices(jl-1)%x,polygon%vertices(jl-1)%y, &
	2379	polygon%vertices(jl)%x,polygon%vertices(jl)%y, &
	2380	polygon%vertices(jl+1)%x,polygon%vertices(jl+1)%y) ) &
	2381	THEN
	2382
	2383	corner_x = polygon%vertices(jl)%x
	2384	corner_y = polygon%vertices(jl)%y
	2385	!
	2386	!-- Create end point for corner navigation
	2387	IF ( add ) THEN
	2388	CALL shift_corner_outward( &
	2389	polygon%vertices(jl-1)%x, polygon%vertices(jl-1)%y,&
	2390	polygon%vertices(jl+1)%x, polygon%vertices(jl+1)%y,&
	2391	polygon%vertices(jl)%x, polygon%vertices(jl)%y, &
	2392	corner_x_e, corner_y_e, 1._wp )
	2393	ENDIF
	2394	!
	2395	!-- Disregard corners outside of the domain
	2396	IF ( corner_x<=(nx+1)*dx .AND. corner_x>=0 .AND. &
	2397	corner_y<=(ny+1)*dy .AND. corner_y>=0) &
	2398	THEN
	2399	cmp = cmp + 1
	2400	IF ( add ) THEN
	2401	CALL set_mesh_point( mesh(cmp), il, jl, &
	2402	corner_x, corner_y, &
	2403	corner_x_e, corner_y_e )
	2404	ENDIF
	2405	ENDIF
	2406	ENDIF
	2407	ENDDO
	2408	ENDDO
	2409	IF ( add ) EXIT
	2410	add = .TRUE.
	2411	ALLOCATE( mesh(1:cmp) )
	2412	ENDDO
[3400]	2413	WRITE(*,'(6X,A,1X,I10,1X,A,/)') 'Done. Added',cmp,'vertices to mesh.'
	2414	WRITE(*,'(1X,A)') 'Establishing connections in mesh ...'
[3159]	2415	!
	2416	!-- CPU measurement
	2417	CALL CPU_TIME(t_start)
	2418	CALL CPU_TIME(t_inter)
	2419	DO il = 1, cmp
	2420	!-- Output status of processing
	2421	CALL CPU_TIME(t_inter1)
	2422	IF ( t_inter1 - t_inter > 4. ) THEN
	2423	t_done = (t_inter1-t_start)/60.
	2424	percent_done = REAL(il)/cmp*100.
	2425	t_left = t_done/percent_done*(100-percent_done)
	2426	WRITE(*,'(3X,2(A,I8),A,F6.2,2(A,F7.1),A,I10)') &
	2427	'Mesh point ',il,' of ' ,cmp, &
	2428	': ' ,percent_done, &
	2429	' % \|\| elapsed time : ' ,t_done, &
	2430	' min \|\| ETA: ' ,t_left, &
	2431	' min \|\| number of connections found: ',number_of_connections
	2432	CALL CPU_TIME(t_inter)
	2433	ENDIF
	2434	x = mesh(il)%x
	2435	y = mesh(il)%y
	2436	pid = mesh(il)%polygon_id
	2437	vid = mesh(il)%vertex_id
	2438	DO jl = 1, cmp
	2439	!
	2440	!-- No mesh point can be connected to itself
	2441	IF ( il == jl ) CYCLE
	2442	x_t = mesh(jl)%x
	2443	y_t = mesh(jl)%y
	2444	pid_t = mesh(jl)%polygon_id
	2445	vid_t = mesh(jl)%vertex_id
	2446	!
	2447	!-- Cycle, if a connection had already been established
	2448	IF ( ANY(mesh(jl)%connected_vertices == il) ) CYCLE
	2449	!
	2450	!-- If the distance between two nodes is larger than 600 m,
	2451	!-- no connection will be made since there will typically no be such
	2452	!-- long, straight ways in a city that a pedestrian will walk
	2453	IF ( SQRT((x_t-x)2 +(y_t-y)2) > 400. ) CYCLE
	2454	!
	2455	!-- If the connecting line between two mesh points points into either
	2456	!-- or both of the corresponding polygons, no connection will be
	2457	!-- established between the two points. This is the case if the
	2458	!-- previous (next) vertex of the polygon is right of the connecting
	2459	!-- line and the next (previous) vertex of the polygon is left of the
	2460	!-- connecting line. This is checked for both polygons.
	2461	IF ( ((is_left(x_t,y_t,x,y,polygons(pid)%vertices(vid-1)%x, &
	2462	polygons(pid)%vertices(vid-1)%y) &
	2463	.AND. is_right(x_t,y_t,x,y,polygons(pid)%vertices(vid+1)%x, &
	2464	polygons(pid)%vertices(vid+1)%y) ) &
	2465	.OR. (is_right(x_t,y_t,x,y,polygons(pid)%vertices(vid-1)%x, &
	2466	polygons(pid)%vertices(vid-1)%y) &
	2467	.AND. is_left(x_t,y_t,x,y,polygons(pid)%vertices(vid+1)%x, &
	2468	polygons(pid)%vertices(vid+1)%y)) ) &
	2469	.OR. ((is_left(x,y,x_t,y_t,polygons(pid_t)%vertices(vid_t-1)%x, &
	2470	polygons(pid_t)%vertices(vid_t-1)%y) &
	2471	.AND. is_right(x,y,x_t,y_t,polygons(pid_t)%vertices(vid_t+1)%x, &
	2472	polygons(pid_t)%vertices(vid_t+1)%y) ) &
	2473	.OR. (is_right(x,y,x_t,y_t,polygons(pid_t)%vertices(vid_t-1)%x, &
	2474	polygons(pid_t)%vertices(vid_t-1)%y) &
	2475	.AND. is_left(x,y,x_t,y_t,polygons(pid_t)%vertices(vid_t+1)%x, &
	2476	polygons(pid_t)%vertices(vid_t+1)%y)) ) ) &
	2477	THEN
	2478	CYCLE
	2479	ENDIF
	2480	!
	2481	!-- For each edge of each polygon, check if it intersects with the
	2482	!-- potential connection. If so, no connection can be made
	2483	!-- THIS IS THE BOTTLENECK OF THE PROGRAM
	2484	intersection_found = .FALSE.
	2485	DO pl = pid, polygon_counter
	2486	DO vl = 1, polygons(pl)%nov
	2487	v1x = polygons(pl)%vertices(vl)%x
	2488	v1y = polygons(pl)%vertices(vl)%y
	2489	v2x = polygons(pl)%vertices(vl+1)%x
	2490	v2y = polygons(pl)%vertices(vl+1)%y
	2491	intersection_found = intersect(x,y,x_t,y_t,v1x,v1y,v2x,v2y)
	2492	IF ( intersection_found ) EXIT
	2493	ENDDO
	2494	IF ( intersection_found ) EXIT
	2495	ENDDO
	2496	IF ( intersection_found ) CYCLE
	2497	DO pl = pid, 1, -1
	2498	IF ( pl == pid ) CYCLE
	2499	DO vl = 1, polygons(pl)%nov
	2500	v1x = polygons(pl)%vertices(vl)%x
	2501	v1y = polygons(pl)%vertices(vl)%y
	2502	v2x = polygons(pl)%vertices(vl+1)%x
	2503	v2y = polygons(pl)%vertices(vl+1)%y
	2504	intersection_found = intersect(x,y,x_t,y_t,v1x,v1y,v2x,v2y)
	2505	IF ( intersection_found ) EXIT
	2506	ENDDO
	2507	IF ( intersection_found ) EXIT
	2508	ENDDO
	2509	IF ( intersection_found ) CYCLE
	2510	!
	2511	!-- If neither of the above two test was true, a connection will be
	2512	!-- established between the two mesh points.
	2513	number_of_connections = number_of_connections + 1
	2514	CALL add_connection(mesh(il),jl, mesh(jl))
	2515	CALL add_connection(mesh(jl),il, mesh(il))
	2516	ENDDO
	2517	ENDDO
	2518	!
	2519	!-- Adapt connected_vertices arrays
	2520	DO il = 1, cmp
	2521	CALL reduce_connections(mesh(il))
	2522	ENDDO
	2523	CALL CPU_TIME(t_end)
	2524	!
	2525	!-- Output to terminal
	2526	WRITE(*,'(6X,A,I10,A)') 'Done. Established ',number_of_connections, &
	2527	' connections in mesh'
	2528	WRITE(*,'(6X,A,F10.1,A)') 'Time needed for calculation: ', &
	2529	t_end-t_start,' seconds'
[3216]	2530
[3159]	2531	END SUBROUTINE create_nav_mesh
	2532
	2533	!------------------------------------------------------------------------------!
	2534	! Description:
	2535	! ------------
	2536	!> Initializes a point of the navigation mesh
	2537	!------------------------------------------------------------------------------!
	2538	SUBROUTINE set_mesh_point (in_mp,pid,vid,x,y,x_s,y_s)
	2539
	2540	IMPLICIT NONE
	2541
	2542	INTEGER(iwp) :: pid !< polygon ID
	2543	INTEGER(iwp) :: vid !< vertex ID
	2544
	2545	REAL(wp) :: x !< x-value of mesh point for path calculation
	2546	REAL(wp) :: x_s !< x-value shifted outward from corner
	2547	REAL(wp) :: y !< y-value of mesh point for path calculation
	2548	REAL(wp) :: y_s !< y-value shifted outward from corner
	2549
	2550	TYPE(mesh_point) :: in_mp !< mesh point to be created
	2551
	2552	in_mp%origin_id = -1
	2553	in_mp%polygon_id = pid
	2554	in_mp%vertex_id = vid
	2555	in_mp%cost_so_far = 1.d12
	2556	in_mp%x = x
	2557	in_mp%y = y
	2558	in_mp%x_s = x_s
	2559	in_mp%y_s = y_s
	2560	in_mp%noc = 0
	2561
	2562	ALLOCATE(in_mp%connected_vertices(1:100), &
	2563	in_mp%distance_to_vertex(1:100))
	2564
	2565	in_mp%connected_vertices = -999
	2566	in_mp%distance_to_vertex = -999.
	2567
	2568	END SUBROUTINE set_mesh_point
	2569
	2570	!------------------------------------------------------------------------------!
	2571	! Description:
	2572	! ------------
	2573	!> Shifts a corner (middle one of three consecutive points a, b and p) outward
	2574	!> by a given length along the angle bisector. Stores the result to res_x/res_y
	2575	!------------------------------------------------------------------------------!
	2576	SUBROUTINE shift_corner_outward ( a_x, a_y, b_x, b_y, p_x, p_y, res_x, &
	2577	res_y, shift )
	2578
	2579	IMPLICIT NONE
	2580
	2581	REAL(wp) :: a_x !< x-value of point A
	2582	REAL(wp) :: a_y !< y-value of point A
	2583	REAL(wp) :: abs_ap !< distance from A to P
	2584	REAL(wp) :: abs_bp !< distance from B to P
	2585	REAL(wp) :: abs_co !< length of angle bisector
	2586	REAL(wp) :: b_x !< x-value of point B
	2587	REAL(wp) :: b_y !< y-value of point B
	2588	REAL(wp) :: eap_x !< x-value of unit vector from A to P
	2589	REAL(wp) :: eap_y !< y-value of unit vector from A to P
	2590	REAL(wp) :: ebp_x !< x-value of unit vector from B to P
	2591	REAL(wp) :: ebp_y !< y-value of unit vector from B to P
	2592	REAL(wp) :: p_x !< x-value of point P
	2593	REAL(wp) :: p_y !< y-value of point P
	2594	REAL(wp) :: res_x !< x-value of result
	2595	REAL(wp) :: res_y !< y-value of result
	2596	REAL(wp) :: shift !< distance of shift in meters
	2597
	2598	!
	2599	!-- Get unit vector from previous to current vertex
	2600	eap_x = p_x - a_x
	2601	eap_y = p_y - a_y
	2602	abs_ap = SQRT(eap_x2+eap_y2)
	2603	eap_x = eap_x/abs_ap
	2604	eap_y = eap_y/abs_ap
	2605	!
	2606	!-- Get unit vector from next to current vertex
	2607	ebp_x = p_x - b_x
	2608	ebp_y = p_y - b_y
	2609	abs_bp = SQRT(ebp_x2+ebp_y2)
	2610	ebp_x = ebp_x/abs_bp
	2611	ebp_y = ebp_y/abs_bp
	2612	!
	2613	!-- Add previous two vectors to get angle bisector of corner.
	2614	!-- Then, set its length to shift and add to original vertex
	2615	!-- vector to shift it outward
	2616	res_x = eap_x + ebp_x
	2617	res_y = eap_y + ebp_y
	2618	abs_co = SQRT(res_x2+res_y2)
	2619	res_x = shift*res_x/abs_co + p_x
	2620	res_y = shift*res_y/abs_co + p_y
	2621
	2622	END SUBROUTINE shift_corner_outward
	2623
	2624	!------------------------------------------------------------------------------!
	2625	! Description:
	2626	! ------------
	2627	!> Adds a connection between two points of the navigation mesh
	2628	!> (one-way: in_mp1 to in_mp2)
	2629	!------------------------------------------------------------------------------!
	2630	SUBROUTINE add_connection (in_mp1,id2,in_mp2)
	2631
	2632	IMPLICIT NONE
	2633
	2634	LOGICAL :: connection_established !< Flag to indicate if connection has already been established
	2635
	2636	INTEGER(iwp) :: id2 !< ID of in_mp2
	2637	INTEGER(iwp) :: il !< local counter
	2638	INTEGER(iwp) :: noc1 !< number of connections in in_mp1
	2639
	2640	INTEGER, DIMENSION(:), ALLOCATABLE :: dum_cv !< dummy array for connected_vertices
	2641
	2642	REAL(wp) :: dist !< Distance between the two points
	2643
	2644	REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_dtv
	2645
	2646	TYPE(mesh_point) :: in_mp1 !< mesh point that gets a new connection
	2647	TYPE(mesh_point) :: in_mp2 !< mesh point in_mp1 will be connected to
	2648
	2649	connection_established = .FALSE.
	2650	!
	2651	!-- Check if connection has already been established
	2652	noc1 = SIZE(in_mp1%connected_vertices)
	2653	DO il = 1, in_mp1%noc
	2654	IF ( in_mp1%connected_vertices(il) == id2 ) THEN
	2655	connection_established = .TRUE.
	2656	EXIT
	2657	ENDIF
	2658	ENDDO
	2659
	2660	IF ( .NOT. connection_established ) THEN
	2661	!
	2662	!-- Resize arrays, if necessary
	2663	IF ( in_mp1%noc >= noc1 ) THEN
	2664	ALLOCATE( dum_cv(1:noc1),dum_dtv(1:noc1) )
	2665	dum_cv = in_mp1%connected_vertices
	2666	dum_dtv = in_mp1%distance_to_vertex
	2667	DEALLOCATE( in_mp1%connected_vertices, in_mp1%distance_to_vertex )
	2668	ALLOCATE( in_mp1%connected_vertices(1:2*noc1), &
	2669	in_mp1%distance_to_vertex(1:2*noc1) )
	2670	in_mp1%connected_vertices = -999
	2671	in_mp1%distance_to_vertex = -999.
	2672	in_mp1%connected_vertices(1:noc1) = dum_cv
	2673	in_mp1%distance_to_vertex(1:noc1) = dum_dtv
	2674	ENDIF
	2675
	2676	!
	2677	!-- Add connection
	2678	in_mp1%noc = in_mp1%noc+1
	2679	dist = SQRT( (in_mp1%x - in_mp2%x)2 + (in_mp1%y - in_mp2%y)2 )
	2680	in_mp1%connected_vertices(in_mp1%noc) = id2
	2681	in_mp1%distance_to_vertex(in_mp1%noc) = dist
	2682	ENDIF
	2683
	2684	END SUBROUTINE add_connection
	2685
	2686	!------------------------------------------------------------------------------!
	2687	! Description:
	2688	! ------------
	2689	!> Reduces the size of connection array to the amount of actual connections
	2690	!> after all connetions were added
	2691	!------------------------------------------------------------------------------!
	2692	SUBROUTINE reduce_connections (in_mp)
	2693
	2694	IMPLICIT NONE
	2695
	2696	INTEGER(iwp) :: noc !< Number of connections
	2697
	2698	INTEGER, DIMENSION(:), ALLOCATABLE :: dum_cv !< dummy: connected_vertices
	2699
	2700	REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_dtv !< dummy: distance_to_vertex
	2701
	2702	TYPE(mesh_point) :: in_mp !< Input mesh point
	2703
	2704	noc = in_mp%noc
	2705	ALLOCATE( dum_cv(1:noc),dum_dtv(1:noc) )
	2706	dum_cv = in_mp%connected_vertices(1:noc)
	2707	dum_dtv = in_mp%distance_to_vertex(1:noc)
	2708	DEALLOCATE( in_mp%connected_vertices, in_mp%distance_to_vertex )
	2709	ALLOCATE( in_mp%connected_vertices(1:noc), &
	2710	in_mp%distance_to_vertex(1:noc) )
	2711	in_mp%connected_vertices(1:noc) = dum_cv(1:noc)
	2712	in_mp%distance_to_vertex(1:noc) = dum_dtv(1:noc)
	2713
	2714	END SUBROUTINE reduce_connections
	2715
	2716	!------------------------------------------------------------------------------!
	2717	! Description:
	2718	! ------------
	2719	!> Writes all NavMesh information into binary file and building data to ASCII
	2720	!------------------------------------------------------------------------------!
	2721	SUBROUTINE bin_out_mesh
	2722
	2723	IMPLICIT NONE
	2724
	2725	INTEGER(iwp) :: il !< local counter
	2726	INTEGER(iwp) :: jl !< local counter
	2727	INTEGER(iwp) :: size_of_mesh !< size of mesh
	2728	INTEGER(iwp) :: size_of_pols !< size of polygon
	2729
[3400]	2730	WRITE(*,'(1X,A)') 'Writing binary output data ...'
[3159]	2731
	2732	OPEN ( 14, FILE= TRIM(runname)//'_nav', FORM='UNFORMATTED', STATUS='replace' )
	2733	!
	2734	!-- Output of mesh data
	2735	size_of_mesh = SIZE(mesh)
	2736	WRITE(14) size_of_mesh
	2737	DO il = 1, size_of_mesh
	2738	WRITE(14) mesh(il)%polygon_id, mesh(il)%vertex_id, mesh(il)%noc, &
	2739	mesh(il)%origin_id, mesh(il)%cost_so_far, mesh(il)%x, &
	2740	mesh(il)%y, mesh(il)%x_s, mesh(il)%y_s
	2741	DO jl = 1, mesh(il)%noc
	2742	WRITE(14) mesh(il)%connected_vertices(jl), &
	2743	mesh(il)%distance_to_vertex(jl)
	2744	ENDDO
	2745	ENDDO
	2746	!
	2747	!-- Output of building polygon data
	2748	size_of_pols = SIZE(polygons)
	2749	WRITE(14) size_of_pols
	2750	DO il = 1, size_of_pols
	2751	WRITE(14) polygons(il)%nov
	2752	DO jl = 0, polygons(il)%nov+1
	2753	WRITE(14) polygons(il)%vertices(jl)%delete, &
	2754	polygons(il)%vertices(jl)%x, polygons(il)%vertices(jl)%y
	2755	ENDDO
	2756	ENDDO
	2757	CLOSE(14)
	2758	!
	2759	!-- Output building data to ASCII file
	2760	OPEN(UNIT=7,FILE='topo.txt',STATUS='replace',ACTION='write')
	2761	DO i = 1, polygon_counter
	2762	IF (polygons(i)%nov == 0) CYCLE
	2763	DO j = 1, polygons(i)%nov
	2764	WRITE(7,150) i,j,polygons(i)%vertices(j)%x, &
	2765	polygons(i)%vertices(j)%y
	2766	ENDDO
	2767	ENDDO
	2768	CLOSE(7)
	2769
	2770	WRITE(*,'(6X,A)') 'Done, tool terminating.', ' ', &
	2771	'Before starting your PALM run, please check the', &
	2772	'ASCII file topo.txt to see if you are satisfied', &
	2773	'with the polygon representation of the building', &
	2774	'data. If not, consider adjusting the parameter', &
	2775	'tolerance_dp accordingly.', ' ', 'Bye, Bye!', ' '
	2776	CALL CPU_TIME(finish)
[3400]	2777	WRITE(*,'(1X,A,F10.4,A)') 'Total runtime: ', finish-start, ' seconds'
[3159]	2778
[3400]	2779	150 FORMAT (2(I7,1X),2(F9.2,1X) )
[3159]	2780
	2781	END SUBROUTINE bin_out_mesh
	2782
	2783	END MODULE mesh_creation
	2784
	2785	PROGRAM nav_mesh
	2786
	2787	USE mesh_creation
	2788	USE polygon_creation
	2789	USE variables
	2790	IMPLICIT NONE
	2791
	2792
	2793	!
	2794	!-- Start CPU mesurement
	2795	CALL CPU_TIME(start)
	2796	!
	2797	!-- Initialization
	2798	CALL init
	2799
	2800	WRITE(,) "Converting building data to polygons ..."
	2801	!
	2802	!-- Convert gridded building data to polygons
	2803	CALL identify_polygons
	2804	!
	2805	!-- Find corners in topography and add them to polygons
	2806	CALL identify_corners
	2807	!
	2808	!-- Sort polygons counter-clockwise, then simplify them
	2809	DO i = 1, polygon_counter
	2810	polygon => polygons(i)
	2811	nov = polygons(i)%nov
	2812	CALL sort_polygon(i)
	2813	!
	2814	!-- Simplify each polygon using douglas-peucker algorithm. If the number
	2815	!-- of vertices would fall below 4 due to this procedure, the tolerance
	2816	!-- for the algorithm is reduced and it is run again.
[3198]	2817	DO i_sc = 0, 2
[3159]	2818	CALL simplify_polygon(1,nov+1,tolerance_dp(i_sc))
	2819	i_cn = 0
	2820	DO j = 1, nov
	2821	IF ( .NOT. polygon%vertices(j)%delete ) i_cn = i_cn + 1
	2822	ENDDO
	2823	IF ( i_cn > 3 ) THEN
	2824	EXIT
	2825	ELSE
	2826	polygon%vertices(:)%delete = .FALSE.
	2827	ENDIF
	2828	ENDDO
	2829	CALL delete_extra_vertices(i)
	2830	ENDDO
	2831	!
	2832	!-- Remove buildings that are surrounded by another building
	2833	IF ( .NOT. internal_buildings ) THEN
	2834	DO i = 1, polygon_counter
	2835	polygon => polygons(i)
	2836	nov = polygons(i)%nov
	2837	CALL inside_other_polygon(i)
	2838	ENDDO
	2839	ENDIF
	2840	!
	2841	!-- Delete vertices that are marked for deletion
	2842	DO i = 1, polygon_counter
	2843	polygon => polygons(i)
	2844	nov = polygons(i)%nov
	2845	CALL delete_extra_vertices(i)
	2846	ENDDO
	2847	!
	2848	!-- Count number of vertices
	2849	vertex_counter = 0
	2850	DO i = 1, polygon_counter
	2851	polygon => polygons(i)
	2852	nov = polygons(i)%nov
	2853	vertex_counter = vertex_counter + nov
	2854	ENDDO
	2855	!
	2856	!-- Delete polygons with no vertices
	2857	CALL delete_empty_polygons
[3400]	2858	WRITE(*,'(2(6X,A,I10,1X,A,/))') &
[3159]	2859	'Done. Created a total of', polygon_counter, 'polygon(s)', &
	2860	' with a total of', vertex_counter, 'vertices'
	2861	!
	2862	!-- Crate Navigation mesh from polygon data
	2863	CALL create_nav_mesh
	2864	!
	2865	!-- Binary mesh output
	2866	CALL bin_out_mesh
	2867
	2868	END PROGRAM nav_mesh

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |