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

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

Added MAS end time, used time_since_reference_point, corrected tolerance_dp in nav_mesh

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