Home

Context Navigation

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

Last change on this file since 3183 was 3168, checked in by sward, 6 years ago
Updated NetCDF variable names
Property svn:keywords set to `Id`
File size: 116.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |