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 knoop $ |
---|
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 | |
---|
413 | 11 WRITE(*,'(2A)') 'errors in file ',input_filename |
---|
414 | |
---|
415 | 12 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 |
---|
1333 | 10 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 |
---|