Home

Context Navigation

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

Last change on this file since 3232 was 3216, checked in by sward, 6 years ago
Bugfix for gfortran: reordering of type definitions
Property svn:keywords set to `Id`
File size: 116.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |