Home

Context Navigation

source: palm/trunk/UTIL/nav_mesh/nav_mesh.f90 @ 3167

Last change on this file since 3167 was 3159, checked in by sward, 6 years ago
Added multi agent system
File size: 115.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |