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

Last change on this file since 3698 was 3665, checked in by raasch, 6 years ago

dummy statements added to avoid compiler warnings about unused variables, unused variables removed, ssh-call for submitting batch jobs on remote systems modified again to avoid output of login messages on specific systems

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