Home

Context Navigation

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

Last change on this file since 3356 was 3259, checked in by sward, 6 years ago
Updated agent preprocessing to avoid gfortran warnings
Property svn:keywords set to `Id`
File size: 116.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |