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