Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |