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

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

Updated NetCDF variable names

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