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

Last change on this file since 3251 was 3216, checked in by sward, 6 years ago

Bugfix for gfortran: reordering of type definitions

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