Home

Context Navigation

nav_mesh.f90 @ 3201

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

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

Context Navigation

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

Download in other formats: