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

Last change on this file since 3301 was 3259, checked in by sward, 6 years ago

Updated agent preprocessing to avoid gfortran warnings

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