Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |