Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |