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

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

Added multi agent system

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