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

Last change on this file since 3208 was 3208, checked in by sward, 3 years ago

Renamed nav_mesh to agent_preprocessing and added it to palmbuild

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