Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |