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

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

Minor bugfix in agent_preprocessing

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