Home

Context Navigation

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

Last change on this file since 3209 was 3208, checked in by sward, 6 years ago
Renamed nav_mesh to agent_preprocessing and added it to palmbuild
Property svn:keywords set to `Id`
File size: 116.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |