[3159] | 1 | !> @file multi_agent_system_mod.f90 |
---|
| 2 | !--------------------------------------------------------------------------------! |
---|
| 3 | ! This file is part of PALM-4U. |
---|
| 4 | ! |
---|
| 5 | ! PALM-4U 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-4U 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 2016 Leibniz Universitaet Hannover |
---|
| 18 | !------------------------------------------------------------------------------! |
---|
| 19 | ! |
---|
| 20 | ! Current revisions: |
---|
| 21 | ! ------------------ |
---|
| 22 | ! |
---|
| 23 | ! |
---|
| 24 | ! Former revisions: |
---|
| 25 | ! ----------------- |
---|
| 26 | ! $Id: multi_agent_system_mod.f90 3268 2018-09-21 13:45:37Z sward $ |
---|
[3268] | 27 | ! Cleaned up agent pointer assignment in output routine |
---|
| 28 | ! |
---|
| 29 | ! 3248 2018-09-14 09:42:06Z sward |
---|
[3248] | 30 | ! Minor formating changes |
---|
| 31 | ! |
---|
| 32 | ! 3246 2018-09-13 15:14:50Z sward |
---|
[3246] | 33 | ! Added error handling for input namelist via parin_fail_message |
---|
| 34 | ! |
---|
| 35 | ! 3241 2018-09-12 15:02:00Z raasch |
---|
[3241] | 36 | ! unused variables removed |
---|
| 37 | ! |
---|
| 38 | ! 3235 2018-09-07 14:06:15Z sward |
---|
[3235] | 39 | ! Bugfix in output, added agent-number dimension and related messages and |
---|
| 40 | ! input parameters, updated cpu logging, added mas_last_actions |
---|
| 41 | ! |
---|
| 42 | ! 3201 2018-08-20 11:45:01Z sward |
---|
[3201] | 43 | ! Bugfix, missing pre-processor directive. Set default |
---|
| 44 | ! read_agents_from_restartfile = .FALSE. restarts not yet implemented. |
---|
| 45 | ! |
---|
| 46 | ! 3198 2018-08-15 09:23:10Z sward |
---|
[3198] | 47 | ! Now using time_since_reference_point; moved multi_agent_system_start and |
---|
| 48 | ! multi_agent_system_end to control_parameters; renamed NAMELIST agents_par to |
---|
| 49 | ! agent_parameters |
---|
| 50 | ! |
---|
| 51 | ! 3187 2018-07-31 10:32:34Z sward |
---|
[3187] | 52 | ! Reworked agent pathfinding to avoid collisions with walls |
---|
| 53 | ! |
---|
| 54 | ! 3165 2018-07-24 13:12:42Z sward |
---|
[3165] | 55 | ! Added agent ID output |
---|
| 56 | ! |
---|
| 57 | ! 3160 2018-07-20 11:52:56Z sward |
---|
[3160] | 58 | ! Changed C_SIZEOF to STORAGE_SIZE |
---|
| 59 | ! |
---|
| 60 | ! 3159 2018-07-20 11:20:01Z sward |
---|
[3159] | 61 | ! Initial revision |
---|
| 62 | ! |
---|
| 63 | ! |
---|
| 64 | ! |
---|
| 65 | ! Authors: |
---|
| 66 | ! -------- |
---|
| 67 | ! @author sward |
---|
| 68 | ! |
---|
| 69 | ! |
---|
| 70 | ! Description: |
---|
| 71 | ! ------------ |
---|
| 72 | !> Multi Agent System for the simulation of pedestrian movement in urban |
---|
| 73 | !> environments |
---|
| 74 | !------------------------------------------------------------------------------! |
---|
| 75 | MODULE multi_agent_system_mod |
---|
| 76 | |
---|
| 77 | USE, INTRINSIC :: ISO_C_BINDING |
---|
| 78 | |
---|
| 79 | USE constants, & |
---|
| 80 | ONLY: pi |
---|
| 81 | |
---|
| 82 | USE control_parameters, & |
---|
[3198] | 83 | ONLY: dt_3d, message_string, time_since_reference_point, dt_write_agent_data |
---|
[3159] | 84 | |
---|
| 85 | USE cpulog, & |
---|
| 86 | ONLY: cpu_log, log_point, log_point_s |
---|
| 87 | |
---|
| 88 | USE grid_variables, & |
---|
| 89 | ONLY: ddx, ddy, dx, dy |
---|
| 90 | |
---|
| 91 | USE indices, & |
---|
| 92 | ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, & |
---|
| 93 | wall_flags_0 |
---|
| 94 | |
---|
| 95 | USE random_function_mod, & |
---|
| 96 | ONLY: random_function |
---|
| 97 | |
---|
| 98 | USE kinds |
---|
| 99 | |
---|
| 100 | USE pegrid |
---|
| 101 | |
---|
| 102 | CHARACTER(LEN=15) :: bc_mas_lr = 'absorb' !< left/right boundary condition |
---|
| 103 | CHARACTER(LEN=15) :: bc_mas_ns = 'absorb' !< north/south boundary condition |
---|
| 104 | |
---|
[3235] | 105 | INTEGER(iwp) :: deleted_agents = 0 !< number of deleted agents per time step |
---|
| 106 | INTEGER(iwp) :: dim_size_agtnum_manual = 9999999 !< namelist parameter (see documentation) |
---|
| 107 | INTEGER(iwp) :: heap_count !< number of items in binary heap (for pathfinding) |
---|
| 108 | INTEGER(iwp) :: ibc_mas_lr !< agent left/right boundary condition dummy |
---|
| 109 | INTEGER(iwp) :: ibc_mas_ns !< agent north/south boundary condition dummy |
---|
[3241] | 110 | ! INTEGER(iwp) :: ind_pm10 = -9 !< chemical species index of PM10 |
---|
| 111 | ! INTEGER(iwp) :: ind_pm25 = -9 !< chemical species index of PM2.5 |
---|
[3235] | 112 | INTEGER(iwp) :: iran_agent = -1234567 !< number for random generator |
---|
| 113 | INTEGER(iwp) :: min_nr_agent = 2 !< namelist parameter (see documentation) |
---|
| 114 | INTEGER(iwp) :: ghla_count_recv !< number of agents in left ghost layer |
---|
| 115 | INTEGER(iwp) :: ghna_count_recv !< number of agents in north ghost layer |
---|
| 116 | INTEGER(iwp) :: ghra_count_recv !< number of agents in right ghost layer |
---|
| 117 | INTEGER(iwp) :: ghsa_count_recv !< number of agents in south ghost layer |
---|
| 118 | INTEGER(iwp) :: maximum_number_of_agents = 0 !< maximum number of agents during run |
---|
| 119 | INTEGER(iwp) :: nr_move_north !< number of agts to move north during exchange_horiz |
---|
| 120 | INTEGER(iwp) :: nr_move_south !< number of agts to move south during exchange_horiz |
---|
| 121 | INTEGER(iwp) :: number_of_agents = 0 !< number of agents for each grid box (3d array is saved on agt_count) |
---|
| 122 | INTEGER(iwp) :: number_of_agent_groups = 1 !< namelist parameter (see documentation) |
---|
| 123 | INTEGER(iwp) :: sort_count_mas = 0 !< counter for sorting agents |
---|
| 124 | INTEGER(iwp) :: agt_path_size = 15 !< size of agent path array |
---|
| 125 | INTEGER(iwp) :: step_dealloc_mas = 100 !< namelist parameter (see documentation) |
---|
| 126 | INTEGER(iwp) :: total_number_of_agents !< total number of agents in the whole model domain |
---|
[3159] | 127 | |
---|
| 128 | INTEGER(iwp), PARAMETER :: NR_2_direction_move = 10000 !< parameter for agent exchange |
---|
| 129 | INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !< phase parameter |
---|
| 130 | INTEGER(iwp), PARAMETER :: PHASE_RELEASE = 2 !< phase parameter |
---|
| 131 | |
---|
| 132 | INTEGER(iwp), PARAMETER :: max_number_of_agent_groups = 100 !< maximum allowed number of agent groups |
---|
| 133 | |
---|
[3235] | 134 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: agt_count !< 3d array of number of agents of every grid box |
---|
[3159] | 135 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: s_measure_height !< k-index(s-grid) for measurement |
---|
[3235] | 136 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: top_top_s !< k-index of first s-gridpoint above topography |
---|
| 137 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: top_top_w !< k-index of first v-gridpoint above topography |
---|
| 138 | INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: obstacle_flags !< flags to identify corners and edges of topography that cannot be crossed by agents |
---|
[3159] | 139 | |
---|
[3235] | 140 | LOGICAL :: deallocate_memory_mas = .TRUE. !< namelist parameter (see documentation) |
---|
| 141 | LOGICAL :: dt_3d_reached_mas !< flag: agent timestep has reached model timestep |
---|
| 142 | LOGICAL :: dt_3d_reached_l_mas !< flag: agent timestep has reached model timestep |
---|
| 143 | LOGICAL :: agents_active = .FALSE. !< flag for agent system |
---|
| 144 | LOGICAL :: random_start_position_agents = .TRUE. !< namelist parameter (see documentation) |
---|
[3201] | 145 | LOGICAL :: read_agents_from_restartfile = .FALSE. !< namelist parameter (see documentation) |
---|
[3235] | 146 | LOGICAL :: agent_own_timestep = .FALSE. !< namelist parameter (see documentation) |
---|
[3159] | 147 | |
---|
| 148 | LOGICAL, DIMENSION(max_number_of_agent_groups) :: a_rand_target = .FALSE. !< namelist parameter (see documentation) |
---|
| 149 | |
---|
[3198] | 150 | REAL(wp) :: agent_maximum_age = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 151 | REAL(wp) :: alloc_factor_mas = 20.0_wp !< namelist parameter (see documentation) |
---|
| 152 | REAL(wp) :: coll_t_0 = 3. !< namelist parameter (see documentation) |
---|
| 153 | REAL(wp) :: corner_gate_start = 0.5_wp !< namelist parameter (see documentation) |
---|
[3235] | 154 | REAL(wp) :: corner_gate_width = 1.0_wp !< namelist parameter (see documentation) |
---|
| 155 | REAL(wp) :: dim_size_factor_agtnum = 1.0_wp !< namelist parameter (see documentation) |
---|
[3198] | 156 | REAL(wp) :: d_sigma_rep_agent !< inverse of sigma_rep_agent |
---|
| 157 | REAL(wp) :: d_sigma_rep_wall !< inverse of sigma_rep_wall |
---|
| 158 | REAL(wp) :: d_tau_accel_agent !< inverse of tau_accel_agent |
---|
| 159 | REAL(wp) :: desired_speed = 1.2_wp !< namelist parameter (see documentation) |
---|
| 160 | REAL(wp) :: des_sp_sig = .2_wp !< namelist parameter (see documentation) |
---|
| 161 | REAL(wp) :: dist_target_reached = 2.0_wp !< distance at which target counts as reached |
---|
| 162 | REAL(wp) :: dist_to_int_target = .25_wp !< namelist parameter (see documentation) |
---|
| 163 | REAL(wp) :: dt_agent = 0.02_wp !< namelist parameter (see documentation) |
---|
| 164 | REAL(wp) :: dt_arel = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 165 | REAL(wp) :: end_time_arel = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 166 | REAL(wp) :: force_x !< dummy value for force on current agent in x-direction |
---|
| 167 | REAL(wp) :: force_y !< dummy value for force on current agent in y-direction |
---|
| 168 | REAL(wp) :: max_dist_from_path = 0.25_wp !< distance from current path at which a new path is calculated |
---|
| 169 | REAL(wp) :: radius_agent = .25_wp !< namelist parameter (see documentation) |
---|
| 170 | REAL(wp) :: repuls_agent = 1.5_wp !< namelist parameter (see documentation) |
---|
| 171 | REAL(wp) :: repuls_wall = 7.0_wp !< namelist parameter (see documentation) |
---|
| 172 | REAL(wp) :: scan_radius_agent = 3.0_wp !< namelist parameter (see documentation) |
---|
| 173 | REAL(wp) :: scan_radius_wall = 2.0_wp !< namelist parameter (see documentation) |
---|
| 174 | REAL(wp) :: sigma_rep_agent = 0.3_wp !< namelist parameter (see documentation) |
---|
| 175 | REAL(wp) :: sigma_rep_wall = 0.1_wp !< namelist parameter (see documentation) |
---|
| 176 | REAL(wp) :: tau_accel_agent = 0.5_wp !< namelist parameter (see documentation) |
---|
| 177 | REAL(wp) :: time_arel = 0.0_wp !< time for agent release |
---|
| 178 | REAL(wp) :: time_write_agent_data = 0.0_wp !< write agent data at current time on file |
---|
| 179 | REAL(wp) :: v_max_agent = 1.3_wp !< namelist parameter (see documentation) |
---|
[3159] | 180 | |
---|
[3187] | 181 | REAL(wp), DIMENSION(:), ALLOCATABLE :: dummy_path_x !< dummy path (x-coordinate) |
---|
| 182 | REAL(wp), DIMENSION(:), ALLOCATABLE :: dummy_path_y !< dummy path (y-coordinate) |
---|
| 183 | |
---|
[3159] | 184 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: adx = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 185 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: ady = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 186 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: asl = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 187 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: asn = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 188 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: asr = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 189 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: ass = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 190 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: at_x = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 191 | REAL(wp), DIMENSION(max_number_of_agent_groups) :: at_y = 9999999.9_wp !< namelist parameter (see documentation) |
---|
| 192 | ! |
---|
| 193 | !-- Type for the definition of an agent |
---|
| 194 | TYPE agent_type |
---|
[3235] | 195 | INTEGER(iwp) :: block_nr !< number for sorting |
---|
| 196 | INTEGER(iwp) :: group !< number of agent group |
---|
| 197 | INTEGER(idp) :: id !< particle ID (64 bit integer) |
---|
| 198 | INTEGER(iwp) :: path_counter !< current target along path (path_x/y) |
---|
| 199 | LOGICAL :: agent_mask !< if this parameter is set to false the agent will be deleted |
---|
| 200 | REAL(wp) :: age !< age of agent |
---|
| 201 | REAL(wp) :: age_m !< age of agent |
---|
| 202 | REAL(wp) :: dt_sum !< sum of agents subtimesteps |
---|
| 203 | REAL(wp) :: clo !< clothing index |
---|
| 204 | REAL(wp) :: energy_storage !< energy stored by agent |
---|
| 205 | REAL(wp) :: force_x !< force term x-direction |
---|
| 206 | REAL(wp) :: force_y !< force term y-direction |
---|
| 207 | REAL(wp) :: origin_x !< origin x-position of agent |
---|
| 208 | REAL(wp) :: origin_y !< origin y-position of agent |
---|
| 209 | REAL(wp) :: pm10 !< PM10 concentration at agent position |
---|
| 210 | REAL(wp) :: pm25 !< PM25 concentration at agent position |
---|
| 211 | REAL(wp) :: speed_abs !< absolute value of agent speed |
---|
| 212 | REAL(wp) :: speed_e_x !< normalized speed of agent in x |
---|
| 213 | REAL(wp) :: speed_e_y !< normalized speed of agent in y |
---|
| 214 | REAL(wp) :: speed_des !< agent's desired speed |
---|
| 215 | REAL(wp) :: speed_x !< speed of agent in x |
---|
| 216 | REAL(wp) :: speed_y !< speed of agent in y |
---|
| 217 | REAL(wp) :: thermal_index !< the dynamic thermal index |
---|
| 218 | REAL(wp) :: windspeed !< absolute value of windspeed at agent position |
---|
| 219 | REAL(wp) :: x !< x-position |
---|
| 220 | REAL(wp) :: y !< y-position |
---|
| 221 | REAL(wp) :: t !< temperature |
---|
| 222 | REAL(wp) :: t_x !< x-position |
---|
| 223 | REAL(wp) :: t_y !< y-position |
---|
| 224 | REAL(wp), DIMENSION(0:15) :: path_x !< agent path to target (x) |
---|
| 225 | REAL(wp), DIMENSION(0:15) :: path_y !< agent path to target (y) |
---|
[3159] | 226 | END TYPE agent_type |
---|
| 227 | |
---|
[3235] | 228 | TYPE(agent_type), DIMENSION(:), POINTER :: agents !< Agent array for this grid cell |
---|
| 229 | TYPE(agent_type) :: zero_agent !< zero agent to avoid weird thing |
---|
| 230 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: move_also_north !< for agent exchange between PEs |
---|
| 231 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: move_also_south !< for agent exchange between PEs |
---|
| 232 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: agt_gh_l !< ghost layer left of pe domain |
---|
| 233 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: agt_gh_n !< ghost layer north of pe domain |
---|
| 234 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: agt_gh_r !< ghost layer right of pe domain |
---|
| 235 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: agt_gh_s !< ghost layer south of pe domain |
---|
[3159] | 236 | ! |
---|
| 237 | !-- Type for 2D grid on which agents are stored |
---|
| 238 | TYPE grid_agent_def |
---|
| 239 | INTEGER(iwp), DIMENSION(0:3) :: start_index !< start agent index for current block |
---|
| 240 | INTEGER(iwp), DIMENSION(0:3) :: end_index !< end agent index for current block |
---|
| 241 | INTEGER(iwp) :: id_counter !< agent id counter (removeable?) |
---|
| 242 | LOGICAL :: time_loop_done !< timestep loop for agent advection |
---|
[3235] | 243 | TYPE(agent_type), POINTER, DIMENSION(:) :: agents !< Particle array for this grid cell |
---|
[3159] | 244 | END TYPE grid_agent_def |
---|
| 245 | |
---|
| 246 | TYPE(grid_agent_def), DIMENSION(:,:), ALLOCATABLE, TARGET :: grid_agents !< 2D grid on which agents are stored |
---|
| 247 | ! |
---|
| 248 | !-- Item in a priority queue (binary heap) |
---|
| 249 | TYPE heap_item |
---|
| 250 | INTEGER(iwp) :: mesh_id !< id of the submitted mesh point |
---|
| 251 | REAL(wp) :: priority !< priority of the mesh point (= distance so far + heuristic to goal) |
---|
| 252 | END TYPE heap_item |
---|
| 253 | |
---|
| 254 | TYPE(heap_item), DIMENSION(:), ALLOCATABLE :: queue !< priority queue realized as binary heap |
---|
| 255 | ! |
---|
| 256 | !-- Type for mesh point in visibility graph |
---|
| 257 | TYPE mesh_point |
---|
| 258 | INTEGER(iwp) :: polygon_id !< Polygon the point belongs to |
---|
| 259 | INTEGER(iwp) :: vertex_id !< Vertex in the polygon |
---|
| 260 | INTEGER(iwp) :: noc !< number of connections |
---|
| 261 | INTEGER(iwp) :: origin_id !< ID of previous mesh point on path (A*) |
---|
| 262 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: connected_vertices !< Index of connected vertices |
---|
| 263 | REAL(wp) :: cost_so_far !< Cost to reach this mesh point (A*) |
---|
| 264 | REAL(wp) :: x !< x-coordinate |
---|
| 265 | REAL(wp) :: y !< y-coordinate |
---|
| 266 | REAL(wp) :: x_s !< corner shifted outward from building by 1m (x) |
---|
| 267 | REAL(wp) :: y_s !< corner shifted outward from building by 1m (y) |
---|
| 268 | REAL(wp), DIMENSION(:), ALLOCATABLE :: distance_to_vertex !< Distance to each vertex |
---|
| 269 | END TYPE mesh_point |
---|
| 270 | |
---|
| 271 | TYPE(mesh_point), DIMENSION(:), ALLOCATABLE :: mesh !< navigation mesh |
---|
| 272 | TYPE(mesh_point), DIMENSION(:), ALLOCATABLE :: tmp_mesh !< temporary navigation mesh |
---|
| 273 | ! |
---|
| 274 | !-- Vertex of a polygon |
---|
| 275 | TYPE vertex_type |
---|
| 276 | LOGICAL :: delete !< Flag to mark vertex for deletion |
---|
| 277 | REAL(wp) :: x !< x-coordinate |
---|
| 278 | REAL(wp) :: y !< y-coordinate |
---|
| 279 | END TYPE vertex_type |
---|
| 280 | ! |
---|
| 281 | !-- Polygon containing a number of vertices |
---|
| 282 | TYPE polygon_type |
---|
| 283 | INTEGER(iwp) :: nov !< Number of vertices in this polygon |
---|
| 284 | TYPE(vertex_type), DIMENSION(:), ALLOCATABLE :: vertices !< Array of vertices |
---|
| 285 | END TYPE polygon_type |
---|
| 286 | |
---|
| 287 | TYPE(polygon_type), DIMENSION(:), ALLOCATABLE :: polygons !< Building data in polygon form |
---|
| 288 | |
---|
| 289 | SAVE |
---|
| 290 | |
---|
| 291 | PRIVATE |
---|
| 292 | ! |
---|
| 293 | !-- Public functions |
---|
[3235] | 294 | PUBLIC mas_init, mas_last_actions, mas_parin, multi_agent_system |
---|
[3159] | 295 | |
---|
| 296 | ! |
---|
| 297 | !-- Public parameters, constants and initial values |
---|
[3198] | 298 | PUBLIC agents_active |
---|
[3159] | 299 | |
---|
| 300 | INTERFACE mas_parin |
---|
| 301 | MODULE PROCEDURE mas_parin |
---|
| 302 | END INTERFACE mas_parin |
---|
| 303 | |
---|
| 304 | INTERFACE mas_init |
---|
| 305 | MODULE PROCEDURE mas_init |
---|
| 306 | END INTERFACE mas_init |
---|
| 307 | |
---|
[3235] | 308 | INTERFACE mas_last_actions |
---|
| 309 | MODULE PROCEDURE mas_last_actions |
---|
| 310 | END INTERFACE mas_last_actions |
---|
| 311 | |
---|
[3159] | 312 | INTERFACE multi_agent_system |
---|
| 313 | MODULE PROCEDURE multi_agent_system |
---|
| 314 | END INTERFACE multi_agent_system |
---|
| 315 | |
---|
| 316 | CONTAINS |
---|
| 317 | |
---|
| 318 | |
---|
| 319 | !------------------------------------------------------------------------------! |
---|
| 320 | ! Description: |
---|
| 321 | ! ------------ |
---|
| 322 | !> Multi Agent System: |
---|
| 323 | !> executes a number of agents sub-timesteps until the model timestep is reached. |
---|
| 324 | !> The agent timestep is usually smaller than the model timestep |
---|
| 325 | !------------------------------------------------------------------------------! |
---|
| 326 | SUBROUTINE multi_agent_system |
---|
| 327 | |
---|
| 328 | ! USE htcm_mod, & |
---|
| 329 | ! ONLY: htcm_dynamic |
---|
| 330 | |
---|
| 331 | IMPLICIT NONE |
---|
| 332 | |
---|
| 333 | INTEGER(iwp) :: i !< counter |
---|
| 334 | INTEGER(iwp) :: ie !< counter |
---|
| 335 | INTEGER(iwp) :: is !< counter |
---|
| 336 | INTEGER(iwp) :: j !< counter |
---|
| 337 | INTEGER(iwp) :: je !< counter |
---|
| 338 | INTEGER(iwp) :: js !< counter |
---|
| 339 | INTEGER(iwp), SAVE :: mas_count = 0 !< counts the mas-calls |
---|
| 340 | |
---|
[3235] | 341 | LOGICAL :: first_loop_stride !< flag for first loop stride of agent sub-timesteps |
---|
| 342 | LOGICAL, SAVE :: first_call = .TRUE. !< first call of mas flag for output |
---|
[3159] | 343 | |
---|
| 344 | CALL cpu_log( log_point(9), 'mas', 'start' ) |
---|
| 345 | ! |
---|
| 346 | !-- Initialize variables for the next (sub-) timestep, i.e., for marking |
---|
| 347 | !-- those agents to be deleted after the timestep |
---|
| 348 | deleted_agents = 0 |
---|
| 349 | ! |
---|
| 350 | !-- If necessary, release new set of agents |
---|
[3198] | 351 | IF ( time_arel >= dt_arel .AND. end_time_arel > time_since_reference_point ) THEN |
---|
[3159] | 352 | |
---|
| 353 | CALL mas_create_agent(PHASE_RELEASE) |
---|
| 354 | ! |
---|
| 355 | !-- The MOD function allows for changes in the output interval with |
---|
| 356 | !-- restart runs. |
---|
| 357 | time_arel = MOD( time_arel, MAX( dt_arel, dt_3d ) ) |
---|
| 358 | |
---|
| 359 | ENDIF |
---|
| 360 | |
---|
| 361 | first_loop_stride = .TRUE. |
---|
| 362 | grid_agents(:,:)%time_loop_done = .TRUE. |
---|
| 363 | ! |
---|
| 364 | !-- Set timestep variable |
---|
| 365 | IF ( .NOT. agent_own_timestep ) dt_agent = dt_3d |
---|
| 366 | ! |
---|
| 367 | !-- Timestep loop for agent transport. |
---|
| 368 | !-- This loop has to be repeated until the transport time of every agent |
---|
| 369 | !-- (within the total domain!) has reached the LES timestep (dt_3d). |
---|
| 370 | !-- Timestep scheme is Euler-forward |
---|
| 371 | DO |
---|
| 372 | ! |
---|
| 373 | !-- Write agent data at current time on file. |
---|
| 374 | time_write_agent_data = time_write_agent_data + dt_agent |
---|
| 375 | IF ( time_write_agent_data >= dt_write_agent_data ) THEN |
---|
| 376 | #if defined( __netcdf ) |
---|
| 377 | IF ( first_loop_stride ) CALL mas_get_prognostic_quantities |
---|
| 378 | CALL mas_data_output_agents ( first_call ) |
---|
| 379 | #else |
---|
| 380 | WRITE( message_string, * ) 'NetCDF is needed for agent output. ', & |
---|
| 381 | 'Set __netcdf in compiler options' |
---|
| 382 | CALL message( 'multi_agent_system', 'PA0071', 1, 2, 0, 6, 0 ) |
---|
| 383 | #endif |
---|
| 384 | IF(first_call) first_call = .FALSE. |
---|
| 385 | time_write_agent_data = time_write_agent_data - dt_write_agent_data |
---|
| 386 | ENDIF |
---|
| 387 | ! |
---|
| 388 | !-- Flag is true by default, will be set to false if an agent has not yet |
---|
| 389 | !-- reached the model timestep |
---|
| 390 | grid_agents(:,:)%time_loop_done = .TRUE. |
---|
| 391 | |
---|
| 392 | ! |
---|
| 393 | !-- First part of agent transport: |
---|
| 394 | !-- Evaluate social forces for all agents at current positions |
---|
[3235] | 395 | CALL cpu_log( log_point_s(9), 'mas_social_forces', 'start' ) |
---|
[3159] | 396 | DO i = nxl, nxr |
---|
| 397 | DO j = nys, nyn |
---|
| 398 | |
---|
| 399 | number_of_agents = agt_count(j,i) |
---|
| 400 | ! |
---|
| 401 | !-- If grid cell is empty, cycle |
---|
| 402 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 403 | |
---|
| 404 | agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 405 | ! |
---|
| 406 | !-- Evaluation of social forces |
---|
| 407 | CALL mas_timestep_forces_call(i,j) |
---|
| 408 | |
---|
| 409 | ENDDO |
---|
| 410 | ENDDO |
---|
[3235] | 411 | CALL cpu_log( log_point_s(9), 'mas_social_forces', 'stop' ) |
---|
[3159] | 412 | ! |
---|
| 413 | !-- Second part of agent transport: |
---|
| 414 | !-- timestep |
---|
[3235] | 415 | CALL cpu_log( log_point_s(16), 'mas_timestep', 'start' ) |
---|
[3159] | 416 | DO i = nxl, nxr |
---|
| 417 | DO j = nys, nyn |
---|
| 418 | |
---|
| 419 | number_of_agents = agt_count(j,i) |
---|
| 420 | ! |
---|
| 421 | !-- If grid cell is empty, flag must be true |
---|
| 422 | IF ( number_of_agents <= 0 ) THEN |
---|
| 423 | grid_agents(j,i)%time_loop_done = .TRUE. |
---|
| 424 | CYCLE |
---|
| 425 | ENDIF |
---|
| 426 | |
---|
| 427 | agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 428 | |
---|
| 429 | agents(1:number_of_agents)%agent_mask = .TRUE. |
---|
| 430 | ! |
---|
| 431 | !-- Initialize the variable storing the total time that an agent |
---|
| 432 | !-- has advanced within the timestep procedure |
---|
| 433 | IF ( first_loop_stride ) THEN |
---|
| 434 | agents(1:number_of_agents)%dt_sum = 0.0_wp |
---|
| 435 | ENDIF |
---|
| 436 | ! |
---|
| 437 | !-- Initialize the switch used for the loop exit condition checked |
---|
| 438 | !-- at the end of this loop. If at least one agent has failed to |
---|
| 439 | !-- reach the LES timestep, this switch will be set false in |
---|
| 440 | !-- mas_transport. |
---|
| 441 | dt_3d_reached_l_mas = .TRUE. |
---|
| 442 | ! |
---|
| 443 | !-- Timestep |
---|
| 444 | CALL mas_timestep |
---|
| 445 | ! |
---|
| 446 | !-- Delete agents that have been simulated longer than allowed |
---|
| 447 | CALL mas_boundary_conds( 'max_sim_time' ) |
---|
| 448 | ! |
---|
| 449 | !-- Delete agents that have reached target area |
---|
| 450 | CALL mas_boundary_conds( 'target_area' ) |
---|
| 451 | ! |
---|
| 452 | !--- If not all agents of the actual grid cell have reached the |
---|
| 453 | !-- LES timestep, this cell has to to another loop iteration. Due to |
---|
| 454 | !-- the fact that agents can move into neighboring grid cell, |
---|
| 455 | !-- these neighbor cells also have to perform another loop iteration |
---|
| 456 | IF ( .NOT. dt_3d_reached_l_mas ) THEN |
---|
| 457 | js = MAX(nys,j-1) |
---|
| 458 | je = MIN(nyn,j+1) |
---|
| 459 | is = MAX(nxl,i-1) |
---|
| 460 | ie = MIN(nxr,i+1) |
---|
| 461 | grid_agents(js:je,is:ie)%time_loop_done = .FALSE. |
---|
| 462 | ENDIF |
---|
| 463 | |
---|
| 464 | ENDDO |
---|
| 465 | ENDDO |
---|
[3235] | 466 | CALL cpu_log( log_point_s(16), 'mas_timestep', 'stop' ) |
---|
[3159] | 467 | |
---|
| 468 | ! |
---|
| 469 | !-- Find out, if all agents on every PE have completed the LES timestep |
---|
| 470 | !-- and set the switch corespondingly |
---|
| 471 | dt_3d_reached_l_mas = ALL(grid_agents(:,:)%time_loop_done) |
---|
| 472 | #if defined( __parallel ) |
---|
| 473 | IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 474 | CALL MPI_ALLREDUCE( dt_3d_reached_l_mas, dt_3d_reached_mas, 1, MPI_LOGICAL, & |
---|
| 475 | MPI_LAND, comm2d, ierr ) |
---|
| 476 | #else |
---|
| 477 | dt_3d_reached_mas = dt_3d_reached_l_mas |
---|
| 478 | #endif |
---|
| 479 | |
---|
| 480 | ! |
---|
| 481 | !-- Increment time since last release |
---|
| 482 | IF ( dt_3d_reached_mas ) time_arel = time_arel + dt_3d |
---|
| 483 | |
---|
| 484 | ! |
---|
| 485 | !-- Move Agents local to PE to a different grid cell |
---|
[3235] | 486 | CALL cpu_log( log_point_s(18), 'mas_move_exch_sort', 'start' ) |
---|
[3159] | 487 | CALL mas_eh_move_agent |
---|
| 488 | ! |
---|
| 489 | !-- Horizontal boundary conditions including exchange between subdmains |
---|
| 490 | CALL mas_eh_exchange_horiz |
---|
| 491 | ! |
---|
| 492 | !-- Pack agents (eliminate those marked for deletion), |
---|
| 493 | !-- determine new number of agents |
---|
| 494 | CALL mas_ps_sort_in_subboxes |
---|
[3235] | 495 | CALL cpu_log( log_point_s(18), 'mas_move_exch_sort', 'stop' ) |
---|
[3159] | 496 | ! |
---|
| 497 | !-- Initialize variables for the next (sub-) timestep, i.e., for marking |
---|
| 498 | !-- those agents to be deleted after the timestep |
---|
| 499 | deleted_agents = 0 |
---|
| 500 | ! |
---|
| 501 | !-- to be included here: call of human thermal comfort mod (and UV exposure) |
---|
| 502 | ! DO i = nxl, nxr |
---|
| 503 | ! DO j = nys, nyn |
---|
| 504 | ! |
---|
| 505 | ! number_of_agents = agt_count(j,i) |
---|
| 506 | ! ! |
---|
| 507 | ! !-- If grid cell gets empty, cycle |
---|
| 508 | ! IF ( number_of_agents <= 0 ) CYCLE |
---|
| 509 | ! |
---|
| 510 | ! agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 511 | ! ! |
---|
| 512 | ! !-- Evaluation of social forces |
---|
| 513 | ! CALL htcm_dynamic(i,j) |
---|
| 514 | ! |
---|
| 515 | ! ENDDO |
---|
| 516 | ! ENDDO |
---|
| 517 | |
---|
| 518 | IF ( dt_3d_reached_mas ) EXIT |
---|
| 519 | |
---|
| 520 | first_loop_stride = .FALSE. |
---|
| 521 | ENDDO ! timestep loop |
---|
| 522 | |
---|
| 523 | ! |
---|
| 524 | !-- Deallocate unused memory |
---|
| 525 | IF ( deallocate_memory_mas .AND. mas_count == step_dealloc_mas ) THEN |
---|
| 526 | CALL mas_eh_dealloc_agents_array |
---|
| 527 | mas_count = 0 |
---|
| 528 | ELSEIF ( deallocate_memory_mas ) THEN |
---|
| 529 | mas_count = mas_count + 1 |
---|
| 530 | ENDIF |
---|
| 531 | |
---|
| 532 | CALL cpu_log( log_point(9), 'mas', 'stop' ) |
---|
| 533 | |
---|
| 534 | END SUBROUTINE multi_agent_system |
---|
| 535 | |
---|
| 536 | !------------------------------------------------------------------------------! |
---|
| 537 | ! Description: |
---|
| 538 | ! ------------ |
---|
| 539 | !> Calculation of the direction vector from each agent to its current |
---|
| 540 | !> intermittent target |
---|
| 541 | !------------------------------------------------------------------------------! |
---|
| 542 | SUBROUTINE mas_agent_direction |
---|
| 543 | |
---|
| 544 | IMPLICIT NONE |
---|
| 545 | |
---|
| 546 | LOGICAL :: path_flag !< true if new path must be calculated |
---|
| 547 | |
---|
| 548 | INTEGER(iwp) :: n !< loop variable over all agents in a grid box |
---|
| 549 | INTEGER(iwp) :: pc !< agent path counter |
---|
| 550 | |
---|
| 551 | REAL(wp) :: abs_dir !< length of direction vector (for normalization) |
---|
[3241] | 552 | ! REAL(wp) :: d_curr_target !< rounding influence expressed as x speed component |
---|
| 553 | ! REAL(wp) :: d_prev_target !< rounding influence expressed as x speed component |
---|
[3159] | 554 | REAL(wp) :: dir_x !< direction of agent (x) |
---|
| 555 | REAL(wp) :: dir_y !< direction of agent (y) |
---|
[3241] | 556 | ! REAL(wp) :: dist_round = 3. !< distance at which agents start rounding a corner |
---|
[3159] | 557 | REAL(wp) :: dtit !< distance to intermittent target |
---|
[3241] | 558 | ! REAL(wp) :: round_fac = 0.2 !< factor for rounding influence |
---|
| 559 | ! REAL(wp) :: speed_round_x !< rounding influence expressed as x speed component |
---|
| 560 | ! REAL(wp) :: speed_round_y !< rounding influence expressed as x speed component |
---|
[3159] | 561 | |
---|
| 562 | ! |
---|
| 563 | !-- loop over all agents in the current grid box |
---|
| 564 | DO n = 1, number_of_agents |
---|
| 565 | path_flag = .FALSE. |
---|
| 566 | pc = agents(n)%path_counter |
---|
| 567 | ! |
---|
| 568 | !-- If no path was calculated for agent yet, do it |
---|
| 569 | IF ( pc >= 999 ) THEN |
---|
| 570 | CALL mas_nav_find_path(n) |
---|
| 571 | pc = agents(n)%path_counter |
---|
| 572 | ! |
---|
| 573 | !-- Check if new path must be calculated and if so, do it |
---|
| 574 | ELSE |
---|
| 575 | ! |
---|
| 576 | !-- Case one: Agent has come close enough to intermittent target. |
---|
| 577 | !-- -> chose new int target and calculate rest of path if no |
---|
| 578 | !-- new intermittent targets are left |
---|
| 579 | dtit = SQRT((agents(n)%x - agents(n)%path_x(pc))**2 & |
---|
| 580 | + (agents(n)%y - agents(n)%path_y(pc))**2) |
---|
| 581 | IF ( dtit < dist_to_int_target ) THEN |
---|
| 582 | agents(n)%path_counter = agents(n)%path_counter + 1 |
---|
| 583 | pc = agents(n)%path_counter |
---|
| 584 | ! |
---|
| 585 | !-- Path counter out of scope (each agent can store a maximum of 15 |
---|
| 586 | !-- intermittent targets on the way to her final target); new path |
---|
| 587 | !-- must be calculated |
---|
| 588 | IF ( pc >= SIZE(agents(n)%path_x) ) THEN |
---|
| 589 | path_flag = .TRUE. |
---|
| 590 | ENDIF |
---|
| 591 | ! |
---|
| 592 | !-- Case two: Agent too far from path |
---|
| 593 | !-- -> set flag for new path to be calculated |
---|
| 594 | ELSEIF ( dist_point_to_edge(agents(n)%path_x(pc-1), & |
---|
| 595 | agents(n)%path_y(pc-1), & |
---|
| 596 | agents(n)%path_x(pc), & |
---|
| 597 | agents(n)%path_y(pc), & |
---|
| 598 | agents(n)%x, agents(n)%y) & |
---|
| 599 | > max_dist_from_path ) & |
---|
| 600 | THEN |
---|
| 601 | path_flag = .TRUE. |
---|
| 602 | ENDIF |
---|
| 603 | ! |
---|
| 604 | !-- If either of the above two cases was true, calculate new path and |
---|
| 605 | !-- reset 0th path point. This point (the last target the agent had) |
---|
| 606 | !-- is needed for the agents rounding of corners and the calculation |
---|
| 607 | !-- of her deviation from her current path |
---|
| 608 | IF ( path_flag ) THEN |
---|
| 609 | CALL mas_nav_find_path(n) |
---|
| 610 | pc = agents(n)%path_counter |
---|
| 611 | ENDIF |
---|
| 612 | ENDIF |
---|
| 613 | ! |
---|
| 614 | !-- Normalize direction vector |
---|
| 615 | abs_dir = 1.0d-12 |
---|
| 616 | dir_x = agents(n)%path_x(pc) - agents(n)%x |
---|
| 617 | dir_y = agents(n)%path_y(pc) - agents(n)%y |
---|
| 618 | abs_dir = SQRT(dir_x**2 + dir_y**2)+1.0d-12 |
---|
| 619 | !-- needed later for corner rounding |
---|
| 620 | ! dir_x = dir_x/abs_dir |
---|
| 621 | ! dir_y = dir_y/abs_dir |
---|
| 622 | ! dir_x = dir_x + speed_round_x |
---|
| 623 | ! dir_y = dir_y + speed_round_y |
---|
| 624 | ! abs_dir = SQRT(dir_x**2 + dir_y**2)+1.0d-12 |
---|
| 625 | agents(n)%speed_e_x = dir_x/abs_dir |
---|
| 626 | agents(n)%speed_e_y = dir_y/abs_dir |
---|
| 627 | ENDDO |
---|
| 628 | |
---|
| 629 | ! |
---|
| 630 | !-- corner rounding; to be added |
---|
| 631 | ! |
---|
| 632 | !-- Calculate direction change due to rounding of corners |
---|
| 633 | |
---|
| 634 | ! speed_round_x = 0. |
---|
| 635 | ! speed_round_y = 0. |
---|
| 636 | ! |
---|
| 637 | ! d_curr_target = SQRT( (agents(n)%path_x(pc) - agents(n)%x)**2 + & |
---|
| 638 | ! (agents(n)%path_y(pc) - agents(n)%y)**2 ) |
---|
| 639 | ! d_prev_target = SQRT( (agents(n)%path_x(pc-1) - agents(n)%x)**2 + & |
---|
| 640 | ! (agents(n)%path_y(pc-1) - agents(n)%y)**2 ) |
---|
| 641 | ! ! |
---|
| 642 | ! !-- Agent is close to next target and that target is not the final one |
---|
| 643 | ! IF ( d_curr_target < dist_round .AND. dist_round < & |
---|
| 644 | ! SQRT( (agents(n)%path_x(pc) - agents(n)%t_x)**2 + & |
---|
| 645 | ! (agents(n)%path_y(pc) - agents(n)%t_y)**2 ) ) & |
---|
| 646 | ! THEN |
---|
| 647 | ! speed_round_x = (agents(n)%path_x(pc+1) - agents(n)%path_x(pc)) / & |
---|
| 648 | ! ABS( agents(n)%path_x(pc) & |
---|
| 649 | ! - agents(n)%path_x(pc+1)) * round_fac * & |
---|
| 650 | ! SIN( pi/dist_round*d_curr_target ) |
---|
| 651 | ! speed_round_y = (agents(n)%path_y(pc+1) - agents(n)%path_y(pc)) / & |
---|
| 652 | ! ABS( agents(n)%path_y(pc) & |
---|
| 653 | ! - agents(n)%path_y(pc+1)) * round_fac * & |
---|
| 654 | ! SIN( pi/dist_round*d_curr_target ) |
---|
| 655 | ! ENDIF |
---|
| 656 | ! |
---|
| 657 | ! IF ( d_prev_target < dist_round ) THEN |
---|
| 658 | ! IF ( agents(n)%path_x(pc) /= agents(n)%path_x(pc+1) ) THEN |
---|
| 659 | ! speed_round_x = speed_round_x + & |
---|
| 660 | ! (agents(n)%path_x(pc) - agents(n)%path_x(pc+1)) / & |
---|
| 661 | ! ABS( agents(n)%path_x(pc) & |
---|
| 662 | ! - agents(n)%path_x(pc+1)) * round_fac * & |
---|
| 663 | ! SIN( pi/dist_round*d_prev_target ) |
---|
| 664 | ! ENDIF |
---|
| 665 | ! |
---|
| 666 | ! IF ( agents(n)%path_y(pc) /= agents(n)%path_y(pc+1) ) THEN |
---|
| 667 | ! speed_round_y = speed_round_y + & |
---|
| 668 | ! (agents(n)%path_y(pc) - agents(n)%path_y(pc+1)) / & |
---|
| 669 | ! ABS( agents(n)%path_y(pc) & |
---|
| 670 | ! - agents(n)%path_y(pc+1)) * round_fac * & |
---|
| 671 | ! SIN( pi/dist_round*d_prev_target ) |
---|
| 672 | ! ENDIF |
---|
| 673 | |
---|
| 674 | ! ENDIF |
---|
| 675 | |
---|
| 676 | |
---|
| 677 | END SUBROUTINE mas_agent_direction |
---|
| 678 | |
---|
| 679 | !------------------------------------------------------------------------------! |
---|
| 680 | ! Description: |
---|
| 681 | ! ------------ |
---|
| 682 | !> Boundary conditions for maximum time, target reached and out of domain |
---|
| 683 | !------------------------------------------------------------------------------! |
---|
| 684 | SUBROUTINE mas_boundary_conds( location ) |
---|
| 685 | |
---|
| 686 | IMPLICIT NONE |
---|
| 687 | |
---|
| 688 | CHARACTER (LEN=*) :: location !< Identifier |
---|
| 689 | |
---|
| 690 | INTEGER(iwp) :: n !< agent number |
---|
| 691 | INTEGER(iwp) :: grp !< agent group |
---|
| 692 | |
---|
| 693 | REAL(wp) :: dist_to_target !< distance to target |
---|
| 694 | |
---|
| 695 | IF ( location == 'max_sim_time' ) THEN |
---|
| 696 | |
---|
| 697 | ! |
---|
| 698 | !-- Delete agents that have been simulated longer than allowed |
---|
| 699 | DO n = 1, number_of_agents |
---|
| 700 | |
---|
| 701 | IF ( agents(n)%age > agent_maximum_age .AND. & |
---|
| 702 | agents(n)%agent_mask ) & |
---|
| 703 | THEN |
---|
| 704 | agents(n)%agent_mask = .FALSE. |
---|
| 705 | deleted_agents = deleted_agents + 1 |
---|
| 706 | ENDIF |
---|
| 707 | |
---|
| 708 | ENDDO |
---|
| 709 | ENDIF |
---|
| 710 | |
---|
| 711 | IF ( location == 'target_area' ) THEN |
---|
| 712 | |
---|
| 713 | ! |
---|
| 714 | !-- Delete agents that entered target region |
---|
| 715 | DO n = 1, number_of_agents |
---|
| 716 | grp = agents(n)%group |
---|
| 717 | dist_to_target = SQRT((agents(n)%x-at_x(grp))**2 & |
---|
| 718 | + (agents(n)%y-at_y(grp))**2) |
---|
| 719 | IF ( dist_to_target < dist_target_reached ) THEN |
---|
| 720 | agents(n)%agent_mask = .FALSE. |
---|
| 721 | deleted_agents = deleted_agents + 1 |
---|
| 722 | ENDIF |
---|
| 723 | |
---|
| 724 | ENDDO |
---|
| 725 | ENDIF |
---|
| 726 | |
---|
| 727 | END SUBROUTINE mas_boundary_conds |
---|
| 728 | |
---|
| 729 | !------------------------------------------------------------------------------! |
---|
| 730 | ! Description: |
---|
| 731 | ! ------------ |
---|
| 732 | !> Release new agents at their respective sources |
---|
| 733 | !------------------------------------------------------------------------------! |
---|
| 734 | SUBROUTINE mas_create_agent (phase) |
---|
| 735 | |
---|
| 736 | IMPLICIT NONE |
---|
| 737 | |
---|
| 738 | INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for agents |
---|
| 739 | INTEGER(iwp) :: i !< loop variable ( agent groups ) |
---|
| 740 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 741 | INTEGER(iwp) :: jp !< index variable along y |
---|
| 742 | INTEGER(iwp) :: loop_stride !< loop variable for initialization |
---|
| 743 | INTEGER(iwp) :: n !< loop variable ( number of agents ) |
---|
| 744 | INTEGER(iwp) :: new_size !< new size of allocated memory for agents |
---|
| 745 | INTEGER(iwp) :: rn_side !< index of agent path |
---|
| 746 | |
---|
| 747 | INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization |
---|
| 748 | |
---|
| 749 | INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg) :: local_count !< start address of new agent |
---|
| 750 | INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg) :: local_start !< start address of new agent |
---|
| 751 | |
---|
| 752 | LOGICAL :: first_stride !< flag for initialization |
---|
| 753 | |
---|
| 754 | REAL(wp) :: pos_x !< increment for agent position in x |
---|
| 755 | REAL(wp) :: pos_y !< increment for agent position in y |
---|
| 756 | REAL(wp) :: rand_contr !< dummy argument for random position |
---|
| 757 | REAL(wp) :: rn_side_dum !< index of agent path |
---|
| 758 | |
---|
| 759 | TYPE(agent_type),TARGET :: tmp_agent !< temporary agent used for initialization |
---|
| 760 | |
---|
| 761 | ! |
---|
| 762 | !-- Calculate agent positions and store agent attributes, if |
---|
| 763 | !-- agent is situated on this PE |
---|
| 764 | DO loop_stride = 1, 2 |
---|
| 765 | first_stride = (loop_stride == 1) |
---|
| 766 | IF ( first_stride ) THEN |
---|
| 767 | local_count = 0 ! count number of agents |
---|
| 768 | ELSE |
---|
| 769 | local_count = agt_count ! Start address of new agents |
---|
| 770 | ENDIF |
---|
| 771 | |
---|
| 772 | DO i = 1, number_of_agent_groups |
---|
| 773 | |
---|
| 774 | pos_y = ass(i) |
---|
| 775 | |
---|
| 776 | DO WHILE ( pos_y <= asn(i) ) |
---|
| 777 | |
---|
| 778 | IF ( pos_y >= nys * dy .AND. & |
---|
| 779 | pos_y < ( nyn + 1 ) * dy ) & |
---|
| 780 | THEN |
---|
| 781 | |
---|
| 782 | pos_x = asl(i) |
---|
| 783 | |
---|
| 784 | xloop: DO WHILE ( pos_x <= asr(i) ) |
---|
| 785 | |
---|
| 786 | IF ( pos_x >= nxl * dx .AND. & |
---|
| 787 | pos_x < ( nxr + 1) * dx ) & |
---|
| 788 | THEN |
---|
| 789 | |
---|
| 790 | tmp_agent%agent_mask = .TRUE. |
---|
| 791 | tmp_agent%group = i |
---|
| 792 | tmp_agent%id = 0_idp |
---|
| 793 | tmp_agent%block_nr = -1 |
---|
| 794 | tmp_agent%path_counter = 999 !SIZE(tmp_agent%path_x) |
---|
| 795 | tmp_agent%age = 0.0_wp |
---|
| 796 | tmp_agent%age_m = 0.0_wp |
---|
| 797 | tmp_agent%dt_sum = 0.0_wp |
---|
| 798 | tmp_agent%clo = 99999.0_wp |
---|
| 799 | tmp_agent%energy_storage= 0.0_wp |
---|
| 800 | tmp_agent%force_x = 0.0_wp |
---|
| 801 | tmp_agent%force_y = 0.0_wp |
---|
| 802 | tmp_agent%origin_x = pos_x |
---|
| 803 | tmp_agent%origin_y = pos_y |
---|
| 804 | tmp_agent%speed_abs = 0.0_wp |
---|
| 805 | tmp_agent%speed_e_x = 0.0_wp |
---|
| 806 | tmp_agent%speed_e_y = 0.0_wp |
---|
| 807 | tmp_agent%speed_des = random_normal(desired_speed,& |
---|
| 808 | des_sp_sig) |
---|
| 809 | tmp_agent%speed_x = 0.0_wp |
---|
| 810 | tmp_agent%speed_y = 0.0_wp |
---|
| 811 | tmp_agent%thermal_index = 99999.0_wp |
---|
| 812 | tmp_agent%x = pos_x |
---|
| 813 | tmp_agent%y = pos_y |
---|
| 814 | tmp_agent%path_x = -1.0_wp |
---|
| 815 | tmp_agent%path_y = -1.0_wp |
---|
| 816 | tmp_agent%t_x = - pi |
---|
| 817 | tmp_agent%t_y = - pi |
---|
| 818 | ! |
---|
| 819 | !-- Determine the grid indices of the agent position |
---|
| 820 | ip = tmp_agent%x * ddx |
---|
| 821 | jp = tmp_agent%y * ddy |
---|
| 822 | ! |
---|
| 823 | !-- Give each agent its target |
---|
| 824 | IF ( a_rand_target(i) ) THEN |
---|
| 825 | ! |
---|
| 826 | !-- Agent shall receive random target just outside |
---|
| 827 | !-- simulated area |
---|
| 828 | rn_side_dum = random_function(iran_agent) |
---|
| 829 | rn_side = FLOOR(4.*rn_side_dum) |
---|
| 830 | IF ( rn_side < 2 ) THEN |
---|
| 831 | IF ( rn_side == 0 ) THEN |
---|
| 832 | tmp_agent%t_y = -2*dy |
---|
| 833 | ELSE |
---|
| 834 | tmp_agent%t_y = (ny+3)*dy |
---|
| 835 | ENDIF |
---|
| 836 | tmp_agent%t_x = random_function(iran_agent) * & |
---|
| 837 | (nx+1)*dx |
---|
| 838 | ELSE |
---|
| 839 | IF ( rn_side == 2 ) THEN |
---|
| 840 | tmp_agent%t_x = -2*dx |
---|
| 841 | ELSE |
---|
| 842 | tmp_agent%t_x = (nx+3)*dx |
---|
| 843 | ENDIF |
---|
| 844 | tmp_agent%t_y = random_function(iran_agent) * & |
---|
| 845 | (ny+1)*dy |
---|
| 846 | ENDIF |
---|
| 847 | ! |
---|
| 848 | !-- Agent gets target of her group |
---|
| 849 | ELSE |
---|
| 850 | tmp_agent%t_x = at_x(i) |
---|
| 851 | tmp_agent%t_y = at_y(i) |
---|
| 852 | ENDIF |
---|
| 853 | |
---|
| 854 | local_count(jp,ip) = local_count(jp,ip) + 1 |
---|
| 855 | |
---|
| 856 | IF ( .NOT. first_stride ) THEN |
---|
| 857 | grid_agents(jp,ip)%agents(local_count(jp,ip)) & |
---|
| 858 | = tmp_agent |
---|
| 859 | ENDIF |
---|
| 860 | |
---|
| 861 | ENDIF |
---|
| 862 | |
---|
| 863 | pos_x = pos_x + adx(i) |
---|
| 864 | |
---|
| 865 | ENDDO xloop |
---|
| 866 | |
---|
| 867 | ENDIF |
---|
| 868 | |
---|
| 869 | pos_y = pos_y + ady(i) |
---|
| 870 | |
---|
| 871 | ENDDO |
---|
| 872 | |
---|
| 873 | ENDDO |
---|
| 874 | |
---|
| 875 | ! |
---|
| 876 | !-- Allocate or reallocate agents array to new size |
---|
| 877 | IF ( first_stride ) THEN |
---|
| 878 | DO ip = nxlg, nxrg |
---|
| 879 | DO jp = nysg, nyng |
---|
| 880 | IF ( phase == PHASE_INIT ) THEN |
---|
| 881 | IF ( local_count(jp,ip) > 0 ) THEN |
---|
| 882 | alloc_size = MAX( INT( local_count(jp,ip) * & |
---|
| 883 | ( 1.0_wp + alloc_factor_mas / 100.0_wp ) ), & |
---|
| 884 | min_nr_agent ) |
---|
| 885 | ELSE |
---|
| 886 | alloc_size = min_nr_agent |
---|
| 887 | ENDIF |
---|
| 888 | ALLOCATE(grid_agents(jp,ip)%agents(1:alloc_size)) |
---|
| 889 | DO n = 1, alloc_size |
---|
| 890 | grid_agents(jp,ip)%agents(n) = zero_agent |
---|
| 891 | ENDDO |
---|
| 892 | ELSEIF ( phase == PHASE_RELEASE ) THEN |
---|
| 893 | IF ( local_count(jp,ip) > 0 ) THEN |
---|
| 894 | new_size = local_count(jp,ip) + agt_count(jp,ip) |
---|
| 895 | alloc_size = MAX( INT( new_size * ( 1.0_wp + & |
---|
| 896 | alloc_factor_mas / 100.0_wp ) ), min_nr_agent ) |
---|
| 897 | IF( alloc_size > SIZE( grid_agents(jp,ip)%agents) ) & |
---|
| 898 | THEN |
---|
| 899 | CALL mas_eh_realloc_agents_array(ip,jp,alloc_size) |
---|
| 900 | ENDIF |
---|
| 901 | ENDIF |
---|
| 902 | ENDIF |
---|
| 903 | ENDDO |
---|
| 904 | ENDDO |
---|
| 905 | ENDIF |
---|
| 906 | |
---|
| 907 | ENDDO |
---|
| 908 | |
---|
| 909 | local_start = agt_count+1 |
---|
| 910 | agt_count = local_count |
---|
| 911 | |
---|
| 912 | ! |
---|
| 913 | !-- Calculate agent IDs |
---|
| 914 | DO ip = nxl, nxr |
---|
| 915 | DO jp = nys, nyn |
---|
| 916 | number_of_agents = agt_count(jp,ip) |
---|
| 917 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 918 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 919 | |
---|
| 920 | DO n = local_start(jp,ip), number_of_agents !only new agents |
---|
| 921 | |
---|
| 922 | agents(n)%id = 10000_idp**2 * grid_agents(jp,ip)%id_counter + & |
---|
| 923 | 10000_idp * jp + ip |
---|
| 924 | ! |
---|
| 925 | !-- Count the number of agents that have been released before |
---|
| 926 | grid_agents(jp,ip)%id_counter = grid_agents(jp,ip)%id_counter & |
---|
| 927 | + 1 |
---|
| 928 | |
---|
| 929 | ENDDO |
---|
| 930 | |
---|
| 931 | ENDDO |
---|
| 932 | ENDDO |
---|
| 933 | |
---|
| 934 | ! |
---|
| 935 | !-- Add random fluctuation to agent positions. |
---|
| 936 | IF ( random_start_position_agents ) THEN |
---|
| 937 | DO ip = nxl, nxr |
---|
| 938 | DO jp = nys, nyn |
---|
| 939 | number_of_agents = agt_count(jp,ip) |
---|
| 940 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 941 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 942 | ! |
---|
| 943 | !-- Move only new agents. Moreover, limit random fluctuation |
---|
| 944 | !-- in order to prevent that agents move more than one grid box, |
---|
| 945 | !-- which would lead to problems concerning agent exchange |
---|
| 946 | !-- between processors in case adx/ady are larger than dx/dy, |
---|
| 947 | !-- respectively. |
---|
| 948 | DO n = local_start(jp,ip), number_of_agents |
---|
| 949 | IF ( asl(agents(n)%group) /= asr(agents(n)%group) ) THEN |
---|
| 950 | rand_contr = ( random_function( iran_agent ) - 0.5_wp ) *& |
---|
| 951 | adx(agents(n)%group) |
---|
| 952 | agents(n)%x = agents(n)%x + & |
---|
| 953 | MERGE( rand_contr, SIGN( dx, rand_contr ), & |
---|
| 954 | ABS( rand_contr ) < dx & |
---|
| 955 | ) |
---|
| 956 | ENDIF |
---|
| 957 | IF ( ass(agents(n)%group) /= asn(agents(n)%group) ) THEN |
---|
| 958 | rand_contr = ( random_function( iran_agent ) - 0.5_wp ) *& |
---|
| 959 | ady(agents(n)%group) |
---|
| 960 | agents(n)%y = agents(n)%y + & |
---|
| 961 | MERGE( rand_contr, SIGN( dy, rand_contr ), & |
---|
| 962 | ABS( rand_contr ) < dy ) |
---|
| 963 | ENDIF |
---|
| 964 | ENDDO |
---|
| 965 | ! |
---|
| 966 | !-- Delete agents that have been simulated longer than allowed |
---|
| 967 | CALL mas_boundary_conds( 'max_sim_time' ) |
---|
| 968 | ! |
---|
| 969 | !-- Delete agents that have reached target area |
---|
| 970 | CALL mas_boundary_conds( 'target_area' ) |
---|
| 971 | |
---|
| 972 | ENDDO |
---|
| 973 | ENDDO |
---|
| 974 | ! |
---|
| 975 | !-- Exchange agents between grid cells and processors |
---|
| 976 | CALL mas_eh_move_agent |
---|
| 977 | CALL mas_eh_exchange_horiz |
---|
| 978 | |
---|
| 979 | ENDIF |
---|
| 980 | ! |
---|
| 981 | !-- In case of random_start_position_agents, delete agents identified by |
---|
| 982 | !-- mas_eh_exchange_horiz and mas_boundary_conds. Then sort agents into |
---|
| 983 | !-- blocks, which is needed for a fast interpolation of the LES fields |
---|
| 984 | !-- on the agent position. |
---|
| 985 | CALL mas_ps_sort_in_subboxes |
---|
| 986 | |
---|
| 987 | ! |
---|
| 988 | !-- Determine the current number of agents |
---|
| 989 | number_of_agents = 0 |
---|
| 990 | DO ip = nxl, nxr |
---|
| 991 | DO jp = nys, nyn |
---|
| 992 | number_of_agents = number_of_agents + agt_count(jp,ip) |
---|
| 993 | ENDDO |
---|
| 994 | ENDDO |
---|
| 995 | ! |
---|
| 996 | !-- Calculate the number of agents of the total domain |
---|
| 997 | #if defined( __parallel ) |
---|
| 998 | IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 999 | CALL MPI_ALLREDUCE( number_of_agents, total_number_of_agents, 1, & |
---|
| 1000 | MPI_INTEGER, MPI_SUM, comm2d, ierr ) |
---|
| 1001 | #else |
---|
| 1002 | total_number_of_agents = number_of_agents |
---|
| 1003 | #endif |
---|
| 1004 | |
---|
| 1005 | RETURN |
---|
| 1006 | |
---|
| 1007 | END SUBROUTINE mas_create_agent |
---|
| 1008 | |
---|
| 1009 | !------------------------------------------------------------------------------! |
---|
| 1010 | ! Description: |
---|
| 1011 | ! ------------ |
---|
| 1012 | !> Creates flags that indicate if a gridbox contains edges or corners. These |
---|
| 1013 | !> flags are used for agents to check if obstacles are close to them. |
---|
| 1014 | !------------------------------------------------------------------------------! |
---|
| 1015 | SUBROUTINE mas_create_obstacle_flags |
---|
| 1016 | |
---|
| 1017 | USE arrays_3d, & |
---|
| 1018 | ONLY: zw |
---|
| 1019 | |
---|
| 1020 | IMPLICIT NONE |
---|
| 1021 | |
---|
| 1022 | INTEGER(iwp) :: il |
---|
| 1023 | INTEGER(iwp) :: jl |
---|
| 1024 | |
---|
| 1025 | ALLOCATE(obstacle_flags(nysg:nyng,nxlg:nxrg)) |
---|
| 1026 | |
---|
| 1027 | obstacle_flags = 0 |
---|
| 1028 | |
---|
| 1029 | DO il = nxlg, nxrg |
---|
| 1030 | DO jl = nysg, nyng |
---|
| 1031 | ! |
---|
[3187] | 1032 | !-- Exclude cyclic topography boundary |
---|
| 1033 | IF ( il < 0 .OR. il > nx .OR. jl < 0 .OR. jl > ny ) CYCLE |
---|
| 1034 | ! |
---|
[3159] | 1035 | !-- North edge |
---|
| 1036 | IF ( jl < nyng ) THEN |
---|
| 1037 | IF ( ( top_top_s(jl,il) - top_top_s(jl+1,il) ) > 1 .AND. & |
---|
| 1038 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1039 | zw( top_top_w(jl+1,il) ) ) > .51_wp ) & |
---|
| 1040 | THEN |
---|
| 1041 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 0 ) |
---|
| 1042 | ENDIF |
---|
| 1043 | ENDIF |
---|
| 1044 | ! |
---|
| 1045 | !-- North right corner |
---|
| 1046 | IF ( jl < nyng .AND. il < nxrg ) THEN |
---|
| 1047 | IF ( ( top_top_s(jl,il) - top_top_s(jl+1,il) ) > 1 .AND. & |
---|
| 1048 | ( top_top_s(jl,il) - top_top_s(jl+1,il+1) ) > 1 .AND. & |
---|
| 1049 | ( top_top_s(jl,il) - top_top_s(jl,il+1) ) > 1 .AND. & |
---|
| 1050 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1051 | zw( top_top_w(jl+1,il+1) ) ) > .51_wp ) & |
---|
| 1052 | THEN |
---|
| 1053 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 1 ) |
---|
| 1054 | ENDIF |
---|
| 1055 | ENDIF |
---|
| 1056 | ! |
---|
| 1057 | !-- Right edge |
---|
| 1058 | IF ( il < nxrg ) THEN |
---|
| 1059 | IF ( ( top_top_s(jl,il) - top_top_s(jl,il+1) ) > 1 .AND. & |
---|
| 1060 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1061 | zw( top_top_w(jl,il+1) ) ) > .51_wp ) & |
---|
| 1062 | THEN |
---|
| 1063 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 2 ) |
---|
| 1064 | ENDIF |
---|
| 1065 | ENDIF |
---|
| 1066 | ! |
---|
| 1067 | !-- South right corner |
---|
| 1068 | IF ( jl > nysg .AND. il < nxrg ) THEN |
---|
| 1069 | IF ( ( top_top_s(jl,il) - top_top_s(jl,il+1) ) > 1 .AND. & |
---|
| 1070 | ( top_top_s(jl,il) - top_top_s(jl-1,il+1) ) > 1 .AND. & |
---|
| 1071 | ( top_top_s(jl,il) - top_top_s(jl-1,il) ) > 1 .AND. & |
---|
| 1072 | ( zw(top_top_w(jl,il)) - & |
---|
| 1073 | zw( top_top_w(jl-1,il+1) ) ) > .51_wp ) & |
---|
| 1074 | THEN |
---|
| 1075 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 3 ) |
---|
| 1076 | ENDIF |
---|
| 1077 | ENDIF |
---|
| 1078 | ! |
---|
| 1079 | !-- South edge |
---|
| 1080 | IF ( jl > nysg ) THEN |
---|
| 1081 | IF ( ( top_top_s(jl,il) - top_top_s(jl-1,il) ) > 1 .AND. & |
---|
| 1082 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1083 | zw( top_top_w(jl-1,il) ) ) > .51_wp ) & |
---|
| 1084 | THEN |
---|
| 1085 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 4 ) |
---|
| 1086 | ENDIF |
---|
| 1087 | ENDIF |
---|
| 1088 | ! |
---|
| 1089 | !-- South left corner |
---|
| 1090 | IF ( jl > nysg .AND. il > nxlg ) THEN |
---|
| 1091 | IF ( ( top_top_s(jl,il) - top_top_s(jl-1,il) ) > 1 .AND. & |
---|
| 1092 | ( top_top_s(jl,il) - top_top_s(jl-1,il-1) ) > 1 .AND. & |
---|
| 1093 | ( top_top_s(jl,il) - top_top_s(jl,il-1) ) > 1 .AND. & |
---|
| 1094 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1095 | zw(top_top_w(jl-1,il-1) ) ) > .51_wp ) & |
---|
| 1096 | THEN |
---|
| 1097 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 5 ) |
---|
| 1098 | ENDIF |
---|
| 1099 | ENDIF |
---|
| 1100 | ! |
---|
| 1101 | !-- Left edge |
---|
| 1102 | IF ( il > nxlg ) THEN |
---|
| 1103 | IF ( ( top_top_s(jl,il) - top_top_s(jl,il-1) ) > 1 .AND. & |
---|
| 1104 | ( zw(top_top_w(jl,il) ) - & |
---|
| 1105 | zw(top_top_w(jl,il-1) ) ) > .51_wp ) & |
---|
| 1106 | THEN |
---|
| 1107 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 6 ) |
---|
| 1108 | ENDIF |
---|
| 1109 | ENDIF |
---|
| 1110 | ! |
---|
| 1111 | !-- North left corner |
---|
| 1112 | IF ( jl < nyng .AND. il > nxlg ) THEN |
---|
| 1113 | IF ( ( top_top_s(jl,il) - top_top_s(jl,il-1) ) > 1 .AND. & |
---|
| 1114 | ( top_top_s(jl,il) - top_top_s(jl+1,il-1) ) > 1 .AND. & |
---|
| 1115 | ( top_top_s(jl,il) - top_top_s(jl+1,il) ) > 1 .AND. & |
---|
| 1116 | ( zw( top_top_w(jl,il) ) - & |
---|
| 1117 | zw( top_top_w(jl+1,il-1) ) ) > .51_wp ) & |
---|
| 1118 | THEN |
---|
| 1119 | obstacle_flags(jl,il) = IBSET( obstacle_flags(jl,il), 7 ) |
---|
| 1120 | ENDIF |
---|
| 1121 | ENDIF |
---|
| 1122 | |
---|
| 1123 | ENDDO |
---|
| 1124 | ENDDO |
---|
| 1125 | |
---|
| 1126 | END SUBROUTINE mas_create_obstacle_flags |
---|
| 1127 | |
---|
| 1128 | !------------------------------------------------------------------------------! |
---|
| 1129 | ! Description: |
---|
| 1130 | ! ------------ |
---|
| 1131 | !> Write agent data in netCDF format |
---|
| 1132 | !------------------------------------------------------------------------------! |
---|
| 1133 | SUBROUTINE mas_data_output_agents( ftest ) |
---|
| 1134 | |
---|
| 1135 | USE control_parameters, & |
---|
[3235] | 1136 | ONLY: agt_time_count, end_time, message_string, & |
---|
| 1137 | multi_agent_system_end, multi_agent_system_start |
---|
[3159] | 1138 | |
---|
| 1139 | USE netcdf_interface, & |
---|
[3241] | 1140 | ONLY: nc_stat, id_set_agt, id_var_time_agt, & |
---|
[3159] | 1141 | id_var_agt, netcdf_handle_error |
---|
| 1142 | |
---|
| 1143 | USE pegrid |
---|
| 1144 | |
---|
| 1145 | #if defined( __netcdf ) |
---|
| 1146 | USE NETCDF |
---|
| 1147 | #endif |
---|
[3235] | 1148 | USE mas_global_attributes, & |
---|
| 1149 | ONLY: dim_size_agtnum |
---|
[3159] | 1150 | |
---|
| 1151 | IMPLICIT NONE |
---|
| 1152 | |
---|
| 1153 | INTEGER(iwp) :: agt_size !< Agent size in bytes |
---|
| 1154 | INTEGER(iwp) :: dummy !< dummy |
---|
[3235] | 1155 | INTEGER(iwp) :: ii !< counter (x) |
---|
[3159] | 1156 | INTEGER(iwp) :: ip !< counter (x) |
---|
| 1157 | INTEGER(iwp) :: jp !< counter (y) |
---|
| 1158 | INTEGER(iwp) :: n !< counter (number of PEs) |
---|
| 1159 | INTEGER(iwp) :: noa !< number of agents |
---|
| 1160 | INTEGER(iwp) :: noa_rcv !< received number of agents |
---|
| 1161 | INTEGER(iwp) :: out_noa !< number of agents for output |
---|
[3235] | 1162 | |
---|
[3159] | 1163 | INTEGER(iwp), DIMENSION(0:numprocs-1) :: noa_arr !< number of agents on each PE |
---|
| 1164 | |
---|
| 1165 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE, TARGET :: trf_agents !< all agents on current PE |
---|
| 1166 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE, TARGET :: out_agents !< all agents in entire domain |
---|
[3235] | 1167 | |
---|
[3159] | 1168 | LOGICAL, INTENT (INOUT) :: ftest |
---|
| 1169 | |
---|
[3235] | 1170 | LOGICAL, SAVE :: agt_dimension_exceeded = .FALSE. |
---|
| 1171 | |
---|
[3159] | 1172 | CALL cpu_log( log_point_s(17), 'mas_data_output', 'start' ) |
---|
| 1173 | ! |
---|
| 1174 | !-- Get total number of agents and put all agents on one PE in one array |
---|
| 1175 | noa = 0 |
---|
| 1176 | DO ip = nxl, nxr |
---|
| 1177 | DO jp = nys, nyn |
---|
| 1178 | noa = noa + agt_count(jp,ip) |
---|
| 1179 | ENDDO |
---|
| 1180 | ENDDO |
---|
| 1181 | IF(noa > 0) THEN |
---|
| 1182 | ALLOCATE(trf_agents(1:noa)) |
---|
| 1183 | dummy = 1 |
---|
| 1184 | DO ip = nxl, nxr |
---|
| 1185 | DO jp = nys, nyn |
---|
| 1186 | IF ( agt_count(jp,ip) == 0 ) CYCLE |
---|
| 1187 | agents => grid_agents(jp,ip)%agents(1:agt_count(jp,ip)) |
---|
| 1188 | trf_agents(dummy:(dummy-1+agt_count(jp,ip))) = agents |
---|
| 1189 | dummy = dummy + agt_count(jp,ip) |
---|
| 1190 | ENDDO |
---|
| 1191 | ENDDO |
---|
| 1192 | ENDIF |
---|
[3201] | 1193 | #if defined( __parallel ) |
---|
[3159] | 1194 | ! |
---|
| 1195 | !-- Gather all agents on PE0 for output |
---|
| 1196 | IF ( myid == 0 ) THEN |
---|
| 1197 | noa_arr(0) = noa |
---|
| 1198 | ! |
---|
| 1199 | !-- Receive data from all other PEs. |
---|
| 1200 | DO n = 1, numprocs-1 |
---|
| 1201 | CALL MPI_RECV( noa_arr(n), 1, MPI_INTEGER, & |
---|
| 1202 | n, 0, comm2d, status, ierr ) |
---|
| 1203 | ENDDO |
---|
| 1204 | ELSE |
---|
| 1205 | CALL MPI_SEND( noa, 1, MPI_INTEGER, 0, 0, comm2d, ierr ) |
---|
| 1206 | ENDIF |
---|
| 1207 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
[3160] | 1208 | agt_size = STORAGE_SIZE(zero_agent)/8 |
---|
[3159] | 1209 | IF ( myid == 0 ) THEN |
---|
| 1210 | ! |
---|
| 1211 | !-- Receive data from all other PEs. |
---|
| 1212 | out_noa = SUM(noa_arr) |
---|
| 1213 | IF ( out_noa > 0 ) THEN |
---|
| 1214 | ALLOCATE( out_agents(1:out_noa) ) |
---|
| 1215 | IF ( noa > 0 ) THEN |
---|
| 1216 | out_agents(1:noa) = trf_agents |
---|
| 1217 | ENDIF |
---|
| 1218 | noa_rcv = noa |
---|
| 1219 | DO n = 1, numprocs-1 |
---|
| 1220 | IF ( noa_arr(n) > 0 ) THEN |
---|
| 1221 | CALL MPI_RECV( out_agents(noa_rcv+1), noa_arr(n)*agt_size, & |
---|
| 1222 | MPI_BYTE, n, 0, comm2d, status, ierr ) |
---|
| 1223 | noa_rcv = noa_rcv + noa_arr(n) |
---|
| 1224 | ENDIF |
---|
| 1225 | ENDDO |
---|
| 1226 | ELSE |
---|
| 1227 | ALLOCATE( out_agents(1:2) ) |
---|
| 1228 | out_agents = zero_agent |
---|
| 1229 | out_noa = 2 |
---|
| 1230 | ENDIF |
---|
| 1231 | ELSE |
---|
| 1232 | IF ( noa > 0 ) THEN |
---|
| 1233 | CALL MPI_SEND( trf_agents(1), noa*agt_size, MPI_BYTE, 0, 0, & |
---|
| 1234 | comm2d, ierr ) |
---|
| 1235 | ENDIF |
---|
| 1236 | ENDIF |
---|
| 1237 | ! |
---|
| 1238 | !-- A barrier has to be set, because otherwise some PEs may |
---|
| 1239 | !-- proceed too fast so that PE0 may receive wrong data on |
---|
| 1240 | !-- tag 0 |
---|
| 1241 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
[3201] | 1242 | #endif |
---|
[3159] | 1243 | IF ( myid == 0 ) THEN |
---|
[3201] | 1244 | #if defined( __parallel ) |
---|
[3241] | 1245 | agents => out_agents |
---|
[3201] | 1246 | #else |
---|
[3241] | 1247 | agents => trf_agents |
---|
[3201] | 1248 | #endif |
---|
[3159] | 1249 | |
---|
| 1250 | #if defined( __netcdf ) |
---|
[3235] | 1251 | ! |
---|
| 1252 | !-- Update maximum number of agents |
---|
| 1253 | maximum_number_of_agents = MAX(maximum_number_of_agents, out_noa) |
---|
| 1254 | ! |
---|
| 1255 | !-- Output in netCDF format |
---|
| 1256 | IF ( ftest ) THEN |
---|
| 1257 | ! |
---|
| 1258 | !-- First, define size of agent number dimension from amount of agents |
---|
| 1259 | !-- released, release interval, time of agent simulation and max |
---|
| 1260 | !-- age of agents |
---|
| 1261 | dim_size_agtnum = MIN( MIN( multi_agent_system_end, end_time ) & |
---|
| 1262 | - multi_agent_system_start, & |
---|
| 1263 | agent_maximum_age) |
---|
[3159] | 1264 | |
---|
[3235] | 1265 | DO ii = 1, number_of_agent_groups |
---|
| 1266 | dim_size_agtnum = dim_size_agtnum & |
---|
| 1267 | + (FLOOR( ( asr(ii)-asl(ii) ) / adx(ii) ) + 1) & |
---|
| 1268 | * (FLOOR( ( asn(ii)-ass(ii) ) / ady(ii) ) + 1) & |
---|
| 1269 | * (FLOOR( dim_size_agtnum / dt_arel ) + 1) & |
---|
| 1270 | * dim_size_factor_agtnum |
---|
| 1271 | dim_size_agtnum = MIN( dim_size_agtnum, dim_size_agtnum_manual ) |
---|
| 1272 | ENDDO |
---|
| 1273 | CALL check_open( 118 ) |
---|
| 1274 | ENDIF |
---|
[3159] | 1275 | |
---|
[3235] | 1276 | ! |
---|
| 1277 | !-- Update the NetCDF time axis |
---|
[3159] | 1278 | agt_time_count = agt_time_count + 1 |
---|
| 1279 | |
---|
[3235] | 1280 | IF ( .NOT. agt_dimension_exceeded ) THEN |
---|
| 1281 | ! |
---|
| 1282 | !-- if number of agents to be output exceeds dimension, set flag and |
---|
| 1283 | !-- print warning |
---|
| 1284 | IF ( out_noa > dim_size_agtnum ) THEN |
---|
| 1285 | |
---|
| 1286 | agt_dimension_exceeded = .TRUE. |
---|
| 1287 | WRITE(message_string,'(A,F11.1,2(A,I8))') & |
---|
| 1288 | 'Number of agents exceeds agent dimension.' // & |
---|
| 1289 | '&Starting at time_since_reference_point = ', & |
---|
| 1290 | time_since_reference_point, & |
---|
| 1291 | ' s, &data may be missing.'// & |
---|
| 1292 | '&Number of agents: ', out_noa, & |
---|
| 1293 | '&Agent dimension size: ', dim_size_agtnum |
---|
| 1294 | |
---|
| 1295 | CALL message( 'mas_data_output_agents', & |
---|
| 1296 | 'PA0420', 0, 1, 0, 6, 0 ) |
---|
| 1297 | |
---|
| 1298 | ENDIF |
---|
| 1299 | ENDIF |
---|
| 1300 | |
---|
| 1301 | ! |
---|
| 1302 | !-- reduce number of output agents to dimension size, if necessary |
---|
| 1303 | IF ( agt_dimension_exceeded ) THEN |
---|
| 1304 | |
---|
| 1305 | out_noa = MIN( out_noa, dim_size_agtnum ) |
---|
| 1306 | |
---|
| 1307 | ENDIF |
---|
| 1308 | |
---|
| 1309 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_time_agt, & |
---|
| 1310 | (/ time_since_reference_point /), & |
---|
| 1311 | start = (/ agt_time_count /), & |
---|
| 1312 | count = (/ 1 /) ) |
---|
[3159] | 1313 | CALL netcdf_handle_error( 'mas_data_output_agents', 1 ) |
---|
| 1314 | |
---|
[3235] | 1315 | ! |
---|
| 1316 | !-- Output agent attributes |
---|
| 1317 | |
---|
| 1318 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(1), agents%id, & |
---|
[3159] | 1319 | start = (/ 1, agt_time_count /), & |
---|
| 1320 | count = (/ out_noa /) ) |
---|
[3235] | 1321 | CALL netcdf_handle_error( 'mas_data_output_agents', 2 ) |
---|
| 1322 | |
---|
| 1323 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(2), agents%x, & |
---|
| 1324 | start = (/ 1, agt_time_count /), & |
---|
| 1325 | count = (/ out_noa /) ) |
---|
[3159] | 1326 | CALL netcdf_handle_error( 'mas_data_output_agents', 3 ) |
---|
| 1327 | |
---|
[3235] | 1328 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(3), agents%y, & |
---|
[3159] | 1329 | start = (/ 1, agt_time_count /), & |
---|
| 1330 | count = (/ out_noa /) ) |
---|
| 1331 | CALL netcdf_handle_error( 'mas_data_output_agents', 4 ) |
---|
| 1332 | |
---|
[3235] | 1333 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(4), agents%windspeed, & |
---|
[3159] | 1334 | start = (/ 1, agt_time_count /), & |
---|
| 1335 | count = (/ out_noa /) ) |
---|
| 1336 | CALL netcdf_handle_error( 'mas_data_output_agents', 5 ) |
---|
| 1337 | |
---|
[3235] | 1338 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(5), agents%t, & |
---|
[3159] | 1339 | start = (/ 1, agt_time_count /), & |
---|
| 1340 | count = (/ out_noa /) ) |
---|
| 1341 | CALL netcdf_handle_error( 'mas_data_output_agents', 6 ) |
---|
| 1342 | |
---|
[3235] | 1343 | nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(6), agents%group, & |
---|
[3159] | 1344 | start = (/ 1, agt_time_count /), & |
---|
| 1345 | count = (/ out_noa /) ) |
---|
| 1346 | CALL netcdf_handle_error( 'mas_data_output_agents', 7 ) |
---|
[3235] | 1347 | CALL cpu_log( log_point_s(17), 'mas_data_output', 'stop' ) |
---|
[3159] | 1348 | |
---|
[3165] | 1349 | |
---|
[3159] | 1350 | ! nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(6), agents%pm10, & |
---|
| 1351 | ! start = (/ 1, agt_time_count /), & |
---|
| 1352 | ! count = (/ out_noa /) ) |
---|
| 1353 | ! CALL netcdf_handle_error( 'mas_data_output_agents', 8 ) |
---|
| 1354 | ! |
---|
| 1355 | ! nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(7), agents%pm25, & |
---|
| 1356 | ! start = (/ 1, agt_time_count /), & |
---|
| 1357 | ! count = (/ out_noa /) ) |
---|
| 1358 | ! CALL netcdf_handle_error( 'mas_data_output_agents', 9 ) |
---|
| 1359 | ! |
---|
| 1360 | ! nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(8), agents%therm_comf,& |
---|
| 1361 | ! start = (/ 1, agt_time_count /), & |
---|
| 1362 | ! count = (/ out_noa /) ) |
---|
| 1363 | ! |
---|
| 1364 | ! nc_stat = NF90_PUT_VAR( id_set_agt, id_var_agt(9), agents%uv, & |
---|
| 1365 | ! start = (/ 1, agt_time_count /), & |
---|
| 1366 | ! count = (/ out_noa /) ) |
---|
[3235] | 1367 | ! CALL netcdf_handle_error( 'mas_data_output_agents', 10 ) |
---|
[3159] | 1368 | |
---|
| 1369 | #endif |
---|
[3268] | 1370 | agents => NULL() |
---|
| 1371 | #if defined( __parallel ) |
---|
| 1372 | IF ( ALLOCATED( out_agents ) ) DEALLOCATE( out_agents ) |
---|
| 1373 | #endif |
---|
[3235] | 1374 | ELSE |
---|
| 1375 | CALL cpu_log( log_point_s(17), 'mas_data_output', 'stop' ) |
---|
[3159] | 1376 | ENDIF |
---|
| 1377 | |
---|
[3268] | 1378 | IF ( ALLOCATED( trf_agents ) ) DEALLOCATE( trf_agents ) |
---|
[3159] | 1379 | |
---|
| 1380 | END SUBROUTINE mas_data_output_agents |
---|
| 1381 | |
---|
| 1382 | !------------------------------------------------------------------------------! |
---|
| 1383 | ! Description: |
---|
| 1384 | ! ------------ |
---|
| 1385 | !> If an agent moves from one processor to another, this subroutine moves |
---|
| 1386 | !> the corresponding elements from the agent arrays of the old grid cells |
---|
| 1387 | !> to the agent arrays of the new grid cells. |
---|
| 1388 | !------------------------------------------------------------------------------! |
---|
| 1389 | SUBROUTINE mas_eh_add_agents_to_gridcell (agent_array) |
---|
| 1390 | |
---|
| 1391 | IMPLICIT NONE |
---|
| 1392 | |
---|
| 1393 | INTEGER(iwp) :: aindex !< dummy argument for new number of agents per grid box |
---|
| 1394 | INTEGER(iwp) :: ip !< grid index (x) of agent |
---|
| 1395 | INTEGER(iwp) :: jp !< grid index (x) of agent |
---|
| 1396 | INTEGER(iwp) :: n !< index variable of agent |
---|
| 1397 | |
---|
| 1398 | LOGICAL :: pack_done !< flag to indicate that packing is done |
---|
| 1399 | |
---|
| 1400 | TYPE(agent_type), DIMENSION(:), INTENT(IN) :: agent_array !< new agents in a grid box |
---|
| 1401 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: temp_ns !< temporary agent array for reallocation |
---|
| 1402 | |
---|
| 1403 | pack_done = .FALSE. |
---|
| 1404 | |
---|
| 1405 | DO n = 1, SIZE(agent_array) |
---|
| 1406 | |
---|
| 1407 | IF ( .NOT. agent_array(n)%agent_mask ) CYCLE |
---|
| 1408 | |
---|
| 1409 | ip = agent_array(n)%x * ddx |
---|
| 1410 | jp = agent_array(n)%y * ddy |
---|
| 1411 | |
---|
| 1412 | IF ( ip >= nxl .AND. ip <= nxr .AND. & |
---|
| 1413 | jp >= nys .AND. jp <= nyn ) & |
---|
| 1414 | THEN ! agent stays on processor |
---|
| 1415 | number_of_agents = agt_count(jp,ip) |
---|
| 1416 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1417 | |
---|
| 1418 | aindex = agt_count(jp,ip)+1 |
---|
| 1419 | IF( aindex > SIZE(grid_agents(jp,ip)%agents) ) THEN |
---|
| 1420 | IF ( pack_done ) THEN |
---|
| 1421 | CALL mas_eh_realloc_agents_array (ip,jp) |
---|
| 1422 | ELSE |
---|
| 1423 | CALL mas_ps_pack |
---|
| 1424 | agt_count(jp,ip) = number_of_agents |
---|
| 1425 | aindex = agt_count(jp,ip)+1 |
---|
| 1426 | IF ( aindex > SIZE(grid_agents(jp,ip)%agents) ) THEN |
---|
| 1427 | CALL mas_eh_realloc_agents_array (ip,jp) |
---|
| 1428 | ENDIF |
---|
| 1429 | pack_done = .TRUE. |
---|
| 1430 | ENDIF |
---|
| 1431 | ENDIF |
---|
| 1432 | grid_agents(jp,ip)%agents(aindex) = agent_array(n) |
---|
| 1433 | agt_count(jp,ip) = aindex |
---|
| 1434 | ELSE |
---|
| 1435 | IF ( jp <= nys - 1 ) THEN |
---|
| 1436 | nr_move_south = nr_move_south+1 |
---|
| 1437 | ! |
---|
| 1438 | !-- Before agent information is swapped to exchange-array, check |
---|
| 1439 | !-- if enough memory is allocated. If required, reallocate exchange |
---|
| 1440 | !-- array. |
---|
| 1441 | IF ( nr_move_south > SIZE(move_also_south) ) THEN |
---|
| 1442 | ! |
---|
| 1443 | !-- At first, allocate further temporary array to swap agent |
---|
| 1444 | !-- information. |
---|
| 1445 | ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move)) |
---|
| 1446 | temp_ns(1:nr_move_south-1) = move_also_south & |
---|
| 1447 | (1:nr_move_south-1) |
---|
| 1448 | DEALLOCATE( move_also_south ) |
---|
| 1449 | ALLOCATE( move_also_south(SIZE(temp_ns)) ) |
---|
| 1450 | move_also_south(1:nr_move_south-1) = temp_ns & |
---|
| 1451 | (1:nr_move_south-1) |
---|
| 1452 | DEALLOCATE( temp_ns ) |
---|
| 1453 | |
---|
| 1454 | ENDIF |
---|
| 1455 | |
---|
| 1456 | move_also_south(nr_move_south) = agent_array(n) |
---|
| 1457 | |
---|
| 1458 | IF ( jp == -1 ) THEN |
---|
| 1459 | ! |
---|
| 1460 | !-- Apply boundary condition along y |
---|
| 1461 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 1462 | move_also_south(nr_move_south)%y = & |
---|
| 1463 | move_also_south(nr_move_south)%y & |
---|
| 1464 | + ( ny + 1 ) * dy |
---|
| 1465 | move_also_south(nr_move_south)%origin_y = & |
---|
| 1466 | move_also_south(nr_move_south)%origin_y & |
---|
| 1467 | + ( ny + 1 ) * dy |
---|
| 1468 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 1469 | ! |
---|
| 1470 | !-- Agent absorption |
---|
| 1471 | move_also_south(nr_move_south)%agent_mask = .FALSE. |
---|
| 1472 | deleted_agents = deleted_agents + 1 |
---|
| 1473 | |
---|
| 1474 | ENDIF |
---|
| 1475 | ENDIF |
---|
| 1476 | ELSEIF ( jp >= nyn+1 ) THEN |
---|
| 1477 | nr_move_north = nr_move_north+1 |
---|
| 1478 | ! |
---|
| 1479 | !-- Before agent information is swapped to exchange-array, check |
---|
| 1480 | !-- if enough memory is allocated. If required, reallocate exchange |
---|
| 1481 | !-- array. |
---|
| 1482 | IF ( nr_move_north > SIZE(move_also_north) ) THEN |
---|
| 1483 | ! |
---|
| 1484 | !-- At first, allocate further temporary array to swap agent |
---|
| 1485 | !-- information. |
---|
| 1486 | ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move)) |
---|
| 1487 | temp_ns(1:nr_move_north-1) = & |
---|
| 1488 | move_also_south(1:nr_move_north-1) |
---|
| 1489 | DEALLOCATE( move_also_north ) |
---|
| 1490 | ALLOCATE( move_also_north(SIZE(temp_ns)) ) |
---|
| 1491 | move_also_north(1:nr_move_north-1) = & |
---|
| 1492 | temp_ns(1:nr_move_north-1) |
---|
| 1493 | DEALLOCATE( temp_ns ) |
---|
| 1494 | |
---|
| 1495 | ENDIF |
---|
| 1496 | |
---|
| 1497 | move_also_north(nr_move_north) = agent_array(n) |
---|
| 1498 | IF ( jp == ny+1 ) THEN |
---|
| 1499 | ! |
---|
| 1500 | !-- Apply boundary condition along y |
---|
| 1501 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 1502 | |
---|
| 1503 | move_also_north(nr_move_north)%y = & |
---|
| 1504 | move_also_north(nr_move_north)%y & |
---|
| 1505 | - ( ny + 1 ) * dy |
---|
| 1506 | move_also_north(nr_move_north)%origin_y = & |
---|
| 1507 | move_also_north(nr_move_north)%origin_y & |
---|
| 1508 | - ( ny + 1 ) * dy |
---|
| 1509 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 1510 | ! |
---|
| 1511 | !-- Agent absorption |
---|
| 1512 | move_also_north(nr_move_north)%agent_mask = .FALSE. |
---|
| 1513 | deleted_agents = deleted_agents + 1 |
---|
| 1514 | |
---|
| 1515 | ENDIF |
---|
| 1516 | ENDIF |
---|
| 1517 | ENDIF |
---|
| 1518 | ENDIF |
---|
| 1519 | ENDDO |
---|
| 1520 | |
---|
| 1521 | RETURN |
---|
| 1522 | |
---|
| 1523 | END SUBROUTINE mas_eh_add_agents_to_gridcell |
---|
| 1524 | |
---|
| 1525 | !------------------------------------------------------------------------------! |
---|
| 1526 | ! Description: |
---|
| 1527 | ! ------------ |
---|
| 1528 | !> After ghost layer agents have been received from neighboring PEs, this |
---|
| 1529 | !> subroutine sorts them into the corresponding grid cells |
---|
| 1530 | !------------------------------------------------------------------------------! |
---|
| 1531 | SUBROUTINE mas_eh_add_ghost_agents_to_gridcell (agent_array) |
---|
| 1532 | |
---|
| 1533 | IMPLICIT NONE |
---|
| 1534 | |
---|
| 1535 | INTEGER(iwp) :: ip !< grid index (x) of agent |
---|
| 1536 | INTEGER(iwp) :: jp !< grid index (x) of agent |
---|
| 1537 | INTEGER(iwp) :: n !< index variable of agent |
---|
| 1538 | INTEGER(iwp) :: aindex !< dummy argument for new number of agents per grid box |
---|
| 1539 | |
---|
| 1540 | LOGICAL :: pack_done !< flag to indicate that packing is done |
---|
| 1541 | |
---|
| 1542 | TYPE(agent_type), DIMENSION(:), INTENT(IN) :: agent_array !< new agents in a grid box |
---|
| 1543 | |
---|
| 1544 | pack_done = .FALSE. |
---|
| 1545 | |
---|
| 1546 | DO n = 1, SIZE(agent_array) |
---|
| 1547 | |
---|
| 1548 | IF ( .NOT. agent_array(n)%agent_mask ) CYCLE |
---|
| 1549 | |
---|
| 1550 | ip = agent_array(n)%x * ddx |
---|
| 1551 | jp = agent_array(n)%y * ddy |
---|
| 1552 | |
---|
| 1553 | IF ( ip < nxl .OR. ip > nxr .OR. jp < nys .OR. jp > nyn ) THEN |
---|
| 1554 | number_of_agents = agt_count(jp,ip) |
---|
| 1555 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1556 | |
---|
| 1557 | aindex = agt_count(jp,ip)+1 |
---|
| 1558 | IF( aindex > SIZE(grid_agents(jp,ip)%agents) ) THEN |
---|
| 1559 | IF ( pack_done ) THEN |
---|
| 1560 | CALL mas_eh_realloc_agents_array (ip,jp) |
---|
| 1561 | ELSE |
---|
| 1562 | CALL mas_ps_pack |
---|
| 1563 | agt_count(jp,ip) = number_of_agents |
---|
| 1564 | aindex = agt_count(jp,ip)+1 |
---|
| 1565 | IF ( aindex > SIZE(grid_agents(jp,ip)%agents) ) THEN |
---|
| 1566 | CALL mas_eh_realloc_agents_array (ip,jp) |
---|
| 1567 | ENDIF |
---|
| 1568 | pack_done = .TRUE. |
---|
| 1569 | ENDIF |
---|
| 1570 | ENDIF |
---|
| 1571 | grid_agents(jp,ip)%agents(aindex) = agent_array(n) |
---|
| 1572 | agt_count(jp,ip) = aindex |
---|
| 1573 | ENDIF |
---|
| 1574 | ENDDO |
---|
| 1575 | END SUBROUTINE mas_eh_add_ghost_agents_to_gridcell |
---|
| 1576 | |
---|
| 1577 | !------------------------------------------------------------------------------! |
---|
| 1578 | ! Description: |
---|
| 1579 | ! ------------ |
---|
| 1580 | !> Resizing of agent arrays |
---|
| 1581 | !------------------------------------------------------------------------------! |
---|
| 1582 | SUBROUTINE mas_eh_dealloc_agents_array |
---|
| 1583 | |
---|
| 1584 | IMPLICIT NONE |
---|
| 1585 | |
---|
| 1586 | INTEGER(iwp) :: i !< grid index (x) of agent |
---|
| 1587 | INTEGER(iwp) :: j !< grid index (y) of agent |
---|
| 1588 | INTEGER(iwp) :: old_size !< old array size |
---|
| 1589 | INTEGER(iwp) :: new_size !< new array size |
---|
| 1590 | INTEGER(iwp) :: noa !< number of agents |
---|
| 1591 | |
---|
| 1592 | LOGICAL :: dealloc !< flag that indicates if reallocation is necessary |
---|
| 1593 | |
---|
| 1594 | TYPE(agent_type), DIMENSION(10) :: tmp_agents_s !< temporary static agent array |
---|
| 1595 | |
---|
| 1596 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: tmp_agents_d !< temporary dynamic agent array |
---|
| 1597 | |
---|
| 1598 | DO i = nxlg, nxrg |
---|
| 1599 | DO j = nysg, nyng |
---|
| 1600 | ! |
---|
| 1601 | !-- Determine number of active agents |
---|
| 1602 | noa = agt_count(j,i) |
---|
| 1603 | ! |
---|
| 1604 | !-- Determine allocated memory size |
---|
| 1605 | old_size = SIZE( grid_agents(j,i)%agents ) |
---|
| 1606 | ! |
---|
| 1607 | !-- Check for large unused memory |
---|
| 1608 | dealloc = ( ( noa < min_nr_agent .AND. old_size > min_nr_agent ) & |
---|
| 1609 | .OR. ( noa > min_nr_agent .AND. old_size - noa * & |
---|
| 1610 | ( 1.0_wp + 0.01_wp * alloc_factor_mas ) > 0.0_wp ) ) |
---|
| 1611 | ! |
---|
| 1612 | !-- If large unused memory was found, resize the corresponding array |
---|
| 1613 | IF ( dealloc ) THEN |
---|
| 1614 | IF ( noa < min_nr_agent ) THEN |
---|
| 1615 | new_size = min_nr_agent |
---|
| 1616 | ELSE |
---|
| 1617 | new_size = INT( noa * ( 1.0_wp + & |
---|
| 1618 | 0.01_wp * alloc_factor_mas ) ) |
---|
| 1619 | ENDIF |
---|
| 1620 | |
---|
| 1621 | IF ( noa <= 10 ) THEN |
---|
| 1622 | |
---|
| 1623 | tmp_agents_s(1:noa) = grid_agents(j,i)%agents(1:noa) |
---|
| 1624 | |
---|
| 1625 | DEALLOCATE(grid_agents(j,i)%agents) |
---|
| 1626 | ALLOCATE(grid_agents(j,i)%agents(1:new_size)) |
---|
| 1627 | |
---|
| 1628 | grid_agents(j,i)%agents(1:noa) = tmp_agents_s(1:noa) |
---|
| 1629 | grid_agents(j,i)%agents(noa+1:new_size) = zero_agent |
---|
| 1630 | |
---|
| 1631 | ELSE |
---|
| 1632 | |
---|
| 1633 | ALLOCATE(tmp_agents_d(noa)) |
---|
| 1634 | tmp_agents_d(1:noa) = grid_agents(j,i)%agents(1:noa) |
---|
| 1635 | |
---|
| 1636 | DEALLOCATE(grid_agents(j,i)%agents) |
---|
| 1637 | ALLOCATE(grid_agents(j,i)%agents(new_size)) |
---|
| 1638 | |
---|
| 1639 | grid_agents(j,i)%agents(1:noa) = tmp_agents_d(1:noa) |
---|
| 1640 | grid_agents(j,i)%agents(noa+1:new_size) = zero_agent |
---|
| 1641 | |
---|
| 1642 | DEALLOCATE(tmp_agents_d) |
---|
| 1643 | |
---|
| 1644 | ENDIF |
---|
| 1645 | |
---|
| 1646 | ENDIF |
---|
| 1647 | ENDDO |
---|
| 1648 | ENDDO |
---|
| 1649 | |
---|
| 1650 | END SUBROUTINE mas_eh_dealloc_agents_array |
---|
| 1651 | |
---|
| 1652 | !------------------------------------------------------------------------------! |
---|
| 1653 | ! Description: |
---|
| 1654 | ! ------------ |
---|
| 1655 | !> Exchange between subdomains. |
---|
| 1656 | !> As soon as one agent has moved beyond the boundary of the domain, it |
---|
| 1657 | !> is included in the relevant transfer arrays and marked for subsequent |
---|
| 1658 | !> deletion on this PE. |
---|
| 1659 | !> First sweep for crossings in x direction. Find out first the number of |
---|
| 1660 | !> agents to be transferred and allocate temporary arrays needed to store |
---|
| 1661 | !> them. |
---|
| 1662 | !> For a one-dimensional decomposition along y, no transfer is necessary, |
---|
| 1663 | !> because the agent remains on the PE, but the agent coordinate has to |
---|
| 1664 | !> be adjusted. |
---|
| 1665 | !------------------------------------------------------------------------------! |
---|
| 1666 | SUBROUTINE mas_eh_exchange_horiz |
---|
| 1667 | |
---|
| 1668 | IMPLICIT NONE |
---|
| 1669 | |
---|
| 1670 | INTEGER(iwp) :: i !< grid index (x) of agent positition |
---|
| 1671 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 1672 | INTEGER(iwp) :: j !< grid index (y) of agent positition |
---|
| 1673 | INTEGER(iwp) :: jp !< index variable along y |
---|
| 1674 | INTEGER(iwp) :: n !< agent index variable |
---|
| 1675 | INTEGER(iwp) :: par_size !< Agent size in bytes |
---|
| 1676 | INTEGER(iwp) :: trla_count !< number of agents send to left PE |
---|
| 1677 | INTEGER(iwp) :: trla_count_recv !< number of agents receive from right PE |
---|
| 1678 | INTEGER(iwp) :: trna_count !< number of agents send to north PE |
---|
| 1679 | INTEGER(iwp) :: trna_count_recv !< number of agents receive from south PE |
---|
| 1680 | INTEGER(iwp) :: trra_count !< number of agents send to right PE |
---|
| 1681 | INTEGER(iwp) :: trra_count_recv !< number of agents receive from left PE |
---|
| 1682 | INTEGER(iwp) :: trsa_count !< number of agents send to south PE |
---|
| 1683 | INTEGER(iwp) :: trsa_count_recv !< number of agents receive from north PE |
---|
| 1684 | |
---|
| 1685 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: rvla !< agents received from right PE |
---|
| 1686 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: rvna !< agents received from south PE |
---|
| 1687 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: rvra !< agents received from left PE |
---|
| 1688 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: rvsa !< agents received from north PE |
---|
| 1689 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: trla !< agents send to left PE |
---|
| 1690 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: trna !< agents send to north PE |
---|
| 1691 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: trra !< agents send to right PE |
---|
| 1692 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: trsa !< agents send to south PE |
---|
| 1693 | |
---|
| 1694 | #if defined( __parallel ) |
---|
| 1695 | |
---|
| 1696 | ! |
---|
| 1697 | !-- Exchange between subdomains. |
---|
| 1698 | !-- As soon as one agent has moved beyond the boundary of the domain, it |
---|
| 1699 | !-- is included in the relevant transfer arrays and marked for subsequent |
---|
| 1700 | !-- deletion on this PE. |
---|
| 1701 | !-- First sweep for crossings in x direction. Find out first the number of |
---|
| 1702 | !-- agents to be transferred and allocate temporary arrays needed to store |
---|
| 1703 | !-- them. |
---|
| 1704 | !-- For a one-dimensional decomposition along y, no transfer is necessary, |
---|
| 1705 | !-- because the agent remains on the PE, but the agent coordinate has to |
---|
| 1706 | !-- be adjusted. |
---|
| 1707 | trla_count = 0 |
---|
| 1708 | trra_count = 0 |
---|
| 1709 | |
---|
| 1710 | trla_count_recv = 0 |
---|
| 1711 | trra_count_recv = 0 |
---|
| 1712 | |
---|
| 1713 | IF ( pdims(1) /= 1 ) THEN |
---|
| 1714 | ! |
---|
| 1715 | !-- First calculate the storage necessary for sending and receiving the data. |
---|
| 1716 | !-- Compute only first (nxl) and last (nxr) loop iterration. |
---|
| 1717 | DO ip = nxl, nxr, nxr - nxl |
---|
| 1718 | DO jp = nys, nyn |
---|
| 1719 | |
---|
| 1720 | number_of_agents = agt_count(jp,ip) |
---|
| 1721 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 1722 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1723 | DO n = 1, number_of_agents |
---|
| 1724 | IF ( agents(n)%agent_mask ) THEN |
---|
| 1725 | i = agents(n)%x * ddx |
---|
| 1726 | ! |
---|
| 1727 | !-- Above calculation does not work for indices less than zero |
---|
| 1728 | IF ( agents(n)%x < 0.0_wp ) i = -1 |
---|
| 1729 | |
---|
| 1730 | IF ( i < nxl ) THEN |
---|
| 1731 | trla_count = trla_count + 1 |
---|
| 1732 | ELSEIF ( i > nxr ) THEN |
---|
| 1733 | trra_count = trra_count + 1 |
---|
| 1734 | ENDIF |
---|
| 1735 | ENDIF |
---|
| 1736 | ENDDO |
---|
| 1737 | |
---|
| 1738 | ENDDO |
---|
| 1739 | ENDDO |
---|
| 1740 | |
---|
| 1741 | IF ( trla_count == 0 ) trla_count = 1 |
---|
| 1742 | IF ( trra_count == 0 ) trra_count = 1 |
---|
| 1743 | |
---|
| 1744 | ALLOCATE( trla(trla_count), trra(trra_count) ) |
---|
| 1745 | |
---|
| 1746 | trla = zero_agent |
---|
| 1747 | trra = zero_agent |
---|
| 1748 | |
---|
| 1749 | trla_count = 0 |
---|
| 1750 | trra_count = 0 |
---|
| 1751 | |
---|
| 1752 | ENDIF |
---|
| 1753 | ! |
---|
| 1754 | !-- Compute only first (nxl) and last (nxr) loop iterration |
---|
| 1755 | DO ip = nxl, nxr, nxr-nxl |
---|
| 1756 | DO jp = nys, nyn |
---|
| 1757 | number_of_agents = agt_count(jp,ip) |
---|
| 1758 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 1759 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1760 | DO n = 1, number_of_agents |
---|
| 1761 | ! |
---|
| 1762 | !-- Only those agents that have not been marked as 'deleted' may |
---|
| 1763 | !-- be moved. |
---|
| 1764 | IF ( agents(n)%agent_mask ) THEN |
---|
| 1765 | |
---|
| 1766 | i = agents(n)%x * ddx |
---|
| 1767 | ! |
---|
| 1768 | !-- Above calculation does not work for indices less than zero |
---|
| 1769 | IF ( agents(n)%x < 0.0_wp ) i = -1 |
---|
| 1770 | |
---|
| 1771 | IF ( i < nxl ) THEN |
---|
| 1772 | IF ( i < 0 ) THEN |
---|
| 1773 | ! |
---|
| 1774 | !-- Apply boundary condition along x |
---|
| 1775 | IF ( ibc_mas_lr == 0 ) THEN |
---|
| 1776 | ! |
---|
| 1777 | !-- Cyclic condition |
---|
| 1778 | IF ( pdims(1) == 1 ) THEN |
---|
| 1779 | agents(n)%x = ( nx + 1 ) * dx + & |
---|
| 1780 | agents(n)%x |
---|
| 1781 | agents(n)%origin_x = ( nx + 1 ) * dx + & |
---|
| 1782 | agents(n)%origin_x |
---|
| 1783 | ELSE |
---|
| 1784 | trla_count = trla_count + 1 |
---|
| 1785 | trla(trla_count) = agents(n) |
---|
| 1786 | trla(trla_count)%x = ( nx + 1 ) * dx + & |
---|
| 1787 | trla(trla_count)%x |
---|
| 1788 | trla(trla_count)%origin_x = & |
---|
| 1789 | trla(trla_count)%origin_x + & |
---|
| 1790 | ( nx + 1 ) * dx |
---|
| 1791 | agents(n)%agent_mask = .FALSE. |
---|
| 1792 | deleted_agents = deleted_agents + 1 |
---|
| 1793 | |
---|
| 1794 | IF ( trla(trla_count)%x >= & |
---|
| 1795 | (nx + 1)* dx - 1.0E-12_wp ) & |
---|
| 1796 | THEN |
---|
| 1797 | trla(trla_count)%x = trla(trla_count)%x - & |
---|
| 1798 | 1.0E-10_wp |
---|
| 1799 | trla(trla_count)%origin_x = & |
---|
| 1800 | trla(trla_count)%origin_x - 1 |
---|
| 1801 | ENDIF |
---|
| 1802 | |
---|
| 1803 | ENDIF |
---|
| 1804 | |
---|
| 1805 | ELSEIF ( ibc_mas_lr == 1 ) THEN |
---|
| 1806 | ! |
---|
| 1807 | !-- Agent absorption |
---|
| 1808 | agents(n)%agent_mask = .FALSE. |
---|
| 1809 | deleted_agents = deleted_agents + 1 |
---|
| 1810 | |
---|
| 1811 | ENDIF |
---|
| 1812 | ELSE |
---|
| 1813 | ! |
---|
| 1814 | !-- Store agent data in the transfer array, which will be |
---|
| 1815 | !-- send to the neighbouring PE |
---|
| 1816 | trla_count = trla_count + 1 |
---|
| 1817 | trla(trla_count) = agents(n) |
---|
| 1818 | agents(n)%agent_mask = .FALSE. |
---|
| 1819 | deleted_agents = deleted_agents + 1 |
---|
| 1820 | |
---|
| 1821 | ENDIF |
---|
| 1822 | |
---|
| 1823 | ELSEIF ( i > nxr ) THEN |
---|
| 1824 | IF ( i > nx ) THEN |
---|
| 1825 | ! |
---|
| 1826 | !-- Apply boundary condition along x |
---|
| 1827 | IF ( ibc_mas_lr == 0 ) THEN |
---|
| 1828 | ! |
---|
| 1829 | !-- Cyclic condition |
---|
| 1830 | IF ( pdims(1) == 1 ) THEN |
---|
| 1831 | agents(n)%x = agents(n)%x - ( nx + 1 ) * dx |
---|
| 1832 | agents(n)%origin_x = agents(n)%origin_x - & |
---|
| 1833 | ( nx + 1 ) * dx |
---|
| 1834 | ELSE |
---|
| 1835 | trra_count = trra_count + 1 |
---|
| 1836 | trra(trra_count) = agents(n) |
---|
| 1837 | trra(trra_count)%x = trra(trra_count)%x - & |
---|
| 1838 | ( nx + 1 ) * dx |
---|
| 1839 | trra(trra_count)%origin_x = & |
---|
| 1840 | trra(trra_count)%origin_x - & |
---|
| 1841 | ( nx + 1 ) * dx |
---|
| 1842 | agents(n)%agent_mask = .FALSE. |
---|
| 1843 | deleted_agents = deleted_agents + 1 |
---|
| 1844 | |
---|
| 1845 | ENDIF |
---|
| 1846 | |
---|
| 1847 | ELSEIF ( ibc_mas_lr == 1 ) THEN |
---|
| 1848 | ! |
---|
| 1849 | !-- Agent absorption |
---|
| 1850 | agents(n)%agent_mask = .FALSE. |
---|
| 1851 | deleted_agents = deleted_agents + 1 |
---|
| 1852 | |
---|
| 1853 | ENDIF |
---|
| 1854 | ELSE |
---|
| 1855 | ! |
---|
| 1856 | !-- Store agent data in the transfer array, which will be send |
---|
| 1857 | !-- to the neighbouring PE |
---|
| 1858 | trra_count = trra_count + 1 |
---|
| 1859 | trra(trra_count) = agents(n) |
---|
| 1860 | agents(n)%agent_mask = .FALSE. |
---|
| 1861 | deleted_agents = deleted_agents + 1 |
---|
| 1862 | |
---|
| 1863 | ENDIF |
---|
| 1864 | |
---|
| 1865 | ENDIF |
---|
| 1866 | ENDIF |
---|
| 1867 | |
---|
| 1868 | ENDDO |
---|
| 1869 | ENDDO |
---|
| 1870 | ENDDO |
---|
| 1871 | |
---|
| 1872 | ! |
---|
| 1873 | !-- Allocate arrays required for north-south exchange, as these |
---|
| 1874 | !-- are used directly after agents are exchange along x-direction. |
---|
| 1875 | ALLOCATE( move_also_north(1:NR_2_direction_move) ) |
---|
| 1876 | ALLOCATE( move_also_south(1:NR_2_direction_move) ) |
---|
| 1877 | |
---|
| 1878 | nr_move_north = 0 |
---|
| 1879 | nr_move_south = 0 |
---|
| 1880 | ! |
---|
| 1881 | !-- Send left boundary, receive right boundary (but first exchange how many |
---|
| 1882 | !-- and chec if agent storage must be extended) |
---|
| 1883 | IF ( pdims(1) /= 1 ) THEN |
---|
| 1884 | |
---|
| 1885 | CALL MPI_SENDRECV( trla_count, 1, MPI_INTEGER, pleft, 0, & |
---|
| 1886 | trra_count_recv, 1, MPI_INTEGER, pright, 0, & |
---|
| 1887 | comm2d, status, ierr ) |
---|
| 1888 | |
---|
| 1889 | ALLOCATE(rvra(MAX(1,trra_count_recv))) |
---|
| 1890 | ! |
---|
| 1891 | !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit |
---|
| 1892 | !-- variables in structure agent_type (due to the calculation of par_size) |
---|
[3160] | 1893 | par_size = STORAGE_SIZE(trla(1))/8 |
---|
[3159] | 1894 | CALL MPI_SENDRECV( trla, max(1,trla_count)*par_size, MPI_BYTE, pleft,& |
---|
| 1895 | 1, rvra, max(1,trra_count_recv)*par_size, MPI_BYTE, pright,& |
---|
| 1896 | 1, comm2d, status, ierr ) |
---|
| 1897 | |
---|
| 1898 | IF ( trra_count_recv > 0 ) THEN |
---|
| 1899 | CALL mas_eh_add_agents_to_gridcell(rvra(1:trra_count_recv)) |
---|
| 1900 | ENDIF |
---|
| 1901 | |
---|
| 1902 | DEALLOCATE(rvra) |
---|
| 1903 | |
---|
| 1904 | ! |
---|
| 1905 | !-- Send right boundary, receive left boundary |
---|
| 1906 | CALL MPI_SENDRECV( trra_count, 1, MPI_INTEGER, pright, 0, & |
---|
| 1907 | trla_count_recv, 1, MPI_INTEGER, pleft, 0, & |
---|
| 1908 | comm2d, status, ierr ) |
---|
| 1909 | |
---|
| 1910 | ALLOCATE(rvla(MAX(1,trla_count_recv))) |
---|
| 1911 | ! |
---|
| 1912 | !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit |
---|
| 1913 | !-- variables in structure agent_type (due to the calculation of par_size) |
---|
[3160] | 1914 | par_size = STORAGE_SIZE(trra(1))/8 |
---|
[3159] | 1915 | CALL MPI_SENDRECV( trra, max(1,trra_count)*par_size, MPI_BYTE, & |
---|
| 1916 | pright, 1, rvla, & |
---|
| 1917 | max(1,trla_count_recv)*par_size, MPI_BYTE, & |
---|
| 1918 | pleft, 1, comm2d, status, ierr ) |
---|
| 1919 | |
---|
| 1920 | IF ( trla_count_recv > 0 ) THEN |
---|
| 1921 | CALL mas_eh_add_agents_to_gridcell(rvla(1:trla_count_recv)) |
---|
| 1922 | ENDIF |
---|
| 1923 | |
---|
| 1924 | DEALLOCATE( rvla ) |
---|
| 1925 | DEALLOCATE( trla, trra ) |
---|
| 1926 | |
---|
| 1927 | ENDIF |
---|
| 1928 | |
---|
| 1929 | ! |
---|
| 1930 | !-- Check whether agents have crossed the boundaries in y direction. Note |
---|
| 1931 | !-- that this case can also apply to agents that have just been received |
---|
| 1932 | !-- from the adjacent right or left PE. |
---|
| 1933 | !-- Find out first the number of agents to be transferred and allocate |
---|
| 1934 | !-- temporary arrays needed to store them. |
---|
| 1935 | !-- For a one-dimensional decomposition along y, no transfer is necessary, |
---|
| 1936 | !-- because the agent remains on the PE. |
---|
| 1937 | trsa_count = nr_move_south |
---|
| 1938 | trna_count = nr_move_north |
---|
| 1939 | |
---|
| 1940 | trsa_count_recv = 0 |
---|
| 1941 | trna_count_recv = 0 |
---|
| 1942 | |
---|
| 1943 | IF ( pdims(2) /= 1 ) THEN |
---|
| 1944 | ! |
---|
| 1945 | !-- First calculate the storage necessary for sending and receiving the |
---|
| 1946 | !-- data |
---|
| 1947 | DO ip = nxl, nxr |
---|
| 1948 | DO jp = nys, nyn, nyn-nys !compute only first (nys) and last (nyn) loop iterration |
---|
| 1949 | number_of_agents = agt_count(jp,ip) |
---|
| 1950 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 1951 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1952 | DO n = 1, number_of_agents |
---|
| 1953 | IF ( agents(n)%agent_mask ) THEN |
---|
| 1954 | j = agents(n)%y * ddy |
---|
| 1955 | ! |
---|
| 1956 | !-- Above calculation does not work for indices less than zero |
---|
| 1957 | IF ( agents(n)%y < 0.0_wp ) j = -1 |
---|
| 1958 | |
---|
| 1959 | IF ( j < nys ) THEN |
---|
| 1960 | trsa_count = trsa_count + 1 |
---|
| 1961 | ELSEIF ( j > nyn ) THEN |
---|
| 1962 | trna_count = trna_count + 1 |
---|
| 1963 | ENDIF |
---|
| 1964 | ENDIF |
---|
| 1965 | ENDDO |
---|
| 1966 | ENDDO |
---|
| 1967 | ENDDO |
---|
| 1968 | |
---|
| 1969 | IF ( trsa_count == 0 ) trsa_count = 1 |
---|
| 1970 | IF ( trna_count == 0 ) trna_count = 1 |
---|
| 1971 | |
---|
| 1972 | ALLOCATE( trsa(trsa_count), trna(trna_count) ) |
---|
| 1973 | |
---|
| 1974 | trsa = zero_agent |
---|
| 1975 | trna = zero_agent |
---|
| 1976 | |
---|
| 1977 | trsa_count = nr_move_south |
---|
| 1978 | trna_count = nr_move_north |
---|
| 1979 | |
---|
| 1980 | trsa(1:nr_move_south) = move_also_south(1:nr_move_south) |
---|
| 1981 | trna(1:nr_move_north) = move_also_north(1:nr_move_north) |
---|
| 1982 | |
---|
| 1983 | ENDIF |
---|
| 1984 | |
---|
| 1985 | DO ip = nxl, nxr |
---|
| 1986 | DO jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration |
---|
| 1987 | number_of_agents = agt_count(jp,ip) |
---|
| 1988 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 1989 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 1990 | DO n = 1, number_of_agents |
---|
| 1991 | ! |
---|
| 1992 | !-- Only those agents that have not been marked as 'deleted' may |
---|
| 1993 | !-- be moved. |
---|
| 1994 | IF ( agents(n)%agent_mask ) THEN |
---|
| 1995 | |
---|
| 1996 | j = agents(n)%y * ddy |
---|
| 1997 | ! |
---|
| 1998 | !-- Above calculation does not work for indices less than zero |
---|
| 1999 | IF ( agents(n)%y < 0.0_wp * dy ) j = -1 |
---|
| 2000 | |
---|
| 2001 | IF ( j < nys ) THEN |
---|
| 2002 | IF ( j < 0 ) THEN |
---|
| 2003 | ! |
---|
| 2004 | !-- Apply boundary condition along y |
---|
| 2005 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 2006 | ! |
---|
| 2007 | !-- Cyclic condition |
---|
| 2008 | IF ( pdims(2) == 1 ) THEN |
---|
| 2009 | agents(n)%y = ( ny + 1 ) * dy + agents(n)%y |
---|
| 2010 | agents(n)%origin_y = ( ny + 1 ) * dy + & |
---|
| 2011 | agents(n)%origin_y |
---|
| 2012 | ELSE |
---|
| 2013 | trsa_count = trsa_count + 1 |
---|
| 2014 | trsa(trsa_count) = agents(n) |
---|
| 2015 | trsa(trsa_count)%y = ( ny + 1 ) * dy + & |
---|
| 2016 | trsa(trsa_count)%y |
---|
| 2017 | trsa(trsa_count)%origin_y = & |
---|
| 2018 | trsa(trsa_count)%origin_y & |
---|
| 2019 | + ( ny + 1 ) * dy |
---|
| 2020 | agents(n)%agent_mask = .FALSE. |
---|
| 2021 | deleted_agents = deleted_agents + 1 |
---|
| 2022 | |
---|
| 2023 | IF ( trsa(trsa_count)%y >= & |
---|
| 2024 | (ny+1)* dy - 1.0E-12_wp ) & |
---|
| 2025 | THEN |
---|
| 2026 | trsa(trsa_count)%y = trsa(trsa_count)%y - & |
---|
| 2027 | 1.0E-10_wp |
---|
| 2028 | trsa(trsa_count)%origin_y = & |
---|
| 2029 | trsa(trsa_count)%origin_y - 1 |
---|
| 2030 | ENDIF |
---|
| 2031 | |
---|
| 2032 | ENDIF |
---|
| 2033 | |
---|
| 2034 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 2035 | ! |
---|
| 2036 | !-- Agent absorption |
---|
| 2037 | agents(n)%agent_mask = .FALSE. |
---|
| 2038 | deleted_agents = deleted_agents + 1 |
---|
| 2039 | |
---|
| 2040 | ENDIF |
---|
| 2041 | ELSE |
---|
| 2042 | ! |
---|
| 2043 | !-- Store agent data in the transfer array, which will |
---|
| 2044 | !-- be send to the neighbouring PE |
---|
| 2045 | trsa_count = trsa_count + 1 |
---|
| 2046 | trsa(trsa_count) = agents(n) |
---|
| 2047 | agents(n)%agent_mask = .FALSE. |
---|
| 2048 | deleted_agents = deleted_agents + 1 |
---|
| 2049 | |
---|
| 2050 | ENDIF |
---|
| 2051 | |
---|
| 2052 | ELSEIF ( j > nyn ) THEN |
---|
| 2053 | IF ( j > ny ) THEN |
---|
| 2054 | ! |
---|
| 2055 | !-- Apply boundary condition along y |
---|
| 2056 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 2057 | ! |
---|
| 2058 | !-- Cyclic condition |
---|
| 2059 | IF ( pdims(2) == 1 ) THEN |
---|
| 2060 | agents(n)%y = agents(n)%y - & |
---|
| 2061 | ( ny + 1 ) * dy |
---|
| 2062 | agents(n)%origin_y = agents(n)%origin_y - & |
---|
| 2063 | ( ny + 1 ) * dy |
---|
| 2064 | ELSE |
---|
| 2065 | trna_count = trna_count + 1 |
---|
| 2066 | trna(trna_count) = agents(n) |
---|
| 2067 | trna(trna_count)%y = & |
---|
| 2068 | trna(trna_count)%y - ( ny + 1 ) * dy |
---|
| 2069 | trna(trna_count)%origin_y = & |
---|
| 2070 | trna(trna_count)%origin_y - & |
---|
| 2071 | ( ny + 1 ) * dy |
---|
| 2072 | agents(n)%agent_mask = .FALSE. |
---|
| 2073 | deleted_agents = deleted_agents + 1 |
---|
| 2074 | ENDIF |
---|
| 2075 | |
---|
| 2076 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 2077 | ! |
---|
| 2078 | !-- Agent absorption |
---|
| 2079 | agents(n)%agent_mask = .FALSE. |
---|
| 2080 | deleted_agents = deleted_agents + 1 |
---|
| 2081 | |
---|
| 2082 | ENDIF |
---|
| 2083 | ELSE |
---|
| 2084 | ! |
---|
| 2085 | !-- Store agent data in the transfer array, which will |
---|
| 2086 | !-- be send to the neighbouring PE |
---|
| 2087 | trna_count = trna_count + 1 |
---|
| 2088 | trna(trna_count) = agents(n) |
---|
| 2089 | agents(n)%agent_mask = .FALSE. |
---|
| 2090 | deleted_agents = deleted_agents + 1 |
---|
| 2091 | |
---|
| 2092 | ENDIF |
---|
| 2093 | |
---|
| 2094 | ENDIF |
---|
| 2095 | ENDIF |
---|
| 2096 | ENDDO |
---|
| 2097 | ENDDO |
---|
| 2098 | ENDDO |
---|
| 2099 | |
---|
| 2100 | ! |
---|
| 2101 | !-- Send front boundary, receive back boundary (but first exchange how many |
---|
| 2102 | !-- and chec if agent storage must be extended) |
---|
| 2103 | IF ( pdims(2) /= 1 ) THEN |
---|
| 2104 | |
---|
| 2105 | CALL MPI_SENDRECV( trsa_count, 1, MPI_INTEGER, psouth, 0, & |
---|
| 2106 | trna_count_recv, 1, MPI_INTEGER, pnorth, 0, & |
---|
| 2107 | comm2d, status, ierr ) |
---|
| 2108 | |
---|
| 2109 | ALLOCATE(rvna(MAX(1,trna_count_recv))) |
---|
| 2110 | ! |
---|
| 2111 | !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit |
---|
| 2112 | !-- variables in structure agent_type (due to the calculation of par_size) |
---|
[3160] | 2113 | par_size = STORAGE_SIZE(trsa(1))/8 |
---|
[3159] | 2114 | CALL MPI_SENDRECV( trsa, trsa_count*par_size, MPI_BYTE, & |
---|
| 2115 | psouth, 1, rvna, & |
---|
| 2116 | trna_count_recv*par_size, MPI_BYTE, pnorth, 1, & |
---|
| 2117 | comm2d, status, ierr ) |
---|
| 2118 | |
---|
| 2119 | IF ( trna_count_recv > 0 ) THEN |
---|
| 2120 | CALL mas_eh_add_agents_to_gridcell(rvna(1:trna_count_recv)) |
---|
| 2121 | ENDIF |
---|
| 2122 | |
---|
| 2123 | DEALLOCATE(rvna) |
---|
| 2124 | |
---|
| 2125 | ! |
---|
| 2126 | !-- Send back boundary, receive front boundary |
---|
| 2127 | CALL MPI_SENDRECV( trna_count, 1, MPI_INTEGER, pnorth, 0, & |
---|
| 2128 | trsa_count_recv, 1, MPI_INTEGER, psouth, 0, & |
---|
| 2129 | comm2d, status, ierr ) |
---|
| 2130 | |
---|
| 2131 | ALLOCATE(rvsa(MAX(1,trsa_count_recv))) |
---|
| 2132 | ! |
---|
| 2133 | !-- This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit |
---|
| 2134 | !-- variables in structure agent_type (due to the calculation of par_size) |
---|
[3160] | 2135 | par_size = STORAGE_SIZE(trna(1))/8 |
---|
[3159] | 2136 | CALL MPI_SENDRECV( trna, trna_count*par_size, MPI_BYTE, & |
---|
| 2137 | pnorth, 1, rvsa, & |
---|
| 2138 | trsa_count_recv*par_size, MPI_BYTE, psouth, 1, & |
---|
| 2139 | comm2d, status, ierr ) |
---|
| 2140 | |
---|
| 2141 | IF ( trsa_count_recv > 0 ) THEN |
---|
| 2142 | CALL mas_eh_add_agents_to_gridcell(rvsa(1:trsa_count_recv)) |
---|
| 2143 | ENDIF |
---|
| 2144 | |
---|
| 2145 | DEALLOCATE(rvsa) |
---|
| 2146 | |
---|
| 2147 | number_of_agents = number_of_agents + trsa_count_recv |
---|
| 2148 | |
---|
| 2149 | DEALLOCATE( trsa, trna ) |
---|
| 2150 | |
---|
| 2151 | ENDIF |
---|
| 2152 | |
---|
| 2153 | DEALLOCATE( move_also_north ) |
---|
| 2154 | DEALLOCATE( move_also_south ) |
---|
| 2155 | |
---|
| 2156 | #else |
---|
| 2157 | |
---|
| 2158 | DO ip = nxl, nxr, nxr-nxl |
---|
| 2159 | DO jp = nys, nyn |
---|
| 2160 | number_of_agents = agt_count(jp,ip) |
---|
| 2161 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2162 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2163 | DO n = 1, number_of_agents |
---|
| 2164 | ! |
---|
| 2165 | !-- Apply boundary conditions |
---|
| 2166 | IF ( agents(n)%x < 0.0_wp ) THEN |
---|
| 2167 | |
---|
| 2168 | IF ( ibc_mas_lr == 0 ) THEN |
---|
| 2169 | ! |
---|
| 2170 | !-- Cyclic boundary. Relevant coordinate has to be changed. |
---|
| 2171 | agents(n)%x = ( nx + 1 ) * dx + agents(n)%x |
---|
| 2172 | agents(n)%origin_x = ( nx + 1 ) * dx + & |
---|
| 2173 | agents(n)%origin_x |
---|
| 2174 | ELSEIF ( ibc_mas_lr == 1 ) THEN |
---|
| 2175 | ! |
---|
| 2176 | !-- Agent absorption |
---|
| 2177 | agents(n)%agent_mask = .FALSE. |
---|
| 2178 | deleted_agents = deleted_agents + 1 |
---|
| 2179 | ENDIF |
---|
| 2180 | |
---|
| 2181 | ELSEIF ( agents(n)%x >= ( nx + 1 ) * dx ) THEN |
---|
| 2182 | |
---|
| 2183 | IF ( ibc_mas_lr == 0 ) THEN |
---|
| 2184 | ! |
---|
| 2185 | !-- Cyclic boundary. Relevant coordinate has to be changed. |
---|
| 2186 | agents(n)%x = agents(n)%x - ( nx + 1 ) * dx |
---|
| 2187 | |
---|
| 2188 | ELSEIF ( ibc_mas_lr == 1 ) THEN |
---|
| 2189 | ! |
---|
| 2190 | !-- Agent absorption |
---|
| 2191 | agents(n)%agent_mask = .FALSE. |
---|
| 2192 | deleted_agents = deleted_agents + 1 |
---|
| 2193 | ENDIF |
---|
| 2194 | |
---|
| 2195 | ENDIF |
---|
| 2196 | ENDDO |
---|
| 2197 | ENDDO |
---|
| 2198 | ENDDO |
---|
| 2199 | |
---|
| 2200 | DO ip = nxl, nxr |
---|
| 2201 | DO jp = nys, nyn, nyn-nys |
---|
| 2202 | number_of_agents = agt_count(jp,ip) |
---|
| 2203 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2204 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2205 | DO n = 1, number_of_agents |
---|
| 2206 | |
---|
| 2207 | IF ( agents(n)%y < 0.0_wp ) THEN |
---|
| 2208 | |
---|
| 2209 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 2210 | ! |
---|
| 2211 | !-- Cyclic boundary. Relevant coordinate has to be changed. |
---|
| 2212 | agents(n)%y = ( ny + 1 ) * dy + agents(n)%y |
---|
| 2213 | agents(n)%origin_y = ( ny + 1 ) * dy + & |
---|
| 2214 | agents(n)%origin_y |
---|
| 2215 | |
---|
| 2216 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 2217 | ! |
---|
| 2218 | !-- Agent absorption |
---|
| 2219 | agents(n)%agent_mask = .FALSE. |
---|
| 2220 | deleted_agents = deleted_agents + 1 |
---|
| 2221 | ENDIF |
---|
| 2222 | |
---|
| 2223 | ELSEIF ( agents(n)%y >= ( ny + 0.5_wp ) * dy ) THEN |
---|
| 2224 | |
---|
| 2225 | IF ( ibc_mas_ns == 0 ) THEN |
---|
| 2226 | ! |
---|
| 2227 | !-- Cyclic boundary. Relevant coordinate has to be changed. |
---|
| 2228 | agents(n)%y = agents(n)%y - ( ny + 1 ) * dy |
---|
| 2229 | |
---|
| 2230 | ELSEIF ( ibc_mas_ns == 1 ) THEN |
---|
| 2231 | ! |
---|
| 2232 | !-- Agent absorption |
---|
| 2233 | agents(n)%agent_mask = .FALSE. |
---|
| 2234 | deleted_agents = deleted_agents + 1 |
---|
| 2235 | ENDIF |
---|
| 2236 | |
---|
| 2237 | ENDIF |
---|
| 2238 | |
---|
| 2239 | ENDDO |
---|
| 2240 | ENDDO |
---|
| 2241 | ENDDO |
---|
| 2242 | #endif |
---|
| 2243 | |
---|
| 2244 | ! |
---|
[3201] | 2245 | !-- Accumulate the number of agents transferred between the subdomains) |
---|
[3159] | 2246 | CALL mas_eh_ghost_exchange |
---|
| 2247 | |
---|
| 2248 | END SUBROUTINE mas_eh_exchange_horiz |
---|
| 2249 | |
---|
| 2250 | !------------------------------------------------------------------------------! |
---|
| 2251 | ! Description: |
---|
| 2252 | ! ------------ |
---|
| 2253 | !> Sends the agents from the three gridcells closest to the |
---|
| 2254 | !> north/south/left/right border of a PE to the corresponding neighbors ghost |
---|
| 2255 | !> layer (which is three grid boxes deep) |
---|
| 2256 | !------------------------------------------------------------------------------! |
---|
| 2257 | SUBROUTINE mas_eh_ghost_exchange |
---|
| 2258 | |
---|
| 2259 | IMPLICIT NONE |
---|
[3201] | 2260 | |
---|
| 2261 | #if defined( __parallel ) |
---|
| 2262 | |
---|
[3159] | 2263 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 2264 | INTEGER(iwp) :: jp !< index variable along y |
---|
| 2265 | INTEGER(iwp) :: agt_size !< Bit size of agent datatype |
---|
| 2266 | INTEGER(iwp) :: ghla_count !< ghost points left agent |
---|
| 2267 | INTEGER(iwp) :: ghna_count !< ghost points north agent |
---|
| 2268 | INTEGER(iwp) :: ghra_count !< ghost points right agent |
---|
| 2269 | INTEGER(iwp) :: ghsa_count !< ghost points south agent |
---|
| 2270 | |
---|
| 2271 | LOGICAL :: ghla_empty !< ghost points left agent |
---|
| 2272 | LOGICAL :: ghla_empty_rcv !< ghost points left agent |
---|
| 2273 | LOGICAL :: ghna_empty !< ghost points north agent |
---|
| 2274 | LOGICAL :: ghna_empty_rcv !< ghost points north agent |
---|
| 2275 | LOGICAL :: ghra_empty !< ghost points right agent |
---|
| 2276 | LOGICAL :: ghra_empty_rcv !< ghost points right agent |
---|
| 2277 | LOGICAL :: ghsa_empty !< ghost points south agent |
---|
| 2278 | LOGICAL :: ghsa_empty_rcv !< ghost points south agent |
---|
| 2279 | |
---|
| 2280 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: ghla !< agents received from right PE |
---|
| 2281 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: ghna !< agents received from south PE |
---|
| 2282 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: ghra !< agents received from left PE |
---|
| 2283 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: ghsa !< agents received from north PE |
---|
| 2284 | |
---|
| 2285 | ghla_empty = .TRUE. |
---|
| 2286 | ghna_empty = .TRUE. |
---|
| 2287 | ghra_empty = .TRUE. |
---|
| 2288 | ghsa_empty = .TRUE. |
---|
| 2289 | ! |
---|
| 2290 | !-- reset ghost layer |
---|
| 2291 | DO ip = nxlg, nxl-1 |
---|
| 2292 | DO jp = nysg, nyng |
---|
| 2293 | agt_count(jp,ip) = 0 |
---|
| 2294 | ENDDO |
---|
| 2295 | ENDDO |
---|
| 2296 | DO ip = nxr+1, nxrg |
---|
| 2297 | DO jp = nysg, nyng |
---|
| 2298 | agt_count(jp,ip) = 0 |
---|
| 2299 | ENDDO |
---|
| 2300 | ENDDO |
---|
| 2301 | DO ip = nxl, nxr |
---|
| 2302 | DO jp = nysg, nys-1 |
---|
| 2303 | agt_count(jp,ip) = 0 |
---|
| 2304 | ENDDO |
---|
| 2305 | ENDDO |
---|
| 2306 | DO ip = nxl, nxr |
---|
| 2307 | DO jp = nyn+1, nyng |
---|
| 2308 | agt_count(jp,ip) = 0 |
---|
| 2309 | ENDDO |
---|
| 2310 | ENDDO |
---|
| 2311 | ! |
---|
| 2312 | !-- Transfer of agents from left to right and vice versa |
---|
| 2313 | IF ( pdims(1) /= 1 ) THEN |
---|
| 2314 | ! |
---|
| 2315 | !-- Reset left and right ghost layers |
---|
| 2316 | ghla_count = 0 |
---|
| 2317 | ghra_count = 0 |
---|
| 2318 | ! |
---|
| 2319 | !-- First calculate the storage necessary for sending |
---|
| 2320 | !-- and receiving the data. |
---|
| 2321 | ghla_count = SUM(agt_count(nys:nyn,nxl:nxl+2)) |
---|
| 2322 | ghra_count = SUM(agt_count(nys:nyn,nxr-2:nxr)) |
---|
| 2323 | ! |
---|
| 2324 | !-- No cyclic boundaries for agents |
---|
| 2325 | IF ( nxl == 0 .OR. ghla_count == 0 ) THEN |
---|
| 2326 | ghla_count = 1 |
---|
| 2327 | ELSE |
---|
| 2328 | ghla_empty = .FALSE. |
---|
| 2329 | ENDIF |
---|
| 2330 | IF ( nxr == nx .OR. ghra_count == 0 ) THEN |
---|
| 2331 | ghra_count = 1 |
---|
| 2332 | ELSE |
---|
| 2333 | ghra_empty = .FALSE. |
---|
| 2334 | ENDIF |
---|
| 2335 | ALLOCATE( ghla(1:ghla_count), ghra(1:ghra_count) ) |
---|
| 2336 | ghla = zero_agent |
---|
| 2337 | ghra = zero_agent |
---|
| 2338 | ! |
---|
| 2339 | !-- Get all agents that will be sent left into one array |
---|
| 2340 | ghla_count = 0 |
---|
| 2341 | IF ( nxl /= 0 ) THEN |
---|
| 2342 | DO ip = nxl, nxl+2 |
---|
| 2343 | DO jp = nys, nyn |
---|
| 2344 | |
---|
| 2345 | number_of_agents = agt_count(jp,ip) |
---|
| 2346 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2347 | ghla(ghla_count+1:ghla_count+number_of_agents) & |
---|
| 2348 | = grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2349 | ghla_count = ghla_count + number_of_agents |
---|
| 2350 | |
---|
| 2351 | ENDDO |
---|
| 2352 | ENDDO |
---|
| 2353 | ENDIF |
---|
| 2354 | IF ( ghla_count == 0 ) ghla_count = 1 |
---|
| 2355 | ! |
---|
| 2356 | !-- Get all agents that will be sent right into one array |
---|
| 2357 | ghra_count = 0 |
---|
| 2358 | IF ( nxr /= nx ) THEN |
---|
| 2359 | DO ip = nxr-2, nxr |
---|
| 2360 | DO jp = nys, nyn |
---|
| 2361 | |
---|
| 2362 | number_of_agents = agt_count(jp,ip) |
---|
| 2363 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2364 | ghra(ghra_count+1:ghra_count+number_of_agents) & |
---|
| 2365 | = grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2366 | ghra_count = ghra_count + number_of_agents |
---|
| 2367 | |
---|
| 2368 | ENDDO |
---|
| 2369 | ENDDO |
---|
| 2370 | ENDIF |
---|
| 2371 | IF ( ghra_count == 0 ) ghra_count = 1 |
---|
| 2372 | ! |
---|
| 2373 | !-- Send/receive number of agents that |
---|
| 2374 | !-- will be transferred to/from left/right neighbor |
---|
| 2375 | CALL MPI_SENDRECV( ghla_count, 1, MPI_INTEGER, pleft, 0, & |
---|
| 2376 | ghra_count_recv, 1, MPI_INTEGER, pright, 0, & |
---|
| 2377 | comm2d, status, ierr ) |
---|
| 2378 | ALLOCATE ( agt_gh_r(1:ghra_count_recv) ) |
---|
| 2379 | ! |
---|
| 2380 | !-- Send/receive number of agents that |
---|
| 2381 | !-- will be transferred to/from right/left neighbor |
---|
| 2382 | CALL MPI_SENDRECV( ghra_count, 1, MPI_INTEGER, pright, 0, & |
---|
| 2383 | ghla_count_recv, 1, MPI_INTEGER, pleft, 0, & |
---|
| 2384 | comm2d, status, ierr ) |
---|
| 2385 | ! |
---|
| 2386 | !-- Send/receive flag that indicates if there are actually any agents |
---|
| 2387 | !-- in ghost layer |
---|
| 2388 | CALL MPI_SENDRECV( ghla_empty, 1, MPI_LOGICAL, pleft, 1, & |
---|
| 2389 | ghra_empty_rcv, 1, MPI_LOGICAL, pright,1, & |
---|
| 2390 | comm2d, status, ierr ) |
---|
| 2391 | CALL MPI_SENDRECV( ghra_empty, 1, MPI_LOGICAL, pright,1, & |
---|
| 2392 | ghla_empty_rcv, 1, MPI_LOGICAL, pleft, 1, & |
---|
| 2393 | comm2d, status, ierr ) |
---|
| 2394 | |
---|
| 2395 | |
---|
| 2396 | ALLOCATE ( agt_gh_l(1:ghla_count_recv) ) |
---|
| 2397 | ! |
---|
| 2398 | !-- Get bit size of one agent |
---|
[3160] | 2399 | agt_size = STORAGE_SIZE(zero_agent)/8 |
---|
[3159] | 2400 | ! |
---|
| 2401 | !-- Send/receive agents to/from left/right neighbor |
---|
| 2402 | CALL MPI_SENDRECV( ghla, ghla_count * agt_size, MPI_BYTE, & |
---|
| 2403 | pleft, 1, & |
---|
| 2404 | agt_gh_r, ghra_count_recv * agt_size, MPI_BYTE, & |
---|
| 2405 | pright,1, & |
---|
| 2406 | comm2d, status, ierr ) |
---|
| 2407 | ! |
---|
| 2408 | !-- Send/receive agents to/from left/right neighbor |
---|
| 2409 | CALL MPI_SENDRECV( ghra, ghra_count * agt_size, MPI_BYTE, & |
---|
| 2410 | pright,1, & |
---|
| 2411 | agt_gh_l, ghla_count_recv * agt_size, MPI_BYTE, & |
---|
| 2412 | pleft, 1, & |
---|
| 2413 | comm2d, status, ierr ) |
---|
| 2414 | ! |
---|
| 2415 | !-- If agents were received, add them to the respective ghost layer cells |
---|
| 2416 | IF ( .NOT. ghra_empty_rcv ) THEN |
---|
| 2417 | CALL mas_eh_add_ghost_agents_to_gridcell(agt_gh_r) |
---|
| 2418 | ENDIF |
---|
| 2419 | |
---|
| 2420 | IF ( .NOT. ghla_empty_rcv ) THEN |
---|
| 2421 | CALL mas_eh_add_ghost_agents_to_gridcell(agt_gh_l) |
---|
| 2422 | ENDIF |
---|
| 2423 | |
---|
| 2424 | DEALLOCATE( ghla, ghra, agt_gh_l, agt_gh_r ) |
---|
| 2425 | |
---|
| 2426 | ENDIF |
---|
| 2427 | |
---|
| 2428 | ! |
---|
| 2429 | !-- Transfer of agents from south to north and vice versa |
---|
| 2430 | IF ( pdims(2) /= 1 ) THEN |
---|
| 2431 | ! |
---|
| 2432 | !-- Reset south and north ghost layers |
---|
| 2433 | ghsa_count = 0 |
---|
| 2434 | ghna_count = 0 |
---|
| 2435 | ! |
---|
| 2436 | !-- First calculate the storage necessary for sending |
---|
| 2437 | !-- and receiving the data. |
---|
| 2438 | ghsa_count = SUM(agt_count(nys:nys+2,nxlg:nxrg)) |
---|
| 2439 | ghna_count = SUM(agt_count(nyn-2:nyn,nxlg:nxrg)) |
---|
| 2440 | ! |
---|
| 2441 | !-- No cyclic boundaries for agents |
---|
| 2442 | IF ( nys == 0 .OR. ghsa_count == 0 ) THEN |
---|
| 2443 | ghsa_count = 1 |
---|
| 2444 | ELSE |
---|
| 2445 | ghsa_empty = .FALSE. |
---|
| 2446 | ENDIF |
---|
| 2447 | IF ( nyn == ny .OR. ghna_count == 0 ) THEN |
---|
| 2448 | ghna_count = 1 |
---|
| 2449 | ELSE |
---|
| 2450 | ghna_empty = .FALSE. |
---|
| 2451 | ENDIF |
---|
| 2452 | ALLOCATE( ghsa(1:ghsa_count), ghna(1:ghna_count) ) |
---|
| 2453 | ghsa = zero_agent |
---|
| 2454 | ghna = zero_agent |
---|
| 2455 | ! |
---|
| 2456 | !-- Get all agents that will be sent south into one array |
---|
| 2457 | ghsa_count = 0 |
---|
| 2458 | IF ( nys /= 0 ) THEN |
---|
| 2459 | DO ip = nxlg, nxrg |
---|
| 2460 | DO jp = nys, nys+2 |
---|
| 2461 | |
---|
| 2462 | number_of_agents = agt_count(jp,ip) |
---|
| 2463 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2464 | ghsa(ghsa_count+1:ghsa_count+number_of_agents) & |
---|
| 2465 | = grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2466 | ghsa_count = ghsa_count + number_of_agents |
---|
| 2467 | |
---|
| 2468 | ENDDO |
---|
| 2469 | ENDDO |
---|
| 2470 | ENDIF |
---|
| 2471 | IF ( ghsa_count == 0 ) ghsa_count = 1 |
---|
| 2472 | ! |
---|
| 2473 | !-- Get all agents that will be sent north into one array |
---|
| 2474 | ghna_count = 0 |
---|
| 2475 | IF ( nyn /= ny ) THEN |
---|
| 2476 | DO ip = nxlg, nxrg |
---|
| 2477 | DO jp = nyn-2, nyn |
---|
| 2478 | |
---|
| 2479 | number_of_agents = agt_count(jp,ip) |
---|
| 2480 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2481 | ghna(ghna_count+1:ghna_count+number_of_agents) & |
---|
| 2482 | = grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 2483 | ghna_count = ghna_count + number_of_agents |
---|
| 2484 | |
---|
| 2485 | ENDDO |
---|
| 2486 | ENDDO |
---|
| 2487 | ENDIF |
---|
| 2488 | IF ( ghna_count == 0 ) ghna_count = 1 |
---|
| 2489 | ! |
---|
| 2490 | !-- Send/receive number of agents that |
---|
| 2491 | !-- will be transferred to/from south/north neighbor |
---|
| 2492 | CALL MPI_SENDRECV( ghsa_count, 1, MPI_INTEGER, psouth, 0, & |
---|
| 2493 | ghna_count_recv, 1, MPI_INTEGER, pnorth, 0, & |
---|
| 2494 | comm2d, status, ierr ) |
---|
| 2495 | ALLOCATE ( agt_gh_n(1:ghna_count_recv) ) |
---|
| 2496 | ! |
---|
| 2497 | !-- Send/receive number of agents that |
---|
| 2498 | !-- will be transferred to/from north/south neighbor |
---|
| 2499 | CALL MPI_SENDRECV( ghna_count, 1, MPI_INTEGER, pnorth, 0, & |
---|
| 2500 | ghsa_count_recv, 1, MPI_INTEGER, psouth, 0, & |
---|
| 2501 | comm2d, status, ierr ) |
---|
| 2502 | ! |
---|
| 2503 | !-- Send/receive flag that indicates if there are actually any agents |
---|
| 2504 | !-- in ghost layer |
---|
| 2505 | CALL MPI_SENDRECV( ghsa_empty, 1, MPI_LOGICAL, psouth, 1, & |
---|
| 2506 | ghna_empty_rcv, 1, MPI_LOGICAL, pnorth, 1, & |
---|
| 2507 | comm2d, status, ierr ) |
---|
| 2508 | CALL MPI_SENDRECV( ghna_empty, 1, MPI_LOGICAL, pnorth, 1, & |
---|
| 2509 | ghsa_empty_rcv, 1, MPI_LOGICAL, psouth, 1, & |
---|
| 2510 | comm2d, status, ierr ) |
---|
| 2511 | |
---|
| 2512 | |
---|
| 2513 | ALLOCATE ( agt_gh_s(1:ghsa_count_recv) ) |
---|
| 2514 | ! |
---|
| 2515 | !-- Get bit size of one agent |
---|
[3160] | 2516 | agt_size = STORAGE_SIZE(zero_agent)/8 |
---|
[3159] | 2517 | ! |
---|
| 2518 | !-- Send/receive agents to/from south/north neighbor |
---|
| 2519 | CALL MPI_SENDRECV( ghsa, ghsa_count * agt_size, MPI_BYTE, & |
---|
| 2520 | psouth,1, & |
---|
| 2521 | agt_gh_n, ghna_count_recv * agt_size, MPI_BYTE, & |
---|
| 2522 | pnorth,1, & |
---|
| 2523 | comm2d, status, ierr ) |
---|
| 2524 | ! |
---|
| 2525 | !-- Send/receive agents to/from south/north neighbor |
---|
| 2526 | CALL MPI_SENDRECV( ghna, ghna_count * agt_size, MPI_BYTE, & |
---|
| 2527 | pnorth,1, & |
---|
| 2528 | agt_gh_s, ghsa_count_recv * agt_size, MPI_BYTE, & |
---|
| 2529 | psouth,1, & |
---|
| 2530 | comm2d, status, ierr ) |
---|
| 2531 | ! |
---|
| 2532 | !-- If agents were received, add them to the respective ghost layer cells |
---|
| 2533 | IF ( .NOT. ghna_empty_rcv ) THEN |
---|
| 2534 | CALL mas_eh_add_ghost_agents_to_gridcell(agt_gh_n) |
---|
| 2535 | ENDIF |
---|
| 2536 | |
---|
| 2537 | IF ( .NOT. ghsa_empty_rcv ) THEN |
---|
| 2538 | CALL mas_eh_add_ghost_agents_to_gridcell(agt_gh_s) |
---|
| 2539 | ENDIF |
---|
| 2540 | |
---|
| 2541 | DEALLOCATE( ghna, ghsa, agt_gh_n, agt_gh_s ) |
---|
| 2542 | |
---|
| 2543 | ENDIF |
---|
| 2544 | |
---|
[3201] | 2545 | #endif |
---|
| 2546 | |
---|
[3159] | 2547 | END SUBROUTINE mas_eh_ghost_exchange |
---|
| 2548 | |
---|
| 2549 | !------------------------------------------------------------------------------! |
---|
| 2550 | ! Description: |
---|
| 2551 | ! ------------ |
---|
| 2552 | !> If an agent moves from one grid cell to another (on the current |
---|
| 2553 | !> processor!), this subroutine moves the corresponding element from the |
---|
| 2554 | !> agent array of the old grid cell to the agent array of the new grid |
---|
| 2555 | !> cell. |
---|
| 2556 | !------------------------------------------------------------------------------! |
---|
| 2557 | SUBROUTINE mas_eh_move_agent |
---|
| 2558 | |
---|
| 2559 | IMPLICIT NONE |
---|
| 2560 | |
---|
| 2561 | INTEGER(iwp) :: i !< grid index (x) of agent position |
---|
| 2562 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 2563 | INTEGER(iwp) :: j !< grid index (y) of agent position |
---|
| 2564 | INTEGER(iwp) :: jp !< index variable along y |
---|
| 2565 | INTEGER(iwp) :: n !< index variable for agent array |
---|
| 2566 | INTEGER(iwp) :: na_before_move !< number of agents per grid box before moving |
---|
| 2567 | INTEGER(iwp) :: aindex !< dummy argument for number of new agent per grid box |
---|
| 2568 | |
---|
| 2569 | TYPE(agent_type), DIMENSION(:), POINTER :: agents_before_move !< agents before moving |
---|
| 2570 | |
---|
| 2571 | DO ip = nxl, nxr |
---|
| 2572 | DO jp = nys, nyn |
---|
| 2573 | |
---|
| 2574 | na_before_move = agt_count(jp,ip) |
---|
| 2575 | IF ( na_before_move <= 0 ) CYCLE |
---|
| 2576 | agents_before_move => grid_agents(jp,ip)%agents(1:na_before_move) |
---|
| 2577 | |
---|
| 2578 | DO n = 1, na_before_move |
---|
| 2579 | i = agents_before_move(n)%x * ddx |
---|
| 2580 | j = agents_before_move(n)%y * ddy |
---|
| 2581 | |
---|
| 2582 | !-- For mas_eh_exchange_horiz to work properly agents need to be |
---|
| 2583 | !-- moved to the outermost gridboxes of the respective processor. |
---|
| 2584 | !-- If the agent index is inside the processor the following |
---|
| 2585 | !-- lines will not change the index |
---|
| 2586 | i = MIN ( i , nxr ) |
---|
| 2587 | i = MAX ( i , nxl ) |
---|
| 2588 | j = MIN ( j , nyn ) |
---|
| 2589 | j = MAX ( j , nys ) |
---|
| 2590 | |
---|
| 2591 | ! |
---|
| 2592 | !-- Check if agent has moved to another grid cell. |
---|
| 2593 | IF ( i /= ip .OR. j /= jp ) THEN |
---|
| 2594 | ! |
---|
| 2595 | !-- If the agent stays on the same processor, the agent |
---|
| 2596 | !-- will be added to the agent array of the new processor. |
---|
| 2597 | number_of_agents = agt_count(j,i) |
---|
| 2598 | agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 2599 | |
---|
| 2600 | aindex = number_of_agents+1 |
---|
| 2601 | IF ( aindex > SIZE(grid_agents(j,i)%agents) ) & |
---|
| 2602 | THEN |
---|
| 2603 | CALL mas_eh_realloc_agents_array(i,j) |
---|
| 2604 | ENDIF |
---|
| 2605 | |
---|
| 2606 | grid_agents(j,i)%agents(aindex) = agents_before_move(n) |
---|
| 2607 | agt_count(j,i) = aindex |
---|
| 2608 | |
---|
| 2609 | agents_before_move(n)%agent_mask = .FALSE. |
---|
| 2610 | ENDIF |
---|
| 2611 | ENDDO |
---|
| 2612 | |
---|
| 2613 | ENDDO |
---|
| 2614 | ENDDO |
---|
| 2615 | |
---|
| 2616 | RETURN |
---|
| 2617 | |
---|
| 2618 | END SUBROUTINE mas_eh_move_agent |
---|
| 2619 | |
---|
| 2620 | !------------------------------------------------------------------------------! |
---|
| 2621 | ! Description: |
---|
| 2622 | ! ------------ |
---|
| 2623 | !> If the allocated memory for the agent array do not suffice to add arriving |
---|
| 2624 | !> agents from neighbour grid cells, this subrouting reallocates the |
---|
| 2625 | !> agent array to assure enough memory is available. |
---|
| 2626 | !------------------------------------------------------------------------------! |
---|
| 2627 | SUBROUTINE mas_eh_realloc_agents_array (i,j,size_in) |
---|
| 2628 | |
---|
| 2629 | IMPLICIT NONE |
---|
| 2630 | |
---|
| 2631 | INTEGER(iwp) :: old_size !< old array size |
---|
| 2632 | INTEGER(iwp) :: new_size !< new array size |
---|
| 2633 | |
---|
| 2634 | INTEGER(iwp), INTENT(in) :: i !< grid index (y) |
---|
| 2635 | INTEGER(iwp), INTENT(in) :: j !< grid index (y) |
---|
| 2636 | |
---|
| 2637 | INTEGER(iwp), INTENT(in), OPTIONAL :: size_in !< size of input array |
---|
| 2638 | |
---|
| 2639 | TYPE(agent_type), DIMENSION(10) :: tmp_agents_s !< temporary static agent array |
---|
| 2640 | |
---|
| 2641 | TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: tmp_agents_d !< temporary dynamic agent array |
---|
| 2642 | |
---|
| 2643 | old_size = SIZE(grid_agents(j,i)%agents) |
---|
| 2644 | |
---|
| 2645 | IF ( PRESENT(size_in) ) THEN |
---|
| 2646 | new_size = size_in |
---|
| 2647 | ELSE |
---|
| 2648 | new_size = old_size * ( 1.0_wp + alloc_factor_mas / 100.0_wp ) |
---|
| 2649 | ENDIF |
---|
| 2650 | |
---|
| 2651 | new_size = MAX( new_size, min_nr_agent, old_size + 1 ) |
---|
| 2652 | |
---|
| 2653 | IF ( old_size <= 10 ) THEN |
---|
| 2654 | |
---|
| 2655 | tmp_agents_s(1:old_size) = grid_agents(j,i)%agents(1:old_size) |
---|
| 2656 | |
---|
| 2657 | DEALLOCATE(grid_agents(j,i)%agents) |
---|
| 2658 | ALLOCATE(grid_agents(j,i)%agents(new_size)) |
---|
| 2659 | |
---|
| 2660 | grid_agents(j,i)%agents(1:old_size) = tmp_agents_s(1:old_size) |
---|
| 2661 | grid_agents(j,i)%agents(old_size+1:new_size) = zero_agent |
---|
| 2662 | |
---|
| 2663 | ELSE |
---|
| 2664 | |
---|
| 2665 | ALLOCATE(tmp_agents_d(new_size)) |
---|
| 2666 | tmp_agents_d(1:old_size) = grid_agents(j,i)%agents |
---|
| 2667 | |
---|
| 2668 | DEALLOCATE(grid_agents(j,i)%agents) |
---|
| 2669 | ALLOCATE(grid_agents(j,i)%agents(new_size)) |
---|
| 2670 | |
---|
| 2671 | grid_agents(j,i)%agents(1:old_size) = tmp_agents_d(1:old_size) |
---|
| 2672 | grid_agents(j,i)%agents(old_size+1:new_size) = zero_agent |
---|
| 2673 | |
---|
| 2674 | DEALLOCATE(tmp_agents_d) |
---|
| 2675 | |
---|
| 2676 | ENDIF |
---|
| 2677 | agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 2678 | |
---|
| 2679 | RETURN |
---|
| 2680 | END SUBROUTINE mas_eh_realloc_agents_array |
---|
| 2681 | |
---|
| 2682 | !------------------------------------------------------------------------------! |
---|
| 2683 | ! Description: |
---|
| 2684 | ! ------------ |
---|
| 2685 | !> Inquires prognostic model quantities at the position of each agent and |
---|
| 2686 | !> stores them in that agent for later output |
---|
| 2687 | !------------------------------------------------------------------------------! |
---|
| 2688 | SUBROUTINE mas_get_prognostic_quantities |
---|
| 2689 | |
---|
| 2690 | USE arrays_3d, & |
---|
| 2691 | ONLY: u, v, pt, hyp |
---|
| 2692 | |
---|
[3241] | 2693 | ! USE chemistry_model_mod, & |
---|
| 2694 | ! ONLY: chem_species |
---|
[3159] | 2695 | |
---|
[3241] | 2696 | ! USE control_parameters, & |
---|
| 2697 | ! ONLY: air_chemistry |
---|
[3159] | 2698 | |
---|
| 2699 | IMPLICIT NONE |
---|
| 2700 | |
---|
| 2701 | INTEGER(iwp) :: i_offset !< index offset for windspeed measurement |
---|
| 2702 | INTEGER(iwp) :: il !< x-index |
---|
| 2703 | INTEGER(iwp) :: is !< subgrid box counter |
---|
| 2704 | INTEGER(iwp) :: j_offset !< index offset for windspeed measurement |
---|
| 2705 | INTEGER(iwp) :: jl !< y-index |
---|
| 2706 | INTEGER(iwp) :: kl !< z-index |
---|
| 2707 | INTEGER(iwp) :: nl !< agent counter |
---|
| 2708 | INTEGER(iwp) :: se !< subgrid box end index |
---|
| 2709 | INTEGER(iwp) :: si !< subgrid box start index |
---|
| 2710 | |
---|
| 2711 | REAL(wp) :: u_a !< windspeed at agent position (x) |
---|
| 2712 | REAL(wp) :: v_a !< windspeed at agent position (y) |
---|
| 2713 | |
---|
| 2714 | DO il = nxl, nxr |
---|
| 2715 | DO jl = nys, nyn |
---|
| 2716 | |
---|
| 2717 | number_of_agents = agt_count(jl,il) |
---|
| 2718 | ! |
---|
| 2719 | !-- If grid cell is empty, cycle |
---|
| 2720 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 2721 | kl = s_measure_height(jl,il) |
---|
| 2722 | |
---|
| 2723 | agents => grid_agents(jl,il)%agents(1:number_of_agents) |
---|
| 2724 | ! |
---|
| 2725 | !-- loop over the four subgrid boxes |
---|
| 2726 | DO is = 0,3 |
---|
| 2727 | ! |
---|
| 2728 | !-- Set indices |
---|
| 2729 | si = grid_agents(jl,il)%start_index(is) |
---|
| 2730 | se = grid_agents(jl,il)%end_index(is) |
---|
| 2731 | DO nl = si, se |
---|
| 2732 | ! |
---|
| 2733 | !-- Calculate index offset in x-direction: |
---|
| 2734 | !-- Left value if wall right of grid box |
---|
| 2735 | !-- Right value if wall left of grid box |
---|
| 2736 | !-- Else the one that is closer to the agent |
---|
| 2737 | IF ( BTEST( obstacle_flags( jl, il+1 ), 6 ) ) THEN |
---|
| 2738 | i_offset = 0 |
---|
| 2739 | ELSEIF ( BTEST( obstacle_flags( jl, il-1 ), 2 ) ) THEN |
---|
| 2740 | i_offset = 1 |
---|
| 2741 | ELSE |
---|
| 2742 | i_offset = MERGE( 0, 1, BTEST(is,1) ) |
---|
| 2743 | ENDIF |
---|
| 2744 | u_a = u( kl, jl, il + i_offset ) |
---|
| 2745 | ! |
---|
| 2746 | !-- Calculate index offset in y-direction: |
---|
| 2747 | !-- South value if wall north of grid box |
---|
| 2748 | !-- North value if wall south of grid box |
---|
| 2749 | !-- Else the one that is closer to the agent |
---|
| 2750 | IF ( BTEST( obstacle_flags( jl+1, il ), 4 ) ) THEN |
---|
| 2751 | j_offset = 0 |
---|
| 2752 | ELSEIF ( BTEST( obstacle_flags( jl-1, il ), 0 ) ) THEN |
---|
| 2753 | j_offset = 1 |
---|
| 2754 | ELSE |
---|
| 2755 | j_offset = MERGE( 0, 1, BTEST(is,0) ) |
---|
| 2756 | ENDIF |
---|
| 2757 | v_a = v( kl, jl + j_offset, il ) |
---|
| 2758 | ! |
---|
| 2759 | !-- Calculate windspeed at agent postion |
---|
| 2760 | agents(nl)%windspeed = SQRT(u_a**2 + v_a**2) |
---|
| 2761 | ! |
---|
| 2762 | !-- Calculate temperature at agent position |
---|
| 2763 | agents(nl)%t = pt(kl,jl,il) * & |
---|
| 2764 | ( hyp(kl) / 100000.0_wp )**0.286_wp |
---|
| 2765 | ! ! |
---|
| 2766 | ! !-- Get PM10 concentration at agent position, if possible |
---|
| 2767 | ! IF ( ind_pm10 == -9 ) THEN |
---|
| 2768 | ! agents(nl)%pm10 = 0.0_wp |
---|
| 2769 | ! ELSE |
---|
| 2770 | ! agents(nl)%pm10 = chem_species(ind_pm10)%conc(kl,jl,il) |
---|
| 2771 | ! ENDIF |
---|
| 2772 | ! ! |
---|
| 2773 | ! !-- Get PM10 concentration at agent position, if possible |
---|
| 2774 | ! IF ( ind_pm25 == -9 ) THEN |
---|
| 2775 | ! agents(nl)%pm25 = 0.0_wp |
---|
| 2776 | ! ELSE |
---|
| 2777 | ! agents(nl)%pm25 = chem_species(ind_pm25)%conc(kl,jl,il) |
---|
| 2778 | ! ENDIF |
---|
| 2779 | ENDDO |
---|
| 2780 | |
---|
| 2781 | ENDDO |
---|
| 2782 | |
---|
| 2783 | ENDDO |
---|
| 2784 | ENDDO |
---|
| 2785 | |
---|
| 2786 | END SUBROUTINE mas_get_prognostic_quantities |
---|
| 2787 | |
---|
| 2788 | !------------------------------------------------------------------------------! |
---|
| 2789 | ! Description: |
---|
| 2790 | ! ------------ |
---|
| 2791 | !> Adds an item to the priority queue (binary heap) at the correct position |
---|
| 2792 | !------------------------------------------------------------------------------! |
---|
| 2793 | SUBROUTINE mas_heap_insert_item( id, priority ) |
---|
| 2794 | |
---|
| 2795 | IMPLICIT NONE |
---|
| 2796 | |
---|
| 2797 | INTEGER(iwp) :: cur_pos !< current position |
---|
| 2798 | INTEGER(iwp) :: id !< mesh ID of item |
---|
| 2799 | |
---|
| 2800 | REAL(wp) :: priority !< item priority |
---|
| 2801 | |
---|
| 2802 | TYPE(heap_item) :: item !< heap item |
---|
| 2803 | |
---|
| 2804 | item%mesh_id = id |
---|
| 2805 | item%priority = priority |
---|
| 2806 | ! |
---|
| 2807 | !-- Extend heap, if necessary |
---|
| 2808 | IF ( heap_count + 1 > SIZE(queue) ) THEN |
---|
| 2809 | CALL mas_heap_extend |
---|
| 2810 | ENDIF |
---|
| 2811 | ! |
---|
| 2812 | !-- Insert item at first unoccupied postion (highest index) of heap |
---|
| 2813 | cur_pos = heap_count |
---|
| 2814 | queue(cur_pos) = item |
---|
| 2815 | ! |
---|
| 2816 | !-- Sort while inserted item is not at top of heap |
---|
| 2817 | DO WHILE ( cur_pos /= 0 ) |
---|
| 2818 | ! |
---|
| 2819 | !-- If priority < its parent's priority, swap them. |
---|
| 2820 | !-- Else, sorting is done. |
---|
| 2821 | IF ( queue(cur_pos)%priority & |
---|
| 2822 | < queue(FLOOR((cur_pos)/2.))%priority ) & |
---|
| 2823 | THEN |
---|
| 2824 | item = queue(cur_pos) |
---|
| 2825 | queue(cur_pos) = queue(FLOOR((cur_pos)/2.)) |
---|
| 2826 | queue(FLOOR((cur_pos)/2.)) = item |
---|
| 2827 | cur_pos = FLOOR((cur_pos)/2.) |
---|
| 2828 | ELSE |
---|
| 2829 | EXIT |
---|
| 2830 | ENDIF |
---|
| 2831 | ENDDO |
---|
| 2832 | ! |
---|
| 2833 | !-- Item was added to heap, so the heap count increases |
---|
| 2834 | heap_count = heap_count + 1 |
---|
| 2835 | |
---|
| 2836 | END SUBROUTINE mas_heap_insert_item |
---|
| 2837 | |
---|
| 2838 | !------------------------------------------------------------------------------! |
---|
| 2839 | ! Description: |
---|
| 2840 | ! ------------ |
---|
| 2841 | !> Extends the size of the priority queue (binary heap) |
---|
| 2842 | !------------------------------------------------------------------------------! |
---|
| 2843 | SUBROUTINE mas_heap_extend |
---|
| 2844 | |
---|
| 2845 | IMPLICIT NONE |
---|
| 2846 | |
---|
| 2847 | INTEGER(iwp) :: soh !< size of heap |
---|
| 2848 | |
---|
| 2849 | TYPE(heap_item), DIMENSION(:), ALLOCATABLE :: dummy_heap !< dummy heap |
---|
| 2850 | |
---|
| 2851 | soh = SIZE(queue)-1 |
---|
| 2852 | ALLOCATE(dummy_heap(0:soh)) |
---|
| 2853 | dummy_heap = queue |
---|
| 2854 | DEALLOCATE(queue) |
---|
| 2855 | ALLOCATE(queue(0:2*soh+1)) |
---|
| 2856 | queue(0:soh) = dummy_heap(0:soh) |
---|
| 2857 | |
---|
| 2858 | END SUBROUTINE mas_heap_extend |
---|
| 2859 | |
---|
| 2860 | !------------------------------------------------------------------------------! |
---|
| 2861 | ! Description: |
---|
| 2862 | ! ------------ |
---|
| 2863 | !> Removes first (smallest) element from the priority queue, reorders the rest |
---|
| 2864 | !> and returns the ID of the removed mesh point |
---|
| 2865 | !------------------------------------------------------------------------------! |
---|
| 2866 | SUBROUTINE mas_heap_extract_item ( id ) |
---|
| 2867 | |
---|
| 2868 | IMPLICIT NONE |
---|
| 2869 | |
---|
| 2870 | INTEGER(iwp) :: id !< ID of item extracted item |
---|
| 2871 | INTEGER(iwp) :: child !< child of item in heap |
---|
| 2872 | INTEGER(iwp) :: cur_pos !< current position of item in heap |
---|
| 2873 | |
---|
| 2874 | TYPE(heap_item) :: dummy |
---|
| 2875 | ! |
---|
| 2876 | !-- Get ID of mesh point with lowest priority (extracted item: top of heap) |
---|
| 2877 | id = queue(0)%mesh_id |
---|
| 2878 | ! |
---|
| 2879 | !-- Put last item in heap at first position |
---|
| 2880 | queue(0) = queue(heap_count-1) |
---|
| 2881 | cur_pos = 0 |
---|
| 2882 | DO |
---|
| 2883 | ! |
---|
| 2884 | !-- If current item has no children, sorting is done |
---|
| 2885 | IF( 2*cur_pos+1 > heap_count - 1 ) THEN |
---|
| 2886 | EXIT |
---|
| 2887 | ! |
---|
| 2888 | !-- If current item has only one child, check if item and its child are |
---|
| 2889 | !-- ordered correctly. Else, swap them. |
---|
| 2890 | ELSEIF ( 2*cur_pos+2 > heap_count - 1 ) THEN |
---|
| 2891 | IF ( queue(cur_pos)%priority > queue(2*cur_pos+1)%priority ) THEN |
---|
| 2892 | dummy = queue(cur_pos) |
---|
| 2893 | queue(cur_pos) = queue(2*cur_pos+1) |
---|
| 2894 | queue(2*cur_pos+1) = dummy |
---|
| 2895 | cur_pos = 2*cur_pos+1 |
---|
| 2896 | ELSE |
---|
| 2897 | EXIT |
---|
| 2898 | ENDIF |
---|
| 2899 | ELSE |
---|
| 2900 | ! |
---|
| 2901 | !-- determine the smaller child |
---|
| 2902 | IF ( queue(2*cur_pos+1)%priority & |
---|
| 2903 | >= queue(2*cur_pos+2)%priority ) & |
---|
| 2904 | THEN |
---|
| 2905 | child = 2 |
---|
| 2906 | ELSE |
---|
| 2907 | child = 1 |
---|
| 2908 | ENDIF |
---|
| 2909 | ! |
---|
| 2910 | !-- Check if item and its smaller child are ordered falsely. If so, |
---|
| 2911 | !-- swap them. Else, sorting is done. |
---|
| 2912 | IF ( queue(cur_pos)%priority > queue(2*cur_pos+child )%priority ) & |
---|
| 2913 | THEN |
---|
| 2914 | dummy = queue(cur_pos) |
---|
| 2915 | queue(cur_pos) = queue(2*cur_pos+child) |
---|
| 2916 | queue(2*cur_pos+child) = dummy |
---|
| 2917 | cur_pos = 2*cur_pos+child |
---|
| 2918 | ELSE |
---|
| 2919 | EXIT |
---|
| 2920 | ENDIF |
---|
| 2921 | ENDIF |
---|
| 2922 | ENDDO |
---|
| 2923 | ! |
---|
| 2924 | !-- Top item was removed from heap, thus, heap_cout decreases by one |
---|
| 2925 | heap_count = heap_count-1 |
---|
| 2926 | |
---|
| 2927 | END SUBROUTINE mas_heap_extract_item |
---|
| 2928 | |
---|
| 2929 | !------------------------------------------------------------------------------! |
---|
| 2930 | ! Description: |
---|
| 2931 | ! ------------ |
---|
| 2932 | !> Initialization of Multi Agent System |
---|
| 2933 | !------------------------------------------------------------------------------! |
---|
| 2934 | SUBROUTINE mas_init |
---|
| 2935 | |
---|
[3241] | 2936 | ! USE chem_gasphase_mod, & |
---|
| 2937 | ! ONLY: nspec |
---|
[3159] | 2938 | |
---|
[3241] | 2939 | ! USE chemistry_model_mod, & |
---|
| 2940 | ! ONLY: chem_species |
---|
[3159] | 2941 | |
---|
| 2942 | USE control_parameters, & |
---|
[3241] | 2943 | ONLY: coupling_char, initializing_actions, io_blocks, io_group |
---|
[3159] | 2944 | |
---|
| 2945 | USE surface_mod, & |
---|
| 2946 | ONLY: get_topography_top_index, get_topography_top_index_ji |
---|
| 2947 | |
---|
| 2948 | USE arrays_3d, & |
---|
| 2949 | ONLY: zu, zw |
---|
| 2950 | |
---|
| 2951 | USE indices, & |
---|
| 2952 | ONLY: nzt |
---|
| 2953 | |
---|
| 2954 | IMPLICIT NONE |
---|
| 2955 | |
---|
| 2956 | INTEGER(iwp) :: i !< grid cell (x) |
---|
| 2957 | INTEGER(iwp) :: ii !< io-block counter |
---|
| 2958 | INTEGER(iwp) :: il !< io-block counter |
---|
| 2959 | INTEGER(iwp) :: jl !< io-block counter |
---|
| 2960 | INTEGER(iwp) :: kl !< io-block counter |
---|
| 2961 | INTEGER(iwp) :: kdum !< io-block counter |
---|
| 2962 | INTEGER(iwp) :: locdum !< io-block counter |
---|
| 2963 | INTEGER(iwp) :: j !< grid cell (y) |
---|
| 2964 | INTEGER(iwp) :: size_of_mesh !< temporary value for read |
---|
| 2965 | INTEGER(iwp) :: size_of_pols !< temporary value for read |
---|
| 2966 | INTEGER(iwp) :: ioerr !< IOSTAT flag for IO-commands ( 0 = no error ) |
---|
| 2967 | |
---|
| 2968 | REAL(wp) :: zdum !< dummy for measurement height |
---|
| 2969 | REAL(wp) :: avg_agt_height = 1.8_wp |
---|
| 2970 | |
---|
| 2971 | |
---|
| 2972 | ! |
---|
| 2973 | !-- Check the number of agent groups. |
---|
| 2974 | IF ( number_of_agent_groups > max_number_of_agent_groups ) THEN |
---|
| 2975 | WRITE( message_string, * ) 'max_number_of_agent_groups =', & |
---|
| 2976 | max_number_of_agent_groups , & |
---|
| 2977 | '&number_of_agent_groups reset to ', & |
---|
| 2978 | max_number_of_agent_groups |
---|
| 2979 | CALL message( 'mas_init', 'PA0072', 0, 1, 0, 6, 0 ) |
---|
| 2980 | number_of_agent_groups = max_number_of_agent_groups |
---|
| 2981 | ENDIF |
---|
| 2982 | |
---|
| 2983 | ! |
---|
| 2984 | !-- Set some parameters |
---|
| 2985 | d_sigma_rep_agent = 1.0_wp/sigma_rep_agent |
---|
| 2986 | d_sigma_rep_wall = 1.0_wp/sigma_rep_wall |
---|
| 2987 | d_tau_accel_agent = 1.0_wp/tau_accel_agent |
---|
| 2988 | IF ( dt_agent /= 999.0_wp ) THEN |
---|
| 2989 | agent_own_timestep = .TRUE. |
---|
| 2990 | ENDIF |
---|
| 2991 | |
---|
| 2992 | ! |
---|
| 2993 | !-- Get index of first grid box above topography |
---|
| 2994 | ALLOCATE( top_top_s(nysg:nyng,nxlg:nxrg), & |
---|
| 2995 | top_top_w(nysg:nyng,nxlg:nxrg), & |
---|
| 2996 | s_measure_height(nys:nyn,nxl:nxr) ) |
---|
| 2997 | ! |
---|
| 2998 | !-- Get first index above topography for scalar grid and last index in |
---|
| 2999 | !-- topography for z-component of wind |
---|
| 3000 | DO il = nxlg, nxrg |
---|
| 3001 | DO jl = nysg, nyng |
---|
| 3002 | top_top_s(jl,il) = get_topography_top_index_ji(jl,il,'s') + 1 |
---|
| 3003 | top_top_w(jl,il) = get_topography_top_index_ji(jl,il,'w') |
---|
| 3004 | ENDDO |
---|
| 3005 | ENDDO |
---|
| 3006 | ! |
---|
| 3007 | !-- Create 2D array containing the index at which measurements are done by |
---|
| 3008 | !-- agents. The height of this measurement is given by avg_agt_height. |
---|
| 3009 | DO il = nxl, nxr |
---|
| 3010 | DO jl = nys, nyn |
---|
| 3011 | |
---|
| 3012 | kdum = top_top_w(jl,il) |
---|
| 3013 | zdum = zw(kdum) |
---|
| 3014 | zdum = zdum + avg_agt_height |
---|
| 3015 | locdum = 0 |
---|
| 3016 | ! |
---|
| 3017 | !-- Locate minimum distance from u-grid to measurement height (zdum) |
---|
| 3018 | DO kl = 1, nzt |
---|
| 3019 | IF ( ABS(zu(kl)-zdum) < ABS(zu(locdum)-zdum) ) locdum = kl |
---|
| 3020 | ENDDO |
---|
| 3021 | s_measure_height(jl,il) = locdum |
---|
| 3022 | |
---|
| 3023 | ENDDO |
---|
| 3024 | ENDDO |
---|
| 3025 | ! ! |
---|
| 3026 | ! !-- Get indices of PM10 and PM2.5 species, if active |
---|
| 3027 | ! IF ( air_chemistry ) THEN |
---|
| 3028 | ! DO il = 1, nspec |
---|
| 3029 | ! print*,chem_species(il)%name |
---|
| 3030 | ! ! IF ( spec_name(1:4) == 'PM10' ) THEN |
---|
| 3031 | ! ! ind_pm10 = il |
---|
| 3032 | ! ! ELSEIF ( spec_name(1:4) == 'PM25' ) THEN |
---|
| 3033 | ! ! ind_pm25 = il |
---|
| 3034 | ! ! ENDIF |
---|
| 3035 | ! |
---|
| 3036 | ! ENDDO |
---|
| 3037 | ! ENDIF |
---|
| 3038 | ! stop |
---|
| 3039 | CALL mas_create_obstacle_flags |
---|
| 3040 | |
---|
| 3041 | ! |
---|
| 3042 | !-- Set default start positions, if necessary |
---|
| 3043 | IF ( asl(1) == 9999999.9_wp ) asl(1) = 0.0_wp |
---|
| 3044 | IF ( asr(1) == 9999999.9_wp ) asr(1) = ( nx + 1 ) * dx |
---|
| 3045 | IF ( ass(1) == 9999999.9_wp ) ass(1) = 0.0_wp |
---|
| 3046 | IF ( asn(1) == 9999999.9_wp ) asn(1) = ( ny + 1 ) * dy |
---|
| 3047 | IF ( adx(1) == 9999999.9_wp .OR. adx(1) == 0.0_wp ) adx(1) = dx |
---|
| 3048 | IF ( ady(1) == 9999999.9_wp .OR. ady(1) == 0.0_wp ) ady(1) = dy |
---|
| 3049 | |
---|
| 3050 | DO j = 2, number_of_agent_groups |
---|
| 3051 | IF ( asl(j) == 9999999.9_wp ) asl(j) = asl(j-1) |
---|
| 3052 | IF ( asr(j) == 9999999.9_wp ) asr(j) = asr(j-1) |
---|
| 3053 | IF ( ass(j) == 9999999.9_wp ) ass(j) = ass(j-1) |
---|
| 3054 | IF ( asn(j) == 9999999.9_wp ) asn(j) = asn(j-1) |
---|
| 3055 | IF ( adx(j) == 9999999.9_wp .OR. adx(j) == 0.0_wp ) adx(j) = adx(j-1) |
---|
| 3056 | IF ( ady(j) == 9999999.9_wp .OR. ady(j) == 0.0_wp ) ady(j) = ady(j-1) |
---|
| 3057 | ENDDO |
---|
| 3058 | |
---|
| 3059 | ! |
---|
| 3060 | !-- Check boundary condition and set internal variables |
---|
| 3061 | SELECT CASE ( bc_mas_lr ) |
---|
| 3062 | |
---|
| 3063 | CASE ( 'cyclic' ) |
---|
| 3064 | ibc_mas_lr = 0 |
---|
| 3065 | |
---|
| 3066 | CASE ( 'absorb' ) |
---|
| 3067 | ibc_mas_lr = 1 |
---|
| 3068 | |
---|
| 3069 | CASE DEFAULT |
---|
| 3070 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 3071 | 'bc_mas_lr = "', TRIM( bc_mas_lr ), '"' |
---|
| 3072 | CALL message( 'mas_init', 'PA0073', 1, 2, 0, 6, 0 ) |
---|
| 3073 | |
---|
| 3074 | END SELECT |
---|
| 3075 | SELECT CASE ( bc_mas_ns ) |
---|
| 3076 | |
---|
| 3077 | CASE ( 'cyclic' ) |
---|
| 3078 | ibc_mas_ns = 0 |
---|
| 3079 | |
---|
| 3080 | CASE ( 'absorb' ) |
---|
| 3081 | ibc_mas_ns = 1 |
---|
| 3082 | |
---|
| 3083 | CASE DEFAULT |
---|
| 3084 | WRITE( message_string, * ) 'unknown boundary condition ', & |
---|
| 3085 | 'bc_mas_ns = "', TRIM( bc_mas_ns ), '"' |
---|
| 3086 | CALL message( 'mas_init', 'PA0074', 1, 2, 0, 6, 0 ) |
---|
| 3087 | |
---|
| 3088 | END SELECT |
---|
| 3089 | |
---|
| 3090 | ! |
---|
| 3091 | !-- For the first model run of a possible job chain initialize the |
---|
| 3092 | !-- agents, otherwise read the agent data from restart file. |
---|
| 3093 | IF ( TRIM( initializing_actions ) == 'read_restart_data' & |
---|
| 3094 | .AND. read_agents_from_restartfile ) THEN |
---|
| 3095 | |
---|
| 3096 | ! CALL mas_read_restart_file |
---|
| 3097 | |
---|
| 3098 | ELSE |
---|
| 3099 | ! |
---|
| 3100 | !-- Read preprocessed data of navigation mesh and building polygons |
---|
| 3101 | !-- for agent pathfinding |
---|
| 3102 | DO ii = 0, io_blocks-1 |
---|
| 3103 | IF ( ii == io_group ) THEN |
---|
| 3104 | OPEN ( 119, FILE='NAVIGATION_DATA'//TRIM( coupling_char ), & |
---|
| 3105 | FORM='UNFORMATTED', IOSTAT=ioerr ) |
---|
| 3106 | ! |
---|
| 3107 | !-- Read mesh data |
---|
| 3108 | READ(119) size_of_mesh |
---|
| 3109 | ALLOCATE( mesh(1:size_of_mesh)) |
---|
| 3110 | DO i = 1, size_of_mesh |
---|
| 3111 | READ(119) mesh(i)%polygon_id, mesh(i)%vertex_id, & |
---|
| 3112 | mesh(i)%noc, mesh(i)%origin_id, & |
---|
| 3113 | mesh(i)%cost_so_far, mesh(i)%x, & |
---|
| 3114 | mesh(i)%y, mesh(i)%x_s, mesh(i)%y_s |
---|
| 3115 | ALLOCATE( mesh(i)%connected_vertices(1:mesh(i)%noc), & |
---|
| 3116 | mesh(i)%distance_to_vertex(1:mesh(i)%noc) ) |
---|
| 3117 | DO j = 1, mesh(i)%noc |
---|
| 3118 | READ(119) mesh(i)%connected_vertices(j), & |
---|
| 3119 | mesh(i)%distance_to_vertex(j) |
---|
| 3120 | ENDDO |
---|
| 3121 | ENDDO |
---|
| 3122 | ! |
---|
| 3123 | !-- Read polygon data |
---|
| 3124 | READ(119) size_of_pols |
---|
| 3125 | ALLOCATE( polygons(1:size_of_pols) ) |
---|
| 3126 | DO i = 1, size_of_pols |
---|
| 3127 | READ(119) polygons(i)%nov |
---|
| 3128 | ALLOCATE( polygons(i)%vertices(0:polygons(i)%nov+1) ) |
---|
| 3129 | DO j = 0, polygons(i)%nov+1 |
---|
| 3130 | READ(119) polygons(i)%vertices(j)%delete, & |
---|
| 3131 | polygons(i)%vertices(j)%x, & |
---|
| 3132 | polygons(i)%vertices(j)%y |
---|
| 3133 | ENDDO |
---|
| 3134 | ENDDO |
---|
| 3135 | CLOSE(119) |
---|
| 3136 | |
---|
| 3137 | ENDIF |
---|
| 3138 | #if defined( __parallel ) && ! defined ( __check ) |
---|
| 3139 | CALL MPI_BARRIER( comm2d, ierr ) |
---|
| 3140 | #endif |
---|
| 3141 | ENDDO |
---|
| 3142 | |
---|
| 3143 | ! |
---|
| 3144 | !-- Allocate agent arrays and set attributes of the initial set of |
---|
| 3145 | !-- agents, which can be also periodically released at later times. |
---|
[3187] | 3146 | ALLOCATE( agt_count (nysg:nyng,nxlg:nxrg), & |
---|
[3159] | 3147 | grid_agents(nysg:nyng,nxlg:nxrg) ) |
---|
[3187] | 3148 | ! |
---|
| 3149 | !-- Allocate dummy arrays for pathfinding |
---|
| 3150 | ALLOCATE( dummy_path_x(0:agt_path_size), & |
---|
| 3151 | dummy_path_y(0:agt_path_size) ) |
---|
[3159] | 3152 | |
---|
| 3153 | number_of_agents = 0 |
---|
| 3154 | sort_count_mas = 0 |
---|
| 3155 | agt_count = 0 |
---|
| 3156 | |
---|
| 3157 | ! |
---|
| 3158 | !-- initialize counter for agent IDs |
---|
| 3159 | grid_agents%id_counter = 1 |
---|
| 3160 | |
---|
| 3161 | ! |
---|
| 3162 | !-- Initialize all agents with dummy values (otherwise errors may |
---|
| 3163 | !-- occur within restart runs). The reason for this is still not clear |
---|
| 3164 | !-- and may be presumably caused by errors in the respective user-interface. |
---|
| 3165 | zero_agent%agent_mask = .FALSE. |
---|
| 3166 | zero_agent%block_nr = -1 |
---|
| 3167 | zero_agent%group = 0 |
---|
| 3168 | zero_agent%id = 0_idp |
---|
[3187] | 3169 | zero_agent%path_counter = agt_path_size |
---|
[3159] | 3170 | zero_agent%age = 0.0_wp |
---|
| 3171 | zero_agent%age_m = 0.0_wp |
---|
| 3172 | zero_agent%dt_sum = 0.0_wp |
---|
| 3173 | zero_agent%clo = 0.0_wp |
---|
| 3174 | zero_agent%energy_storage= 0.0_wp |
---|
| 3175 | zero_agent%force_x = 0.0_wp |
---|
| 3176 | zero_agent%force_y = 0.0_wp |
---|
| 3177 | zero_agent%origin_x = 0.0_wp |
---|
| 3178 | zero_agent%origin_y = 0.0_wp |
---|
| 3179 | zero_agent%speed_abs = 0.0_wp |
---|
| 3180 | zero_agent%speed_e_x = 0.0_wp |
---|
| 3181 | zero_agent%speed_e_y = 0.0_wp |
---|
| 3182 | zero_agent%speed_des = random_normal(desired_speed, des_sp_sig) |
---|
| 3183 | zero_agent%speed_x = 0.0_wp |
---|
| 3184 | zero_agent%speed_y = 0.0_wp |
---|
| 3185 | zero_agent%thermal_index = 0.0_wp |
---|
| 3186 | zero_agent%x = 0.0_wp |
---|
| 3187 | zero_agent%y = 0.0_wp |
---|
| 3188 | zero_agent%path_x = 0.0_wp |
---|
| 3189 | zero_agent%path_y = 0.0_wp |
---|
| 3190 | zero_agent%t_x = 0.0_wp |
---|
| 3191 | zero_agent%t_y = 0.0_wp |
---|
| 3192 | |
---|
| 3193 | ! |
---|
| 3194 | !-- Set a seed value for the random number generator to be exclusively |
---|
| 3195 | !-- used for the agent code. The generated random numbers should be |
---|
| 3196 | !-- different on the different PEs. |
---|
| 3197 | iran_agent = iran_agent + myid |
---|
| 3198 | |
---|
| 3199 | CALL mas_create_agent (PHASE_INIT) |
---|
| 3200 | |
---|
| 3201 | ENDIF |
---|
| 3202 | |
---|
| 3203 | ! |
---|
| 3204 | !-- To avoid programm abort, assign agents array to the local version of |
---|
| 3205 | !-- first grid cell |
---|
| 3206 | number_of_agents = agt_count(nys,nxl) |
---|
| 3207 | agents => grid_agents(nys,nxl)%agents(1:number_of_agents) |
---|
| 3208 | |
---|
| 3209 | END SUBROUTINE mas_init |
---|
| 3210 | |
---|
[3187] | 3211 | !------------------------------------------------------------------------------! |
---|
| 3212 | ! Description: |
---|
| 3213 | ! ------------ |
---|
[3235] | 3214 | !> Output of informative message about maximum agent number |
---|
| 3215 | !------------------------------------------------------------------------------! |
---|
| 3216 | SUBROUTINE mas_last_actions |
---|
| 3217 | |
---|
| 3218 | USE control_parameters, & |
---|
| 3219 | ONLY: message_string |
---|
| 3220 | |
---|
| 3221 | IMPLICIT NONE |
---|
| 3222 | |
---|
| 3223 | WRITE(message_string,'(A,I8,A)') & |
---|
| 3224 | 'The maximumn number of agents durin this run was', & |
---|
| 3225 | maximum_number_of_agents, & |
---|
| 3226 | '&Consider adjusting the INPUT parameter'// & |
---|
| 3227 | '&dim_size_agtnum_manual accordingly for the next run.' |
---|
| 3228 | |
---|
| 3229 | CALL message( 'mas_data_output_agents', 'PA0457', 0, 0, 0, 6, 0 ) |
---|
| 3230 | |
---|
| 3231 | END SUBROUTINE mas_last_actions |
---|
| 3232 | |
---|
| 3233 | !------------------------------------------------------------------------------! |
---|
| 3234 | ! Description: |
---|
| 3235 | ! ------------ |
---|
[3187] | 3236 | !> Finds the shortest path from a start position to a target position using the |
---|
| 3237 | !> A*-algorithm |
---|
| 3238 | !------------------------------------------------------------------------------! |
---|
| 3239 | SUBROUTINE mas_nav_a_star( start_x, start_y, target_x, target_y, nsteps ) |
---|
| 3240 | |
---|
| 3241 | IMPLICIT NONE |
---|
| 3242 | |
---|
| 3243 | LOGICAL :: target_reached !< flag |
---|
| 3244 | |
---|
| 3245 | INTEGER(iwp) :: cur_node !< current node of binary heap |
---|
| 3246 | INTEGER(iwp) :: il !< counter (x) |
---|
| 3247 | INTEGER(iwp) :: neigh_node !< neighbor node |
---|
| 3248 | INTEGER(iwp) :: node_counter !< binary heap node counter |
---|
| 3249 | INTEGER(iwp) :: path_ag !< index of agent path |
---|
| 3250 | INTEGER(iwp) :: som !< size of mesh |
---|
| 3251 | INTEGER(iwp) :: steps !< steps along the path |
---|
| 3252 | INTEGER(iwp) :: nsteps !< number of steps |
---|
| 3253 | |
---|
| 3254 | REAL(wp) :: start_x !< x-coordinate agent |
---|
| 3255 | REAL(wp) :: start_y !< y-coordinate agent |
---|
| 3256 | REAL(wp) :: new_cost !< updated cost to reach node |
---|
| 3257 | REAL(wp) :: new_priority !< priority of node to be added to queue |
---|
| 3258 | REAL(wp) :: rn_gate !< random number for corner gate |
---|
| 3259 | REAL(wp) :: target_x !< x-coordinate target |
---|
| 3260 | REAL(wp) :: target_y !< y-coordinate target |
---|
| 3261 | ! |
---|
| 3262 | !-- Coordinate Type |
---|
| 3263 | TYPE coord |
---|
| 3264 | REAL(wp) :: x !< x-coordinate |
---|
| 3265 | REAL(wp) :: x_s !< x-coordinate (shifted) |
---|
| 3266 | REAL(wp) :: y !< y-coordinate |
---|
| 3267 | REAL(wp) :: y_s !< y-coordinate (shifted) |
---|
| 3268 | END TYPE coord |
---|
| 3269 | |
---|
| 3270 | TYPE(coord), DIMENSION(:), ALLOCATABLE, TARGET :: path !< path array |
---|
| 3271 | TYPE(coord), DIMENSION(:), ALLOCATABLE, TARGET :: tmp_path !< temporary path for resizing |
---|
| 3272 | |
---|
| 3273 | node_counter = 0 |
---|
| 3274 | ! |
---|
| 3275 | !-- Create temporary navigation mesh including agent and target positions |
---|
| 3276 | CALL mas_nav_create_tmp_mesh( start_x, start_y, target_x, target_y, som ) |
---|
| 3277 | tmp_mesh(som)%cost_so_far = 0.0_wp |
---|
| 3278 | ! |
---|
| 3279 | !-- Initialize priority queue |
---|
| 3280 | heap_count = 0_iwp |
---|
| 3281 | ALLOCATE(queue(0:100)) |
---|
| 3282 | target_reached = .FALSE. |
---|
| 3283 | ! |
---|
| 3284 | !-- Add starting point (agent position) to frontier (the frontier consists |
---|
| 3285 | !-- of all the nodes that are to be visited. The node with the smallest |
---|
| 3286 | !-- priority will be visited first. The priority consists of the distance |
---|
| 3287 | !-- from the start node to this node plus a minimal guess (direct distance) |
---|
| 3288 | !-- from this node to the goal). For the starting node, the priority is set |
---|
| 3289 | !-- to 0, as it's the only node thus far |
---|
| 3290 | CALL mas_heap_insert_item(som,0.0_wp) |
---|
| 3291 | cur_node = som |
---|
| 3292 | DO WHILE ( heap_count > 0 ) |
---|
| 3293 | ! |
---|
| 3294 | !-- Step one: Pick lowest priority item from queue |
---|
| 3295 | node_counter = node_counter + 1 |
---|
| 3296 | CALL mas_heap_extract_item(cur_node) |
---|
| 3297 | ! |
---|
| 3298 | !-- Node 0 is the goal node |
---|
| 3299 | IF ( cur_node == 0 ) THEN |
---|
| 3300 | EXIT |
---|
| 3301 | ENDIF |
---|
| 3302 | ! |
---|
| 3303 | !-- Loop over all of cur_node's neighbors |
---|
| 3304 | DO il = 1, tmp_mesh(cur_node)%noc |
---|
| 3305 | neigh_node = tmp_mesh(cur_node)%connected_vertices(il) |
---|
| 3306 | ! |
---|
| 3307 | !-- Check, if the way from the start node to this neigh_node via |
---|
| 3308 | !-- cur_node is shorter than the previously found shortest path to it. |
---|
| 3309 | !-- If so, replace said cost and add neigh_node to the frontier. |
---|
| 3310 | !-- cost_so_far is initialized as 1.d12 so that all found distances |
---|
| 3311 | !-- should be smaller. |
---|
| 3312 | new_cost = tmp_mesh(cur_node)%cost_so_far & |
---|
| 3313 | + tmp_mesh(cur_node)%distance_to_vertex(il) |
---|
| 3314 | IF ( new_cost < tmp_mesh(neigh_node)%cost_so_far ) THEN |
---|
| 3315 | tmp_mesh(neigh_node)%cost_so_far = new_cost |
---|
| 3316 | tmp_mesh(neigh_node)%origin_id = cur_node |
---|
| 3317 | ! |
---|
| 3318 | !-- Priority in the queue is cost_so_far + heuristic to goal |
---|
| 3319 | new_priority = new_cost & |
---|
| 3320 | + heuristic(tmp_mesh(neigh_node)%x, & |
---|
| 3321 | tmp_mesh(neigh_node)%y, tmp_mesh(0)%x, & |
---|
| 3322 | tmp_mesh(0)%y) |
---|
| 3323 | CALL mas_heap_insert_item(neigh_node,new_priority) |
---|
| 3324 | ENDIF |
---|
| 3325 | ENDDO |
---|
| 3326 | ENDDO |
---|
| 3327 | ! |
---|
| 3328 | !-- Add nodes to a path array. To do this, we must backtrack from the target |
---|
| 3329 | !-- node to its origin to its origin and so on until an node is reached that |
---|
| 3330 | !-- has no origin (%origin_id == -1). This is the starting node. |
---|
| 3331 | DEALLOCATE(queue) |
---|
| 3332 | cur_node = 0 |
---|
| 3333 | steps = 0 |
---|
| 3334 | ALLOCATE(path(1:100)) |
---|
| 3335 | DO WHILE ( cur_node /= -1 ) |
---|
| 3336 | steps = steps + 1 |
---|
| 3337 | ! |
---|
| 3338 | !-- Resize path array if necessary |
---|
| 3339 | IF ( steps > SIZE(path) ) THEN |
---|
| 3340 | ALLOCATE(tmp_path(1:steps-1)) |
---|
| 3341 | tmp_path(1:steps-1) = path(1:steps-1) |
---|
| 3342 | DEALLOCATE(path) |
---|
| 3343 | ALLOCATE(path(1:2*(steps-1))) |
---|
| 3344 | path(1:steps-1) = tmp_path(1:steps-1) |
---|
| 3345 | DEALLOCATE(tmp_path) |
---|
| 3346 | ENDIF |
---|
| 3347 | path(steps)%x = tmp_mesh(cur_node)%x |
---|
| 3348 | path(steps)%y = tmp_mesh(cur_node)%y |
---|
| 3349 | path(steps)%x_s = tmp_mesh(cur_node)%x_s |
---|
| 3350 | path(steps)%y_s = tmp_mesh(cur_node)%y_s |
---|
| 3351 | cur_node = tmp_mesh(cur_node)%origin_id |
---|
| 3352 | ENDDO |
---|
| 3353 | ! |
---|
| 3354 | !-- Add calculated intermittent targets to the path until either the |
---|
| 3355 | !-- target or the maximum number of intermittent targets is reached. |
---|
| 3356 | !-- Ignore starting point (reduce index by one), it is agent position. |
---|
| 3357 | dummy_path_x = -1 |
---|
| 3358 | dummy_path_y = -1 |
---|
| 3359 | path_ag = 1 |
---|
| 3360 | steps = steps - 1 |
---|
| 3361 | nsteps = 0 |
---|
| 3362 | DO WHILE( steps > 0 .AND. path_ag <= agt_path_size ) |
---|
| 3363 | ! |
---|
| 3364 | !-- Each target point is randomly chosen along a line target along the |
---|
| 3365 | !-- bisector of the building corner that starts at corner_gate_start |
---|
| 3366 | !-- and has a width of corner_gate_width. This is to avoid clustering |
---|
| 3367 | !-- when opposing agent groups try to reach the same corner target. |
---|
| 3368 | rn_gate = random_function(iran_agent) * corner_gate_width & |
---|
| 3369 | + corner_gate_start |
---|
| 3370 | dummy_path_x(path_ag) = path(steps)%x + rn_gate & |
---|
| 3371 | * (path(steps)%x_s - path(steps)%x) |
---|
| 3372 | dummy_path_y(path_ag) = path(steps)%y + rn_gate & |
---|
| 3373 | * (path(steps)%y_s - path(steps)%y) |
---|
| 3374 | steps = steps - 1 |
---|
| 3375 | path_ag = path_ag + 1 |
---|
| 3376 | nsteps = nsteps + 1 |
---|
| 3377 | ENDDO |
---|
| 3378 | ! |
---|
| 3379 | !-- Set current intermittent target of this agent |
---|
| 3380 | DEALLOCATE(tmp_mesh, path) |
---|
| 3381 | |
---|
| 3382 | END SUBROUTINE mas_nav_a_star |
---|
| 3383 | |
---|
[3159] | 3384 | !------------------------------------------------------------------------------! |
---|
| 3385 | ! Description: |
---|
| 3386 | ! ------------ |
---|
| 3387 | !> Adds a connection between two points of the navigation mesh |
---|
| 3388 | !> (one-way: in_mp1 to in_mp2) |
---|
| 3389 | !------------------------------------------------------------------------------! |
---|
| 3390 | SUBROUTINE mas_nav_add_connection ( in_mp1, id2, in_mp2 ) |
---|
| 3391 | |
---|
| 3392 | IMPLICIT NONE |
---|
| 3393 | |
---|
| 3394 | LOGICAL :: connection_established !< Flag to indicate if connection has already been established |
---|
| 3395 | |
---|
| 3396 | INTEGER(iwp) :: id2 !< ID of in_mp2 |
---|
| 3397 | INTEGER(iwp) :: il !< local counter |
---|
| 3398 | INTEGER(iwp) :: noc1 !< number of connections in in_mp1 |
---|
| 3399 | |
---|
| 3400 | INTEGER, DIMENSION(:), ALLOCATABLE :: dum_cv !< dummy array for connected_vertices |
---|
| 3401 | |
---|
| 3402 | REAL(wp) :: dist !< Distance between the two points |
---|
| 3403 | |
---|
| 3404 | REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_dtv |
---|
| 3405 | |
---|
| 3406 | TYPE(mesh_point) :: in_mp1 !< mesh point that gets a new connection |
---|
| 3407 | TYPE(mesh_point) :: in_mp2 !< mesh point in_mp1 will be connected to |
---|
| 3408 | |
---|
| 3409 | connection_established = .FALSE. |
---|
| 3410 | ! |
---|
| 3411 | !-- Check if connection has already been established |
---|
| 3412 | noc1 = SIZE(in_mp1%connected_vertices) |
---|
| 3413 | DO il = 1, in_mp1%noc |
---|
| 3414 | IF ( in_mp1%connected_vertices(il) == id2 ) THEN |
---|
| 3415 | connection_established = .TRUE. |
---|
| 3416 | EXIT |
---|
| 3417 | ENDIF |
---|
| 3418 | ENDDO |
---|
| 3419 | |
---|
| 3420 | IF ( .NOT. connection_established ) THEN |
---|
| 3421 | ! |
---|
| 3422 | !-- Resize arrays, if necessary |
---|
| 3423 | IF ( in_mp1%noc >= noc1 ) THEN |
---|
| 3424 | ALLOCATE( dum_cv(1:noc1),dum_dtv(1:noc1) ) |
---|
| 3425 | dum_cv = in_mp1%connected_vertices |
---|
| 3426 | dum_dtv = in_mp1%distance_to_vertex |
---|
| 3427 | DEALLOCATE( in_mp1%connected_vertices, in_mp1%distance_to_vertex ) |
---|
| 3428 | ALLOCATE( in_mp1%connected_vertices(1:2*noc1), & |
---|
| 3429 | in_mp1%distance_to_vertex(1:2*noc1) ) |
---|
| 3430 | in_mp1%connected_vertices = -999 |
---|
| 3431 | in_mp1%distance_to_vertex = -999. |
---|
| 3432 | in_mp1%connected_vertices(1:noc1) = dum_cv |
---|
| 3433 | in_mp1%distance_to_vertex(1:noc1) = dum_dtv |
---|
| 3434 | ENDIF |
---|
| 3435 | |
---|
| 3436 | ! |
---|
| 3437 | !-- Add connection |
---|
| 3438 | in_mp1%noc = in_mp1%noc+1 |
---|
| 3439 | dist = SQRT( (in_mp1%x - in_mp2%x)**2 + (in_mp1%y - in_mp2%y)**2 ) |
---|
| 3440 | in_mp1%connected_vertices(in_mp1%noc) = id2 |
---|
| 3441 | in_mp1%distance_to_vertex(in_mp1%noc) = dist |
---|
| 3442 | ENDIF |
---|
| 3443 | |
---|
| 3444 | END SUBROUTINE mas_nav_add_connection |
---|
| 3445 | |
---|
| 3446 | !------------------------------------------------------------------------------! |
---|
| 3447 | ! Description: |
---|
| 3448 | ! ------------ |
---|
| 3449 | !> Adds a vertex (curren position of agent or target) to the existing tmp_mesh |
---|
| 3450 | !------------------------------------------------------------------------------! |
---|
| 3451 | SUBROUTINE mas_nav_add_vertex_to_mesh ( in_mp, in_id ) |
---|
| 3452 | |
---|
| 3453 | IMPLICIT NONE |
---|
| 3454 | |
---|
| 3455 | LOGICAL :: intersection_found !< flag |
---|
| 3456 | |
---|
| 3457 | INTEGER(iwp) :: jl !< mesh point counter |
---|
| 3458 | INTEGER(iwp) :: pl !< polygon counter |
---|
| 3459 | INTEGER(iwp) :: vl !< vertex counter |
---|
| 3460 | INTEGER(iwp) :: pid_t !< polygon id of tested mesh point |
---|
| 3461 | INTEGER(iwp) :: vid_t !< vertex id of tested mesh point |
---|
| 3462 | INTEGER(iwp) :: in_id !< vertex id of tested mesh point |
---|
| 3463 | |
---|
| 3464 | REAL(wp) :: v1x !< x-coordinate of test vertex 1 for intersection test |
---|
| 3465 | REAL(wp) :: v1y !< y-coordinate of test vertex 1 for intersection test |
---|
| 3466 | REAL(wp) :: v2x !< x-coordinate of test vertex 2 for intersection test |
---|
| 3467 | REAL(wp) :: v2y !< y-coordinate of test vertex 2 for intersection test |
---|
| 3468 | REAL(wp) :: x !< x-coordinate of current mesh point |
---|
| 3469 | REAL(wp) :: x_t !< x-coordinate of tested mesh point |
---|
| 3470 | REAL(wp) :: y !< y-coordinate of current mesh point |
---|
| 3471 | REAL(wp) :: y_t !< y-coordinate of tested mesh point |
---|
| 3472 | |
---|
| 3473 | TYPE(mesh_point) :: in_mp !< Input mesh point |
---|
| 3474 | ! |
---|
| 3475 | !-- |
---|
| 3476 | x = in_mp%x |
---|
| 3477 | y = in_mp%y |
---|
| 3478 | DO jl = 0, SIZE(tmp_mesh)-2 |
---|
| 3479 | IF ( in_id == jl ) CYCLE |
---|
| 3480 | ! |
---|
| 3481 | !-- Ignore mesh points with 0 connections |
---|
| 3482 | IF ( tmp_mesh(jl)%polygon_id /= -1 ) THEN |
---|
| 3483 | IF ( tmp_mesh(jl)%noc == 0 ) CYCLE |
---|
| 3484 | ENDIF |
---|
| 3485 | x_t = tmp_mesh(jl)%x |
---|
| 3486 | y_t = tmp_mesh(jl)%y |
---|
| 3487 | pid_t = tmp_mesh(jl)%polygon_id |
---|
| 3488 | vid_t = tmp_mesh(jl)%vertex_id |
---|
| 3489 | ! |
---|
| 3490 | !-- If the connecting line between the target and a mesh point points |
---|
| 3491 | !-- into the mesh point's polygon, no connection will be |
---|
| 3492 | !-- established between the two points. This is the case if the |
---|
| 3493 | !-- previous (next) vertex of the polygon is right of the connecting |
---|
| 3494 | !-- line and the next (previous) vertex of the polygon is left of the |
---|
| 3495 | !-- connecting line. |
---|
| 3496 | IF ( pid_t > 0 .AND. pid_t <= SIZE(polygons) ) THEN |
---|
| 3497 | IF ( (((is_left(x,y,x_t,y_t,polygons(pid_t)%vertices(vid_t-1)%x, & |
---|
| 3498 | polygons(pid_t)%vertices(vid_t-1)%y) & |
---|
| 3499 | .AND. is_right(x,y,x_t,y_t, & |
---|
| 3500 | polygons(pid_t)%vertices(vid_t+1)%x, & |
---|
| 3501 | polygons(pid_t)%vertices(vid_t+1)%y) ) & |
---|
| 3502 | .OR. (is_right(x,y,x_t,y_t, & |
---|
| 3503 | polygons(pid_t)%vertices(vid_t-1)%x, & |
---|
| 3504 | polygons(pid_t)%vertices(vid_t-1)%y) & |
---|
| 3505 | .AND. is_left(x,y,x_t,y_t, & |
---|
| 3506 | polygons(pid_t)%vertices(vid_t+1)%x, & |
---|
| 3507 | polygons(pid_t)%vertices(vid_t+1)%y)))))& |
---|
| 3508 | THEN |
---|
| 3509 | CYCLE |
---|
| 3510 | ENDIF |
---|
| 3511 | ENDIF |
---|
| 3512 | ! |
---|
| 3513 | !-- For each edge of each polygon, check if it intersects with the |
---|
| 3514 | !-- potential connection. If at least one intersection is found, |
---|
| 3515 | !-- no connection can be made |
---|
| 3516 | intersection_found = .FALSE. |
---|
| 3517 | DO pl = 1, SIZE(polygons) |
---|
| 3518 | DO vl = 1, polygons(pl)%nov |
---|
| 3519 | v1x = polygons(pl)%vertices(vl)%x |
---|
| 3520 | v1y = polygons(pl)%vertices(vl)%y |
---|
| 3521 | v2x = polygons(pl)%vertices(vl+1)%x |
---|
| 3522 | v2y = polygons(pl)%vertices(vl+1)%y |
---|
| 3523 | intersection_found = intersect(x,y,x_t,y_t,v1x,v1y,v2x,v2y) |
---|
| 3524 | IF ( intersection_found ) THEN |
---|
| 3525 | EXIT |
---|
| 3526 | ENDIF |
---|
| 3527 | ENDDO |
---|
| 3528 | IF ( intersection_found ) EXIT |
---|
| 3529 | ENDDO |
---|
| 3530 | IF ( intersection_found ) CYCLE |
---|
| 3531 | ! |
---|
| 3532 | !-- If neither of the above two test was true, a connection will be |
---|
| 3533 | !-- established between the two mesh points. |
---|
| 3534 | CALL mas_nav_add_connection(in_mp,jl, tmp_mesh(jl)) |
---|
| 3535 | CALL mas_nav_add_connection(tmp_mesh(jl),in_id, in_mp) |
---|
| 3536 | ENDDO |
---|
| 3537 | CALL mas_nav_reduce_connections(in_mp) |
---|
| 3538 | |
---|
| 3539 | END SUBROUTINE mas_nav_add_vertex_to_mesh |
---|
| 3540 | |
---|
| 3541 | !------------------------------------------------------------------------------! |
---|
| 3542 | ! Description: |
---|
| 3543 | ! ------------ |
---|
| 3544 | !> Creates a temporary copy of the navigation mesh to be used for pathfinding |
---|
| 3545 | !------------------------------------------------------------------------------! |
---|
| 3546 | SUBROUTINE mas_nav_create_tmp_mesh( a_x, a_y, t_x, t_y, som ) |
---|
| 3547 | |
---|
| 3548 | IMPLICIT NONE |
---|
| 3549 | |
---|
| 3550 | INTEGER(iwp) :: som !< size of mesh |
---|
| 3551 | INTEGER(iwp) :: noc !< number of connetions |
---|
| 3552 | INTEGER(iwp) :: im !< local mesh point counter |
---|
| 3553 | |
---|
| 3554 | REAL(wp) :: a_x !< x-coordinate agent |
---|
| 3555 | REAL(wp) :: a_y !< y-coordinate agent |
---|
| 3556 | REAL(wp) :: t_x !< x-coordinate target |
---|
| 3557 | REAL(wp) :: t_y !< y-coordinate target |
---|
| 3558 | ! |
---|
| 3559 | !-- give tmp_mesh the size of mesh |
---|
| 3560 | som = SIZE(mesh)+1 |
---|
| 3561 | ALLOCATE(tmp_mesh(0:som)) |
---|
| 3562 | ! |
---|
| 3563 | !-- give the allocatable variables in tmp_mesh their respctive sizes |
---|
| 3564 | DO im = 1, som-1 |
---|
| 3565 | noc = mesh(im)%noc |
---|
| 3566 | ALLOCATE(tmp_mesh(im)%connected_vertices(1:noc)) |
---|
| 3567 | ALLOCATE(tmp_mesh(im)%distance_to_vertex(1:noc)) |
---|
| 3568 | ENDDO |
---|
| 3569 | ! |
---|
| 3570 | !-- copy mesh to tmp_mesh |
---|
| 3571 | tmp_mesh(1:som-1) = mesh(1:som-1) |
---|
| 3572 | ! |
---|
| 3573 | !-- Add target point ... |
---|
| 3574 | CALL mas_nav_init_mesh_point(tmp_mesh(0),-1_iwp,-1_iwp,t_x, t_y) |
---|
| 3575 | CALL mas_nav_add_vertex_to_mesh(tmp_mesh(0),0_iwp) |
---|
| 3576 | ! |
---|
| 3577 | !-- ... and start point to temp mesh |
---|
| 3578 | CALL mas_nav_init_mesh_point(tmp_mesh(som),-1_iwp,-1_iwp,a_x, a_y) |
---|
| 3579 | CALL mas_nav_add_vertex_to_mesh(tmp_mesh(som),som) |
---|
| 3580 | |
---|
| 3581 | END SUBROUTINE mas_nav_create_tmp_mesh |
---|
[3187] | 3582 | |
---|
[3159] | 3583 | |
---|
| 3584 | !------------------------------------------------------------------------------! |
---|
| 3585 | ! Description: |
---|
| 3586 | ! ------------ |
---|
[3187] | 3587 | !> Finds the shortest path from an agents' position to her target. As the |
---|
| 3588 | !> actual pathfinding algorithm uses the obstacle corners and then shifts them |
---|
| 3589 | !> outward after pathfinding, cases can uccur in which the connection between |
---|
| 3590 | !> these intermittent targets then intersect with obstacles. To remedy this |
---|
| 3591 | !> the pathfinding algorithm is then run on every two subsequent intermittent |
---|
| 3592 | !> targets iteratively and new intermittent targets may be added to the path |
---|
| 3593 | !> this way. |
---|
[3159] | 3594 | !------------------------------------------------------------------------------! |
---|
[3187] | 3595 | SUBROUTINE mas_nav_find_path( nl ) |
---|
[3159] | 3596 | |
---|
| 3597 | IMPLICIT NONE |
---|
| 3598 | |
---|
[3187] | 3599 | INTEGER(iwp) :: nl !< local agent counter |
---|
| 3600 | INTEGER(iwp) :: il !< local counter |
---|
| 3601 | INTEGER(iwp) :: jl !< local counter |
---|
| 3602 | INTEGER(iwp) :: kl !< local counter |
---|
| 3603 | INTEGER(iwp) :: nsteps_total !< number of steps on path |
---|
| 3604 | INTEGER(iwp) :: nsteps_dummy !< number of steps on path |
---|
| 3605 | |
---|
| 3606 | REAL(wp), DIMENSION(0:30) :: ld_path_x !< local dummy agent path to target (x) |
---|
| 3607 | REAL(wp), DIMENSION(0:30) :: ld_path_y !< local dummy agent path to target (y) |
---|
[3159] | 3608 | ! |
---|
[3187] | 3609 | !-- Initialize agent path arrays |
---|
| 3610 | agents(nl)%path_x = -1 |
---|
| 3611 | agents(nl)%path_y = -1 |
---|
| 3612 | agents(nl)%path_x(0) = agents(nl)%x |
---|
| 3613 | agents(nl)%path_y(0) = agents(nl)%y |
---|
[3159] | 3614 | ! |
---|
[3187] | 3615 | !-- Calculate initial path |
---|
| 3616 | CALL mas_nav_a_star( agents(nl)%x, agents(nl)%y, & |
---|
| 3617 | agents(nl)%t_x, agents(nl)%t_y, nsteps_total ) |
---|
[3159] | 3618 | ! |
---|
[3187] | 3619 | !-- Set the rest of the agent path that was just calculated |
---|
| 3620 | agents(nl)%path_x(1:nsteps_total) = dummy_path_x(1:nsteps_total) |
---|
| 3621 | agents(nl)%path_y(1:nsteps_total) = dummy_path_y(1:nsteps_total) |
---|
[3159] | 3622 | ! |
---|
[3187] | 3623 | !-- Iterate through found path and check more intermittent targets need |
---|
| 3624 | !-- to be added. For this, run pathfinding between every two consecutive |
---|
| 3625 | !-- intermittent targets. |
---|
| 3626 | DO il = 0, MIN(agt_path_size-1, nsteps_total-1) |
---|
[3159] | 3627 | ! |
---|
[3187] | 3628 | !-- pathfinding between two consecutive intermittent targets |
---|
| 3629 | CALL mas_nav_a_star( agents(nl)%path_x(il), agents(nl)%path_y(il), & |
---|
| 3630 | agents(nl)%path_x(il+1), agents(nl)%path_y(il+1),& |
---|
| 3631 | nsteps_dummy ) |
---|
| 3632 | nsteps_dummy = nsteps_dummy - 1 |
---|
[3159] | 3633 | ! |
---|
[3187] | 3634 | !-- If additional intermittent targets are found, add them to the path |
---|
| 3635 | IF ( nsteps_dummy > 0 ) THEN |
---|
| 3636 | ld_path_x = -1 |
---|
| 3637 | ld_path_y = -1 |
---|
| 3638 | ld_path_x(il+1:il+nsteps_dummy) = dummy_path_x(1:nsteps_dummy) |
---|
| 3639 | ld_path_y(il+1:il+nsteps_dummy) = dummy_path_y(1:nsteps_dummy) |
---|
| 3640 | kl = 1 |
---|
| 3641 | DO jl = il+1,nsteps_total |
---|
| 3642 | ld_path_x( il+nsteps_dummy+kl ) = agents(nl)%path_x(jl) |
---|
| 3643 | ld_path_y( il+nsteps_dummy+kl ) = agents(nl)%path_y(jl) |
---|
| 3644 | kl = kl + 1 |
---|
| 3645 | IF ( kl > agt_path_size ) EXIT |
---|
| 3646 | ENDDO |
---|
| 3647 | nsteps_total = MIN(nsteps_total + nsteps_dummy, agt_path_size) |
---|
| 3648 | agents(nl)%path_x(il+1:nsteps_total) = ld_path_x(il+1:nsteps_total) |
---|
| 3649 | agents(nl)%path_y(il+1:nsteps_total) = ld_path_y(il+1:nsteps_total) |
---|
[3159] | 3650 | ENDIF |
---|
[3187] | 3651 | |
---|
[3159] | 3652 | ENDDO |
---|
| 3653 | ! |
---|
[3187] | 3654 | !-- reset path counter to first intermittent target |
---|
| 3655 | agents(nl)%path_counter = 1 |
---|
[3159] | 3656 | |
---|
| 3657 | END SUBROUTINE mas_nav_find_path |
---|
| 3658 | |
---|
| 3659 | !------------------------------------------------------------------------------! |
---|
| 3660 | ! Description: |
---|
| 3661 | ! ------------ |
---|
| 3662 | !> Reduces the size of connection array to the amount of actual connections |
---|
| 3663 | !> after all connetions were added to a mesh point |
---|
| 3664 | !------------------------------------------------------------------------------! |
---|
| 3665 | SUBROUTINE mas_nav_reduce_connections ( in_mp ) |
---|
| 3666 | |
---|
| 3667 | IMPLICIT NONE |
---|
| 3668 | |
---|
| 3669 | INTEGER(iwp) :: noc !< number of connections |
---|
| 3670 | |
---|
| 3671 | INTEGER, DIMENSION(:), ALLOCATABLE :: dum_cv !< dummy connected_vertices |
---|
| 3672 | |
---|
| 3673 | REAL(wp), DIMENSION(:), ALLOCATABLE :: dum_dtv !< dummy distance_to_vertex |
---|
| 3674 | |
---|
| 3675 | TYPE(mesh_point) :: in_mp |
---|
| 3676 | |
---|
| 3677 | noc = in_mp%noc |
---|
| 3678 | ALLOCATE( dum_cv(1:noc),dum_dtv(1:noc) ) |
---|
| 3679 | dum_cv = in_mp%connected_vertices(1:noc) |
---|
| 3680 | dum_dtv = in_mp%distance_to_vertex(1:noc) |
---|
| 3681 | DEALLOCATE( in_mp%connected_vertices, in_mp%distance_to_vertex ) |
---|
| 3682 | ALLOCATE( in_mp%connected_vertices(1:noc), & |
---|
| 3683 | in_mp%distance_to_vertex(1:noc) ) |
---|
| 3684 | in_mp%connected_vertices(1:noc) = dum_cv(1:noc) |
---|
| 3685 | in_mp%distance_to_vertex(1:noc) = dum_dtv(1:noc) |
---|
| 3686 | |
---|
| 3687 | END SUBROUTINE mas_nav_reduce_connections |
---|
| 3688 | |
---|
| 3689 | !------------------------------------------------------------------------------! |
---|
| 3690 | ! Description: |
---|
| 3691 | ! ------------ |
---|
| 3692 | !> Initializes a point of the navigation mesh |
---|
| 3693 | !------------------------------------------------------------------------------! |
---|
| 3694 | SUBROUTINE mas_nav_init_mesh_point ( in_mp, pid, vid, x, y ) |
---|
| 3695 | |
---|
| 3696 | IMPLICIT NONE |
---|
| 3697 | |
---|
| 3698 | INTEGER(iwp) :: pid !< polygon ID |
---|
| 3699 | INTEGER(iwp) :: vid !< vertex ID |
---|
| 3700 | |
---|
| 3701 | REAL(wp) :: x !< x-coordinate |
---|
| 3702 | REAL(wp) :: y !< y-coordinate |
---|
| 3703 | |
---|
| 3704 | TYPE(mesh_point) :: in_mp !< mesh point to be initialized |
---|
| 3705 | |
---|
| 3706 | in_mp%origin_id = -1 |
---|
| 3707 | in_mp%polygon_id = pid |
---|
| 3708 | in_mp%vertex_id = vid |
---|
| 3709 | in_mp%cost_so_far = 1.d12 |
---|
| 3710 | in_mp%x = x |
---|
| 3711 | in_mp%y = y |
---|
| 3712 | in_mp%x_s = x |
---|
| 3713 | in_mp%y_s = y |
---|
| 3714 | ALLOCATE(in_mp%connected_vertices(1:100), & |
---|
| 3715 | in_mp%distance_to_vertex(1:100)) |
---|
| 3716 | in_mp%connected_vertices = -999 |
---|
| 3717 | in_mp%distance_to_vertex = -999. |
---|
| 3718 | in_mp%noc = 0 |
---|
| 3719 | |
---|
| 3720 | END SUBROUTINE mas_nav_init_mesh_point |
---|
| 3721 | |
---|
| 3722 | !------------------------------------------------------------------------------! |
---|
| 3723 | ! Description: |
---|
| 3724 | ! ------------ |
---|
| 3725 | !> Reading of namlist from parin file |
---|
| 3726 | !------------------------------------------------------------------------------! |
---|
| 3727 | SUBROUTINE mas_parin |
---|
| 3728 | |
---|
| 3729 | USE control_parameters, & |
---|
[3198] | 3730 | ONLY: agent_time_unlimited, multi_agent_system_end, & |
---|
| 3731 | multi_agent_system_start |
---|
[3159] | 3732 | |
---|
| 3733 | IMPLICIT NONE |
---|
| 3734 | |
---|
| 3735 | CHARACTER (LEN=80) :: line !< |
---|
| 3736 | |
---|
[3198] | 3737 | NAMELIST /agent_parameters/ a_rand_target, & |
---|
| 3738 | adx, & |
---|
| 3739 | ady, & |
---|
| 3740 | agent_maximum_age, & |
---|
| 3741 | agent_time_unlimited, & |
---|
| 3742 | alloc_factor_mas, & |
---|
| 3743 | asl, & |
---|
| 3744 | asn, & |
---|
| 3745 | asr, & |
---|
| 3746 | ass, & |
---|
| 3747 | at_x, & |
---|
| 3748 | at_y, & |
---|
| 3749 | bc_mas_lr, & |
---|
| 3750 | bc_mas_ns, & |
---|
| 3751 | coll_t_0, & |
---|
| 3752 | corner_gate_start, & |
---|
| 3753 | corner_gate_width, & |
---|
[3235] | 3754 | dim_size_agtnum_manual, & |
---|
| 3755 | dim_size_factor_agtnum, & |
---|
[3198] | 3756 | deallocate_memory_mas, & |
---|
| 3757 | dist_to_int_target, & |
---|
| 3758 | dt_agent, & |
---|
| 3759 | dt_arel, & |
---|
| 3760 | dt_write_agent_data, & |
---|
| 3761 | end_time_arel, & |
---|
| 3762 | max_dist_from_path, & |
---|
| 3763 | min_nr_agent, & |
---|
| 3764 | multi_agent_system_end, & |
---|
| 3765 | multi_agent_system_start, & |
---|
| 3766 | number_of_agent_groups, & |
---|
| 3767 | radius_agent, & |
---|
| 3768 | random_start_position_agents, & |
---|
| 3769 | read_agents_from_restartfile, & |
---|
| 3770 | repuls_agent, & |
---|
| 3771 | repuls_wall, & |
---|
| 3772 | scan_radius_agent, & |
---|
| 3773 | sigma_rep_agent, & |
---|
| 3774 | sigma_rep_wall, & |
---|
| 3775 | step_dealloc_mas, & |
---|
| 3776 | tau_accel_agent |
---|
[3159] | 3777 | |
---|
| 3778 | ! |
---|
| 3779 | !-- Try to find agent package |
---|
| 3780 | REWIND ( 11 ) |
---|
| 3781 | line = ' ' |
---|
[3248] | 3782 | DO WHILE ( INDEX( line, '&agent_parameters' ) == 0 ) |
---|
[3159] | 3783 | READ ( 11, '(A)', END=20 ) line |
---|
| 3784 | ENDDO |
---|
| 3785 | BACKSPACE ( 11 ) |
---|
| 3786 | |
---|
| 3787 | ! |
---|
| 3788 | !-- Read user-defined namelist |
---|
[3246] | 3789 | READ ( 11, agent_parameters, ERR = 10, END = 20 ) |
---|
[3159] | 3790 | |
---|
| 3791 | ! |
---|
| 3792 | !-- Set flag that indicates that agents are switched on |
---|
| 3793 | agents_active = .TRUE. |
---|
[3246] | 3794 | GOTO 20 |
---|
[3159] | 3795 | |
---|
[3246] | 3796 | 10 BACKSPACE( 11 ) |
---|
[3248] | 3797 | READ( 11 , '(A)') line |
---|
| 3798 | CALL parin_fail_message( 'agent_parameters', line ) |
---|
[3246] | 3799 | |
---|
[3159] | 3800 | 20 CONTINUE |
---|
| 3801 | |
---|
| 3802 | END SUBROUTINE mas_parin |
---|
| 3803 | |
---|
| 3804 | !------------------------------------------------------------------------------! |
---|
| 3805 | ! Description: |
---|
| 3806 | ! ------------ |
---|
| 3807 | !> Routine for the whole processor |
---|
| 3808 | !> Sort all agents into the 4 respective subgrid boxes |
---|
| 3809 | !------------------------------------------------------------------------------! |
---|
| 3810 | SUBROUTINE mas_ps_sort_in_subboxes |
---|
| 3811 | |
---|
| 3812 | IMPLICIT NONE |
---|
| 3813 | |
---|
| 3814 | INTEGER(iwp) :: i !< grid box (x) |
---|
| 3815 | INTEGER(iwp) :: ip !< counter (x) |
---|
| 3816 | INTEGER(iwp) :: is !< box counter |
---|
| 3817 | INTEGER(iwp) :: j !< grid box (y) |
---|
| 3818 | INTEGER(iwp) :: jp !< counter (y) |
---|
| 3819 | INTEGER(iwp) :: m !< sorting index |
---|
| 3820 | INTEGER(iwp) :: n !< agent index |
---|
| 3821 | INTEGER(iwp) :: nn !< agent counter |
---|
| 3822 | INTEGER(iwp) :: sort_index !< sorting index |
---|
| 3823 | |
---|
| 3824 | INTEGER(iwp), DIMENSION(0:3) :: sort_count !< number of agents in one subbox |
---|
| 3825 | |
---|
| 3826 | TYPE(agent_type), DIMENSION(:,:), ALLOCATABLE :: sort_agents !< sorted agent array |
---|
| 3827 | |
---|
| 3828 | DO ip = nxl, nxr |
---|
| 3829 | DO jp = nys, nyn |
---|
| 3830 | number_of_agents = agt_count(jp,ip) |
---|
| 3831 | IF ( number_of_agents <= 0 ) CYCLE |
---|
| 3832 | agents => grid_agents(jp,ip)%agents(1:number_of_agents) |
---|
| 3833 | |
---|
| 3834 | nn = 0 |
---|
| 3835 | sort_count = 0 |
---|
| 3836 | ALLOCATE( sort_agents(number_of_agents, 0:3) ) |
---|
| 3837 | |
---|
| 3838 | DO n = 1, number_of_agents |
---|
| 3839 | sort_index = 0 |
---|
| 3840 | |
---|
| 3841 | IF ( agents(n)%agent_mask ) THEN |
---|
| 3842 | nn = nn + 1 |
---|
| 3843 | ! |
---|
| 3844 | !-- Sorting agents with a binary scheme |
---|
| 3845 | !-- sort_index=11_2=3_10 -> agent at the left,south subgridbox |
---|
| 3846 | !-- sort_index=10_2=2_10 -> agent at the left,north subgridbox |
---|
| 3847 | !-- sort_index=01_2=1_10 -> agent at the right,south subgridbox |
---|
| 3848 | !-- sort_index=00_2=0_10 -> agent at the right,north subgridbox |
---|
| 3849 | !-- For this the center of the gridbox is calculated |
---|
| 3850 | i = (agents(n)%x + 0.5_wp * dx) * ddx |
---|
| 3851 | j = (agents(n)%y + 0.5_wp * dy) * ddy |
---|
| 3852 | |
---|
| 3853 | IF ( i == ip ) sort_index = sort_index + 2 |
---|
| 3854 | IF ( j == jp ) sort_index = sort_index + 1 |
---|
| 3855 | |
---|
| 3856 | sort_count(sort_index) = sort_count(sort_index) + 1 |
---|
| 3857 | m = sort_count(sort_index) |
---|
| 3858 | sort_agents(m,sort_index) = agents(n) |
---|
| 3859 | sort_agents(m,sort_index)%block_nr = sort_index |
---|
| 3860 | ENDIF |
---|
| 3861 | ENDDO |
---|
| 3862 | |
---|
| 3863 | nn = 0 |
---|
| 3864 | DO is = 0,3 |
---|
| 3865 | grid_agents(jp,ip)%start_index(is) = nn + 1 |
---|
| 3866 | DO n = 1,sort_count(is) |
---|
| 3867 | nn = nn + 1 |
---|
| 3868 | agents(nn) = sort_agents(n,is) |
---|
| 3869 | ENDDO |
---|
| 3870 | grid_agents(jp,ip)%end_index(is) = nn |
---|
| 3871 | ENDDO |
---|
| 3872 | |
---|
| 3873 | number_of_agents = nn |
---|
| 3874 | agt_count(jp,ip) = number_of_agents |
---|
| 3875 | DEALLOCATE(sort_agents) |
---|
| 3876 | ENDDO |
---|
| 3877 | ENDDO |
---|
| 3878 | |
---|
| 3879 | END SUBROUTINE mas_ps_sort_in_subboxes |
---|
| 3880 | |
---|
| 3881 | !------------------------------------------------------------------------------! |
---|
| 3882 | ! Description: |
---|
| 3883 | ! ------------ |
---|
| 3884 | !> Move all agents not marked for deletion to lowest indices (packing) |
---|
| 3885 | !------------------------------------------------------------------------------! |
---|
| 3886 | SUBROUTINE mas_ps_pack |
---|
| 3887 | |
---|
| 3888 | IMPLICIT NONE |
---|
| 3889 | |
---|
| 3890 | INTEGER(iwp) :: n !< agent counter |
---|
| 3891 | INTEGER(iwp) :: nn !< number of agents |
---|
| 3892 | ! |
---|
| 3893 | !-- Find out elements marked for deletion and move data from highest index |
---|
| 3894 | !-- values to these free indices |
---|
| 3895 | nn = number_of_agents |
---|
| 3896 | |
---|
| 3897 | DO WHILE ( .NOT. agents(nn)%agent_mask ) |
---|
| 3898 | nn = nn-1 |
---|
| 3899 | IF ( nn == 0 ) EXIT |
---|
| 3900 | ENDDO |
---|
| 3901 | |
---|
| 3902 | IF ( nn > 0 ) THEN |
---|
| 3903 | DO n = 1, number_of_agents |
---|
| 3904 | IF ( .NOT. agents(n)%agent_mask ) THEN |
---|
| 3905 | agents(n) = agents(nn) |
---|
| 3906 | nn = nn - 1 |
---|
| 3907 | DO WHILE ( .NOT. agents(nn)%agent_mask ) |
---|
| 3908 | nn = nn-1 |
---|
| 3909 | IF ( n == nn ) EXIT |
---|
| 3910 | ENDDO |
---|
| 3911 | ENDIF |
---|
| 3912 | IF ( n == nn ) EXIT |
---|
| 3913 | ENDDO |
---|
| 3914 | ENDIF |
---|
| 3915 | |
---|
| 3916 | ! |
---|
| 3917 | !-- The number of deleted agents has been determined in routines |
---|
| 3918 | !-- mas_boundary_conds, mas_droplet_collision, and mas_eh_exchange_horiz |
---|
| 3919 | number_of_agents = nn |
---|
| 3920 | |
---|
| 3921 | END SUBROUTINE mas_ps_pack |
---|
| 3922 | |
---|
| 3923 | !------------------------------------------------------------------------------! |
---|
| 3924 | ! Description: |
---|
| 3925 | ! ------------ |
---|
| 3926 | !> Sort agents in each sub-grid box into two groups: agents that already |
---|
| 3927 | !> completed the LES timestep, and agents that need further timestepping to |
---|
| 3928 | !> complete the LES timestep. |
---|
| 3929 | !------------------------------------------------------------------------------! |
---|
[3241] | 3930 | ! SUBROUTINE mas_ps_sort_timeloop_done |
---|
[3159] | 3931 | ! |
---|
[3241] | 3932 | ! IMPLICIT NONE |
---|
| 3933 | ! |
---|
| 3934 | ! INTEGER(iwp) :: end_index !< agent end index for each sub-box |
---|
| 3935 | ! INTEGER(iwp) :: i !< index of agent grid box in x-direction |
---|
| 3936 | ! INTEGER(iwp) :: j !< index of agent grid box in y-direction |
---|
| 3937 | ! INTEGER(iwp) :: n !< running index for number of agents |
---|
| 3938 | ! INTEGER(iwp) :: nb !< index of subgrid boux |
---|
| 3939 | ! INTEGER(iwp) :: nf !< indices for agents in each sub-box that already finalized their substeps |
---|
| 3940 | ! INTEGER(iwp) :: nnf !< indices for agents in each sub-box that need further treatment |
---|
| 3941 | ! INTEGER(iwp) :: num_finalized !< number of agents in each sub-box that already finalized their substeps |
---|
| 3942 | ! INTEGER(iwp) :: start_index !< agent start index for each sub-box |
---|
| 3943 | ! |
---|
| 3944 | ! TYPE(agent_type), DIMENSION(:), ALLOCATABLE :: sort_agents !< temporary agent array |
---|
| 3945 | ! |
---|
| 3946 | ! DO i = nxl, nxr |
---|
| 3947 | ! DO j = nys, nyn |
---|
| 3948 | ! |
---|
| 3949 | ! number_of_agents = agt_count(j,i) |
---|
| 3950 | ! IF ( number_of_agents <= 0 ) CYCLE |
---|
| 3951 | ! |
---|
| 3952 | ! agents => grid_agents(j,i)%agents(1:number_of_agents) |
---|
| 3953 | ! |
---|
| 3954 | ! DO nb = 0, 3 |
---|
| 3955 | ! |
---|
[3159] | 3956 | !-- Obtain start and end index for each subgrid box |
---|
[3241] | 3957 | ! start_index = grid_agents(j,i)%start_index(nb) |
---|
| 3958 | ! end_index = grid_agents(j,i)%end_index(nb) |
---|
[3159] | 3959 | ! |
---|
| 3960 | !-- Allocate temporary array used for sorting |
---|
[3241] | 3961 | ! ALLOCATE( sort_agents(start_index:end_index) ) |
---|
[3159] | 3962 | ! |
---|
| 3963 | !-- Determine number of agents already completed the LES |
---|
| 3964 | !-- timestep, and write them into a temporary array |
---|
[3241] | 3965 | ! nf = start_index |
---|
| 3966 | ! num_finalized = 0 |
---|
| 3967 | ! DO n = start_index, end_index |
---|
| 3968 | ! IF ( dt_3d - agents(n)%dt_sum < 1E-8_wp ) THEN |
---|
| 3969 | ! sort_agents(nf) = agents(n) |
---|
| 3970 | ! nf = nf + 1 |
---|
| 3971 | ! num_finalized = num_finalized + 1 |
---|
| 3972 | ! ENDIF |
---|
| 3973 | ! ENDDO |
---|
[3159] | 3974 | ! |
---|
| 3975 | !-- Determine number of agents that not completed the LES |
---|
| 3976 | !-- timestep, and write them into a temporary array |
---|
[3241] | 3977 | ! nnf = nf |
---|
| 3978 | ! DO n = start_index, end_index |
---|
| 3979 | ! IF ( dt_3d - agents(n)%dt_sum > 1E-8_wp ) THEN |
---|
| 3980 | ! sort_agents(nnf) = agents(n) |
---|
| 3981 | ! nnf = nnf + 1 |
---|
| 3982 | ! ENDIF |
---|
| 3983 | ! ENDDO |
---|
[3159] | 3984 | ! |
---|
| 3985 | !-- Write back sorted agents |
---|
[3241] | 3986 | ! agents(start_index:end_index) = & |
---|
| 3987 | ! sort_agents(start_index:end_index) |
---|
[3159] | 3988 | ! |
---|
| 3989 | !-- Determine updated start_index, used to masked already |
---|
| 3990 | !-- completed agents. |
---|
[3241] | 3991 | ! grid_agents(j,i)%start_index(nb) = & |
---|
| 3992 | ! grid_agents(j,i)%start_index(nb) & |
---|
| 3993 | ! + num_finalized |
---|
[3159] | 3994 | ! |
---|
| 3995 | !-- Deallocate dummy array |
---|
[3241] | 3996 | ! DEALLOCATE ( sort_agents ) |
---|
[3159] | 3997 | ! |
---|
| 3998 | !-- Finally, if number of non-completed agents is non zero |
---|
| 3999 | !-- in any of the sub-boxes, set control flag appropriately. |
---|
[3241] | 4000 | ! IF ( nnf > nf ) & |
---|
| 4001 | ! grid_agents(j,i)%time_loop_done = .FALSE. |
---|
| 4002 | ! |
---|
| 4003 | ! ENDDO |
---|
| 4004 | ! ENDDO |
---|
| 4005 | ! ENDDO |
---|
| 4006 | ! |
---|
| 4007 | ! END SUBROUTINE mas_ps_sort_timeloop_done |
---|
[3159] | 4008 | |
---|
| 4009 | !------------------------------------------------------------------------------! |
---|
| 4010 | ! Description: |
---|
| 4011 | ! ------------ |
---|
| 4012 | !> Calls social forces calculations |
---|
| 4013 | !------------------------------------------------------------------------------! |
---|
| 4014 | SUBROUTINE mas_timestep_forces_call ( ip, jp ) |
---|
| 4015 | |
---|
| 4016 | IMPLICIT NONE |
---|
| 4017 | |
---|
| 4018 | INTEGER(iwp) :: ip !< counter, x-direction |
---|
| 4019 | INTEGER(iwp) :: jp !< counter, y-direction |
---|
| 4020 | INTEGER(iwp) :: n !< loop variable over all agents in a grid box |
---|
| 4021 | |
---|
| 4022 | ! |
---|
| 4023 | !-- Get direction for all agents in current grid cell |
---|
| 4024 | CALL mas_agent_direction |
---|
| 4025 | |
---|
| 4026 | DO n = 1, number_of_agents |
---|
| 4027 | |
---|
| 4028 | force_x = 0.0_wp |
---|
| 4029 | force_y = 0.0_wp |
---|
| 4030 | |
---|
| 4031 | CALL mas_timestep_social_forces ( 'acceleration', n, ip, jp ) |
---|
| 4032 | |
---|
| 4033 | CALL mas_timestep_social_forces ( 'other_agents', n, ip, jp ) |
---|
| 4034 | |
---|
| 4035 | CALL mas_timestep_social_forces ( 'walls', n, ip, jp ) |
---|
| 4036 | ! |
---|
| 4037 | !-- Update forces |
---|
| 4038 | agents(n)%force_x = force_x |
---|
| 4039 | agents(n)%force_y = force_y |
---|
| 4040 | ENDDO |
---|
| 4041 | |
---|
| 4042 | END SUBROUTINE mas_timestep_forces_call |
---|
| 4043 | |
---|
| 4044 | !------------------------------------------------------------------------------! |
---|
| 4045 | ! Description: |
---|
| 4046 | ! ------------ |
---|
| 4047 | !> Euler timestep of agent transport |
---|
| 4048 | !------------------------------------------------------------------------------! |
---|
| 4049 | SUBROUTINE mas_timestep |
---|
| 4050 | |
---|
| 4051 | IMPLICIT NONE |
---|
| 4052 | |
---|
| 4053 | INTEGER(iwp) :: n !< loop variable over all agents in a grid box |
---|
| 4054 | |
---|
| 4055 | REAL(wp) :: abs_v !< absolute value of velocity |
---|
| 4056 | REAL(wp) :: abs_f !< absolute value of force |
---|
| 4057 | |
---|
| 4058 | DO n = 1, number_of_agents |
---|
| 4059 | ! |
---|
| 4060 | !-- Limit absolute force to a maximum to prevent unrealistic acceleration |
---|
| 4061 | abs_f = SQRT((agents(n)%force_x)**2 + (agents(n)%force_y)**2) |
---|
| 4062 | IF ( abs_f > 20. ) THEN |
---|
| 4063 | agents(n)%force_x = agents(n)%force_x * 20. / abs_f |
---|
| 4064 | agents(n)%force_y = agents(n)%force_y * 20. / abs_f |
---|
| 4065 | ENDIF |
---|
| 4066 | ! |
---|
| 4067 | !-- Update agent speed |
---|
| 4068 | agents(n)%speed_x = agents(n)%speed_x + agents(n)%force_x * dt_agent |
---|
| 4069 | agents(n)%speed_y = agents(n)%speed_y + agents(n)%force_y * dt_agent |
---|
| 4070 | ! |
---|
| 4071 | !-- Reduction of agent speed to maximum agent speed |
---|
| 4072 | abs_v = SQRT((agents(n)%speed_x)**2 + (agents(n)%speed_y)**2) |
---|
| 4073 | IF ( abs_v > v_max_agent ) THEN |
---|
| 4074 | agents(n)%speed_x = agents(n)%speed_x * v_max_agent / abs_v |
---|
| 4075 | agents(n)%speed_y = agents(n)%speed_y * v_max_agent / abs_v |
---|
| 4076 | ENDIF |
---|
| 4077 | ! |
---|
| 4078 | !-- Update agent position |
---|
| 4079 | agents(n)%x = agents(n)%x + agents(n)%speed_x * dt_agent |
---|
| 4080 | agents(n)%y = agents(n)%y + agents(n)%speed_y * dt_agent |
---|
| 4081 | ! |
---|
| 4082 | !-- Update absolute value of agent speed |
---|
| 4083 | agents(n)%speed_abs = abs_v |
---|
| 4084 | ! |
---|
| 4085 | !-- Increment the agent age and the total time that the agent |
---|
| 4086 | !-- has advanced within the agent timestep procedure |
---|
| 4087 | agents(n)%age_m = agents(n)%age |
---|
| 4088 | agents(n)%age = agents(n)%age + dt_agent |
---|
| 4089 | agents(n)%dt_sum = agents(n)%dt_sum + dt_agent |
---|
| 4090 | ! |
---|
| 4091 | !-- Check whether there is still an agent that has not yet completed |
---|
| 4092 | !-- the total LES timestep |
---|
| 4093 | IF ( ( dt_3d - agents(n)%dt_sum ) > 1E-8_wp ) THEN |
---|
| 4094 | dt_3d_reached_l_mas = .FALSE. |
---|
| 4095 | ENDIF |
---|
| 4096 | |
---|
| 4097 | ENDDO |
---|
| 4098 | |
---|
| 4099 | END SUBROUTINE mas_timestep |
---|
| 4100 | |
---|
| 4101 | !------------------------------------------------------------------------------! |
---|
| 4102 | ! Description: |
---|
| 4103 | ! ------------ |
---|
| 4104 | !> Calculates the Social Forces (Helbing and Molnar, 1995) that the agent |
---|
| 4105 | !> experiences due to acceleration towards target and repulsion by obstacles |
---|
| 4106 | !------------------------------------------------------------------------------! |
---|
| 4107 | SUBROUTINE mas_timestep_social_forces ( mode, nl, ip, jp ) |
---|
| 4108 | |
---|
[3241] | 4109 | ! USE constants, & |
---|
| 4110 | ! ONLY: pi |
---|
[3159] | 4111 | |
---|
| 4112 | IMPLICIT NONE |
---|
| 4113 | |
---|
| 4114 | CHARACTER (LEN=*) :: mode !< identifier for the mode of calculation |
---|
| 4115 | |
---|
| 4116 | INTEGER(iwp) :: ij_dum !< index of nearest wall |
---|
| 4117 | INTEGER(iwp) :: il !< index variable along x |
---|
| 4118 | INTEGER(iwp) :: ip !< index variable along x |
---|
| 4119 | INTEGER(iwp) :: jl !< index variable along y |
---|
| 4120 | INTEGER(iwp) :: jp !< index variable along y |
---|
| 4121 | INTEGER(iwp) :: nl !< loop variable over all agents in a grid box |
---|
| 4122 | INTEGER(iwp) :: no !< loop variable over all agents in a grid box |
---|
| 4123 | INTEGER(iwp) :: noa !< amount of agents in a grid box |
---|
| 4124 | INTEGER(iwp) :: sc_x_end !< index for scan for topography/other agents |
---|
| 4125 | INTEGER(iwp) :: sc_x_start !< index for scan for topography/other agents |
---|
| 4126 | INTEGER(iwp) :: sc_y_end !< index for scan for topography/other agents |
---|
| 4127 | INTEGER(iwp) :: sc_y_start !< index for scan for topography/other agents |
---|
| 4128 | |
---|
| 4129 | LOGICAL :: corner_found !< flag that indicates a corner has been found near agent |
---|
| 4130 | |
---|
| 4131 | REAL(wp) :: a_pl !< factor for collision avoidance |
---|
| 4132 | REAL(wp) :: ax_semimaj !< semiminor axis of repulsive ellipse |
---|
| 4133 | REAL(wp) :: b_pl !< factor for collision avoidance |
---|
| 4134 | REAL(wp) :: c_pl !< factor for collision avoidance |
---|
| 4135 | REAL(wp) :: coll_t !< time at which the next collision would happen |
---|
| 4136 | REAL(wp) :: d_coll_t_0 !< inverse of collision cutoff time |
---|
| 4137 | REAL(wp) :: d_pl !< factor for collision avoidance |
---|
| 4138 | REAL(wp) :: ddum_f !< dummy devisor collision avoidance |
---|
| 4139 | REAL(wp) :: dist !< distance to obstacle |
---|
| 4140 | REAL(wp) :: dist_sq !< distance to obstacle squared |
---|
| 4141 | REAL(wp) :: pos_rel_x !< relative position of two agents (x) |
---|
| 4142 | REAL(wp) :: pos_rel_y !< relative position of two agents (y) |
---|
| 4143 | REAL(wp) :: r_sq !< y-position |
---|
| 4144 | REAL(wp) :: sra !< scan radius (agents) |
---|
| 4145 | REAL(wp) :: srw !< local variable for scan radius (walls) |
---|
| 4146 | REAL(wp) :: v_rel_x !< relative velocity (x); collision avoidance |
---|
| 4147 | REAL(wp) :: v_rel_y !< relative velocity (y); collision avoidance |
---|
| 4148 | REAL(wp) :: x_a !< x-position |
---|
| 4149 | REAL(wp) :: x_wall !< x-position of wall |
---|
| 4150 | REAL(wp) :: y_a !< y-position |
---|
| 4151 | REAL(wp) :: y_wall !< y-position of wall |
---|
| 4152 | |
---|
| 4153 | REAL(wp), PARAMETER :: k_pl = 1.5 !< factor for collision avoidance |
---|
| 4154 | |
---|
| 4155 | TYPE(agent_type), DIMENSION(:), POINTER :: l_agts !< agents that repulse current agent |
---|
| 4156 | |
---|
| 4157 | ! |
---|
| 4158 | !-- Initialization |
---|
| 4159 | x_a = agents(nl)%x |
---|
| 4160 | y_a = agents(nl)%y |
---|
| 4161 | |
---|
| 4162 | SELECT CASE ( TRIM( mode ) ) |
---|
| 4163 | ! |
---|
| 4164 | !-- Calculation of force due to agent trying to approach desired velocity |
---|
| 4165 | CASE ( 'acceleration' ) |
---|
| 4166 | |
---|
| 4167 | force_x = force_x + d_tau_accel_agent & |
---|
| 4168 | * ( agents(nl)%speed_des*agents(nl)%speed_e_x & |
---|
| 4169 | -agents(nl)%speed_x ) |
---|
| 4170 | |
---|
| 4171 | force_y = force_y + d_tau_accel_agent & |
---|
| 4172 | * ( agents(nl)%speed_des*agents(nl)%speed_e_y & |
---|
| 4173 | -agents(nl)%speed_y ) |
---|
| 4174 | |
---|
| 4175 | ! |
---|
| 4176 | !-- Calculation of repulsive forces by other agents in a radius around the |
---|
| 4177 | !-- current one |
---|
| 4178 | CASE ( 'other_agents' ) |
---|
| 4179 | |
---|
| 4180 | sra = scan_radius_agent |
---|
| 4181 | d_coll_t_0 = 1./coll_t_0 |
---|
| 4182 | ! |
---|
| 4183 | !-- Find relevant gridboxes (those that could contain agents within |
---|
| 4184 | !-- scan radius) |
---|
| 4185 | sc_x_start = FLOOR( (x_a - sra) * ddx ) |
---|
| 4186 | sc_x_end = FLOOR( (x_a + sra) * ddx ) |
---|
| 4187 | sc_y_start = FLOOR( (y_a - sra) * ddx ) |
---|
| 4188 | sc_y_end = FLOOR( (y_a + sra) * ddx ) |
---|
| 4189 | IF ( sc_x_start < nxlg ) sc_x_start = nxlg |
---|
| 4190 | IF ( sc_x_end > nxrg ) sc_x_end = nxrg |
---|
| 4191 | IF ( sc_y_start < nysg ) sc_y_start = nysg |
---|
| 4192 | IF ( sc_y_end > nyng ) sc_y_end = nyng |
---|
| 4193 | |
---|
| 4194 | sra = sra**2 |
---|
| 4195 | ! |
---|
| 4196 | !-- Loop over all previously found relevant gridboxes |
---|
| 4197 | DO il = sc_x_start, sc_x_end |
---|
| 4198 | DO jl = sc_y_start, sc_y_end |
---|
| 4199 | noa = agt_count(jl,il) |
---|
| 4200 | IF ( noa <= 0 ) CYCLE |
---|
| 4201 | l_agts => grid_agents(jl,il)%agents(1:noa) |
---|
| 4202 | DO no = 1, noa |
---|
| 4203 | ! |
---|
| 4204 | !-- Skip self |
---|
| 4205 | IF ( jl == jp .AND. il == ip .AND. no == nl ) CYCLE |
---|
| 4206 | pos_rel_x = l_agts(no)%x - x_a |
---|
| 4207 | pos_rel_y = l_agts(no)%y - y_a |
---|
| 4208 | dist_sq = pos_rel_x**2 + pos_rel_y**2 |
---|
| 4209 | IF ( dist_sq > sra ) CYCLE |
---|
| 4210 | r_sq = (2*radius_agent)**2 |
---|
| 4211 | v_rel_x = agents(nl)%speed_x - l_agts(no)%speed_x |
---|
| 4212 | v_rel_y = agents(nl)%speed_y - l_agts(no)%speed_y |
---|
| 4213 | ! |
---|
| 4214 | !-- Collision is already occuring, default to standard |
---|
| 4215 | !-- social forces |
---|
| 4216 | IF ( dist_sq <= r_sq ) THEN |
---|
| 4217 | dist = SQRT(dist_sq) + 1.0d-12 |
---|
| 4218 | ax_semimaj = .5_wp*SQRT( dist ) |
---|
| 4219 | |
---|
| 4220 | force_x = force_x - 0.125_wp * repuls_agent & |
---|
| 4221 | * d_sigma_rep_agent / ax_semimaj & |
---|
| 4222 | * EXP( -ax_semimaj*d_sigma_rep_agent ) & |
---|
| 4223 | * (pos_rel_x/dist) |
---|
| 4224 | |
---|
| 4225 | force_y = force_y - 0.125_wp * repuls_agent & |
---|
| 4226 | * d_sigma_rep_agent / ax_semimaj & |
---|
| 4227 | * EXP( -ax_semimaj*d_sigma_rep_agent ) & |
---|
| 4228 | * (pos_rel_y/dist) |
---|
| 4229 | ! |
---|
| 4230 | !-- Currently no collision, calculate collision avoidance |
---|
| 4231 | !-- force according to Karamouzas et al (2014, PRL 113,238701) |
---|
| 4232 | ELSE |
---|
| 4233 | ! |
---|
| 4234 | !-- factors |
---|
| 4235 | a_pl = v_rel_x**2 + v_rel_y**2 |
---|
| 4236 | b_pl = pos_rel_x*v_rel_x + pos_rel_y*v_rel_y |
---|
| 4237 | c_pl = dist_sq - r_sq |
---|
| 4238 | d_pl = b_pl**2 - a_pl*c_pl |
---|
| 4239 | ! |
---|
| 4240 | !-- If the two agents are moving non-parallel, calculate |
---|
| 4241 | !-- collision avoidance social force |
---|
| 4242 | IF ( d_pl > 0.0_wp .AND. & |
---|
| 4243 | ( a_pl < -0.00001 .OR. a_pl > 0.00001 ) ) & |
---|
| 4244 | THEN |
---|
| 4245 | |
---|
| 4246 | d_pl = SQRT(d_pl) |
---|
| 4247 | coll_t = (b_pl - d_pl)/a_pl |
---|
| 4248 | IF ( coll_t > 0.0_wp ) THEN |
---|
| 4249 | ! |
---|
| 4250 | !-- Dummy factor |
---|
| 4251 | ddum_f = 1. / ( a_pl * coll_t**2 ) & |
---|
| 4252 | * ( 2. / coll_t + 1.0 * d_coll_t_0 ) |
---|
| 4253 | ! |
---|
| 4254 | !-- x-component of social force |
---|
| 4255 | force_x = force_x - k_pl * & |
---|
| 4256 | EXP( -coll_t * d_coll_t_0 ) * & |
---|
| 4257 | ( v_rel_x - & |
---|
| 4258 | ( b_pl * v_rel_x - & |
---|
| 4259 | a_pl * pos_rel_x ) / d_pl ) * & |
---|
| 4260 | ddum_f |
---|
| 4261 | ! |
---|
| 4262 | !-- y-component of social force |
---|
| 4263 | force_y = force_y - k_pl * & |
---|
| 4264 | EXP( -coll_t * d_coll_t_0 ) * & |
---|
| 4265 | ( v_rel_y - & |
---|
| 4266 | ( b_pl * v_rel_y - & |
---|
| 4267 | a_pl * pos_rel_y ) / d_pl ) * & |
---|
| 4268 | ddum_f |
---|
| 4269 | |
---|
| 4270 | ENDIF |
---|
| 4271 | ENDIF |
---|
| 4272 | ENDIF |
---|
| 4273 | ENDDO |
---|
| 4274 | ENDDO |
---|
| 4275 | ENDDO |
---|
| 4276 | |
---|
| 4277 | CASE ( 'walls' ) |
---|
| 4278 | |
---|
| 4279 | srw = scan_radius_wall |
---|
| 4280 | corner_found = .FALSE. |
---|
| 4281 | ! |
---|
| 4282 | !-- find relevant grid boxes (those that could contain topography |
---|
| 4283 | !-- within radius) |
---|
| 4284 | sc_x_start = (x_a - srw) * ddx |
---|
| 4285 | sc_x_end = (x_a + srw) * ddx |
---|
| 4286 | sc_y_start = (y_a - srw) * ddx |
---|
| 4287 | sc_y_end = (y_a + srw) * ddx |
---|
| 4288 | IF ( sc_x_start < nxlg ) sc_x_start = nxlg |
---|
| 4289 | IF ( sc_x_end > nxrg ) sc_x_end = nxrg |
---|
| 4290 | IF ( sc_y_start < nysg ) sc_y_start = nysg |
---|
| 4291 | IF ( sc_y_end > nyng ) sc_y_end = nyng |
---|
| 4292 | ! |
---|
| 4293 | !-- Find "walls" ( i.e. topography steps (up or down) higher than one |
---|
| 4294 | !-- grid box ) that are perpendicular to the agent within the defined |
---|
| 4295 | !-- search radius. Such obstacles cannot be passed and a social force |
---|
| 4296 | !-- to that effect is applied. |
---|
| 4297 | !-- Walls only apply a force perpendicular to the wall to the agent. |
---|
| 4298 | !-- There is therefore a search for walls directly right, left, south |
---|
| 4299 | !-- and north of the agent. All other walls are ignored. |
---|
| 4300 | !-- |
---|
| 4301 | !-- Check for wall left of current agent |
---|
| 4302 | ij_dum = 0 |
---|
| 4303 | IF ( sc_x_start < ip ) THEN |
---|
| 4304 | DO il = ip - 1, sc_x_start, -1 |
---|
| 4305 | ! |
---|
| 4306 | !-- Going left from the agent, check for a right wall |
---|
| 4307 | IF ( BTEST( obstacle_flags(jp,il), 2 ) ) THEN |
---|
| 4308 | ! |
---|
| 4309 | !-- obstacle found in grid box il, wall at right side |
---|
| 4310 | x_wall = (il+1)*dx |
---|
| 4311 | ! |
---|
| 4312 | !-- Calculate force of found wall on agent |
---|
[3241] | 4313 | CALL mas_timestep_wall_corner_force( x_a, x_wall, y_a, & |
---|
| 4314 | y_a ) |
---|
[3159] | 4315 | ! |
---|
| 4316 | !-- calculate new x starting index for later scan for corners |
---|
| 4317 | ij_dum = il + 1 |
---|
| 4318 | EXIT |
---|
| 4319 | ENDIF |
---|
| 4320 | ENDDO |
---|
| 4321 | ENDIF |
---|
| 4322 | IF ( ij_dum /= 0 ) sc_x_start = ij_dum |
---|
| 4323 | |
---|
| 4324 | ! |
---|
| 4325 | !-- Check for wall right of current agent |
---|
| 4326 | ij_dum = 0 |
---|
| 4327 | IF ( sc_x_end > ip ) THEN |
---|
| 4328 | DO il = ip + 1, sc_x_end |
---|
| 4329 | ! |
---|
| 4330 | !-- Going right from the agent, check for a left wall |
---|
| 4331 | IF ( BTEST( obstacle_flags(jp,il), 6 ) ) THEN |
---|
| 4332 | ! |
---|
| 4333 | !-- obstacle found in grid box il, wall at left side |
---|
| 4334 | x_wall = il*dx |
---|
| 4335 | ! |
---|
| 4336 | !-- Calculate force of found wall on agent |
---|
[3241] | 4337 | CALL mas_timestep_wall_corner_force( x_a, x_wall, y_a, & |
---|
| 4338 | y_a ) |
---|
[3159] | 4339 | ! |
---|
| 4340 | !-- calculate new x end index for later scan for corners |
---|
| 4341 | ij_dum = il - 1 |
---|
| 4342 | EXIT |
---|
| 4343 | ENDIF |
---|
| 4344 | ENDDO |
---|
| 4345 | ENDIF |
---|
| 4346 | IF ( ij_dum /= 0 ) sc_x_end = ij_dum |
---|
| 4347 | |
---|
| 4348 | ! |
---|
| 4349 | !-- Check for wall south of current agent |
---|
| 4350 | ij_dum = 0 |
---|
| 4351 | IF ( sc_y_start < jp ) THEN |
---|
| 4352 | DO jl = jp - 1, sc_y_start, -1 |
---|
| 4353 | ! |
---|
| 4354 | !-- Going south from the agent, check for a north wall |
---|
| 4355 | IF ( BTEST( obstacle_flags(jl,ip), 0 ) ) THEN |
---|
| 4356 | ! |
---|
| 4357 | !-- obstacle found in grid box jl, wall at left side |
---|
| 4358 | y_wall = (jl+1)*dy |
---|
| 4359 | |
---|
[3241] | 4360 | CALL mas_timestep_wall_corner_force( x_a, x_a, y_a, & |
---|
| 4361 | y_wall ) |
---|
[3159] | 4362 | ! |
---|
| 4363 | !-- calculate new y starting index for later scan for corners |
---|
| 4364 | ij_dum = jl + 1 |
---|
| 4365 | EXIT |
---|
| 4366 | ENDIF |
---|
| 4367 | ENDDO |
---|
| 4368 | ENDIF |
---|
| 4369 | IF ( ij_dum /= 0 ) sc_y_start = ij_dum |
---|
| 4370 | |
---|
| 4371 | ! |
---|
| 4372 | !-- Check for wall north of current agent |
---|
| 4373 | ij_dum = 0 |
---|
| 4374 | IF ( sc_y_end > jp ) THEN |
---|
| 4375 | DO jl = jp + 1, sc_y_end |
---|
| 4376 | ! |
---|
| 4377 | !-- Going north from the agent, check for a south wall |
---|
| 4378 | IF ( BTEST( obstacle_flags(jl,ip), 4 ) ) THEN |
---|
| 4379 | ! |
---|
| 4380 | !-- obstacle found in grid box jl, wall at left side |
---|
| 4381 | y_wall = jl*dy |
---|
| 4382 | |
---|
[3241] | 4383 | CALL mas_timestep_wall_corner_force( x_a, x_a, y_a, & |
---|
| 4384 | y_wall ) |
---|
[3159] | 4385 | ! |
---|
| 4386 | !-- calculate new y end index for later scan for corners |
---|
| 4387 | ij_dum = jl - 1 |
---|
| 4388 | ENDIF |
---|
| 4389 | ENDDO |
---|
| 4390 | ENDIF |
---|
| 4391 | IF ( ij_dum /= 0 ) sc_y_end = ij_dum |
---|
| 4392 | |
---|
| 4393 | ! |
---|
| 4394 | !-- Scan for corners surrounding current agent. |
---|
| 4395 | !-- Only gridcells that are closer than the closest wall in each |
---|
| 4396 | !-- direction (n,s,r,l) are considered in the search since those |
---|
| 4397 | !-- further away would have a significantly smaller resulting force |
---|
| 4398 | !-- than the closer wall. |
---|
| 4399 | DO il = sc_x_start, sc_x_end |
---|
| 4400 | DO jl = sc_y_start, sc_y_end |
---|
| 4401 | IF ( il == ip .OR. jl == jp ) CYCLE |
---|
| 4402 | ! |
---|
| 4403 | !-- corners left of agent |
---|
| 4404 | IF ( il < ip ) THEN |
---|
| 4405 | ! |
---|
| 4406 | !-- south left quadrant: look for north right corner |
---|
| 4407 | IF ( jl < jp ) THEN |
---|
| 4408 | IF ( BTEST( obstacle_flags(jl,il), 1 ) ) THEN |
---|
| 4409 | ! |
---|
| 4410 | !-- calculate coordinates of the found corner |
---|
| 4411 | x_wall = (il+1)*dx |
---|
| 4412 | y_wall = (jl+1)*dy |
---|
| 4413 | |
---|
[3241] | 4414 | CALL mas_timestep_wall_corner_force( x_a, x_wall, & |
---|
| 4415 | y_a, y_wall ) |
---|
[3159] | 4416 | |
---|
| 4417 | ENDIF |
---|
| 4418 | ! |
---|
| 4419 | !-- north left quadrant: look for south right corner |
---|
| 4420 | ELSEIF ( jl > jp ) THEN |
---|
| 4421 | IF ( BTEST( obstacle_flags(jl,il), 3 ) ) THEN |
---|
| 4422 | ! |
---|
| 4423 | !-- calculate coordinates of the corner of said gridcell |
---|
| 4424 | !-- that is closest to the current agent |
---|
| 4425 | x_wall = (il+1)*dx |
---|
| 4426 | y_wall = jl*dy |
---|
| 4427 | |
---|
[3241] | 4428 | CALL mas_timestep_wall_corner_force( x_a, x_wall, & |
---|
| 4429 | y_a, y_wall ) |
---|
[3159] | 4430 | |
---|
| 4431 | ENDIF |
---|
| 4432 | ENDIF |
---|
| 4433 | ELSEIF ( il > ip ) THEN |
---|
| 4434 | ! |
---|
| 4435 | !-- south right quadrant: look for north left corner |
---|
| 4436 | IF ( jl < jp ) THEN |
---|
| 4437 | IF ( BTEST( obstacle_flags(jl,il), 7 ) ) THEN |
---|
| 4438 | ! |
---|
| 4439 | !-- calculate coordinates of the corner of said gridcell |
---|
| 4440 | !-- that is closest to the current agent |
---|
| 4441 | x_wall = il*dx |
---|
| 4442 | y_wall = (jl+1)*dy |
---|
| 4443 | |
---|
[3241] | 4444 | CALL mas_timestep_wall_corner_force( x_a, x_wall, & |
---|
| 4445 | y_a, y_wall ) |
---|
[3159] | 4446 | |
---|
| 4447 | ENDIF |
---|
| 4448 | ! |
---|
| 4449 | !-- north right quadrant: look for south left corner |
---|
| 4450 | ELSEIF ( jl > jp ) THEN |
---|
| 4451 | IF ( BTEST( obstacle_flags(jl,il), 5 ) ) THEN |
---|
| 4452 | ! |
---|
| 4453 | !-- calculate coordinates of the corner of said gridcell |
---|
| 4454 | !-- that is closest to the current agent |
---|
| 4455 | x_wall = il*dx |
---|
| 4456 | y_wall = jl*dy |
---|
| 4457 | |
---|
[3241] | 4458 | CALL mas_timestep_wall_corner_force( x_a, x_wall, & |
---|
| 4459 | y_a, y_wall ) |
---|
[3159] | 4460 | |
---|
| 4461 | ENDIF |
---|
| 4462 | ENDIF |
---|
| 4463 | ENDIF |
---|
| 4464 | ENDDO |
---|
| 4465 | ENDDO |
---|
| 4466 | |
---|
| 4467 | CASE DEFAULT |
---|
| 4468 | |
---|
| 4469 | END SELECT |
---|
| 4470 | |
---|
| 4471 | END SUBROUTINE mas_timestep_social_forces |
---|
| 4472 | |
---|
| 4473 | !------------------------------------------------------------------------------! |
---|
| 4474 | ! Description: |
---|
| 4475 | ! ------------ |
---|
| 4476 | !> Given a distance to the current agent, calculates the force a found corner |
---|
| 4477 | !> or wall exerts on that agent |
---|
| 4478 | !------------------------------------------------------------------------------! |
---|
[3241] | 4479 | SUBROUTINE mas_timestep_wall_corner_force( xa, xw, ya, yw ) |
---|
[3159] | 4480 | |
---|
| 4481 | IMPLICIT NONE |
---|
| 4482 | |
---|
| 4483 | REAL(wp) :: dist_l !< distance to obstacle |
---|
| 4484 | REAL(wp) :: force_d_x !< increment of social force, x-direction |
---|
| 4485 | REAL(wp) :: force_d_y !< increment of social force, x-direction |
---|
| 4486 | REAL(wp) :: xa !< x-position of agent |
---|
| 4487 | REAL(wp) :: xw !< x-position of wall |
---|
| 4488 | REAL(wp) :: ya !< x-position of agent |
---|
| 4489 | REAL(wp) :: yw !< y-position of wall |
---|
| 4490 | |
---|
| 4491 | force_d_x = 0.0_wp |
---|
| 4492 | force_d_y = 0.0_wp |
---|
| 4493 | ! |
---|
| 4494 | !-- calculate coordinates of corner relative to agent |
---|
| 4495 | !-- postion and distance between corner and agent |
---|
| 4496 | xw = xa - xw |
---|
| 4497 | yw = ya - yw |
---|
| 4498 | dist_l = SQRT( (xw)**2 + (yw)**2 ) |
---|
| 4499 | ! |
---|
| 4500 | !-- calculate x and y component of repulsive force |
---|
| 4501 | !-- induced by previously found corner |
---|
| 4502 | IF ( dist_l > 0 ) THEN |
---|
| 4503 | force_d_x = repuls_wall * d_sigma_rep_wall & |
---|
| 4504 | * EXP( -dist_l * d_sigma_rep_wall ) & |
---|
| 4505 | * xw / (dist_l) |
---|
| 4506 | force_d_y = repuls_wall * d_sigma_rep_wall & |
---|
| 4507 | * EXP( -dist_l * d_sigma_rep_wall ) & |
---|
| 4508 | * yw / (dist_l) |
---|
| 4509 | ENDIF |
---|
| 4510 | |
---|
| 4511 | ! !-- forces that are located outside of a sight radius of |
---|
| 4512 | ! !-- 200° (-> COS(100./180.*pi) = COS(.555*pi)) of |
---|
| 4513 | ! !-- current agent are considered to have an effect of 50% |
---|
| 4514 | ! IF ( force_d_x * agents(nl)%speed_e_x + & |
---|
| 4515 | ! force_d_y * agents(nl)%speed_e_y < & |
---|
| 4516 | ! SQRT(force_d_x**2 + force_d_y**2) * & |
---|
| 4517 | ! COS( .55555555 * 3.1415 ) ) & |
---|
| 4518 | ! THEN |
---|
| 4519 | ! force_d_x = force_d_x * .5_wp |
---|
| 4520 | ! force_d_y = force_d_y * .5_wp |
---|
| 4521 | ! ENDIF |
---|
| 4522 | |
---|
| 4523 | ! |
---|
| 4524 | !-- add force increment to total force of current agent |
---|
| 4525 | force_x = force_x + force_d_x |
---|
| 4526 | force_y = force_y + force_d_y |
---|
| 4527 | |
---|
| 4528 | END SUBROUTINE mas_timestep_wall_corner_force |
---|
| 4529 | |
---|
| 4530 | ! |
---|
| 4531 | !-- Calculates distance of point P to edge (A,B). If A = B, calculates |
---|
| 4532 | !-- point-to-point distance from A/B to P |
---|
| 4533 | FUNCTION dist_point_to_edge ( a_x, a_y, b_x, b_y, p_x, p_y ) |
---|
| 4534 | |
---|
| 4535 | IMPLICIT NONE |
---|
| 4536 | |
---|
| 4537 | REAL(wp) :: ab_x !< x-coordinate of vector from A to B |
---|
| 4538 | REAL(wp) :: ab_y !< y-coordinate of vector from A to B |
---|
| 4539 | REAL(wp) :: ab_d !< inverse length of vector from A to B |
---|
| 4540 | REAL(wp) :: ab_u_x !< x-coordinate of vector with direction of ab and length 1 |
---|
| 4541 | REAL(wp) :: ab_u_y !< y-coordinate of vector with direction of ab and length 1 |
---|
| 4542 | REAL(wp) :: ba_x !< x-coordinate of vector from B to A |
---|
| 4543 | REAL(wp) :: ba_y !< y-coordinate of vector from B to A |
---|
| 4544 | REAL(wp) :: ap_x !< x-coordinate of vector from A to P |
---|
| 4545 | REAL(wp) :: ap_y !< y-coordinate of vector from A to P |
---|
| 4546 | REAL(wp) :: bp_x !< x-coordinate of vector from B to P |
---|
| 4547 | REAL(wp) :: bp_y !< y-coordinate of vector from B to P |
---|
| 4548 | REAL(wp) :: a_x !< x-coordinate of point A of edge |
---|
| 4549 | REAL(wp) :: a_y !< y-coordinate of point A of edge |
---|
| 4550 | REAL(wp) :: b_x !< x-coordinate of point B of edge |
---|
| 4551 | REAL(wp) :: b_y !< y-coordinate of point B of edge |
---|
| 4552 | REAL(wp) :: p_x !< x-coordinate of point P |
---|
| 4553 | REAL(wp) :: p_y !< y-coordinate of point P |
---|
| 4554 | REAL(wp) :: dist_x !< x-coordinate of point P |
---|
| 4555 | REAL(wp) :: dist_y !< y-coordinate of point P |
---|
| 4556 | REAL(wp) :: dist_point_to_edge !< y-coordinate of point P |
---|
| 4557 | |
---|
| 4558 | ab_x = - a_x + b_x |
---|
| 4559 | ab_y = - a_y + b_y |
---|
| 4560 | ba_x = - b_x + a_x |
---|
| 4561 | ba_y = - b_y + a_y |
---|
| 4562 | ap_x = - a_x + p_x |
---|
| 4563 | ap_y = - a_y + p_y |
---|
| 4564 | bp_x = - b_x + p_x |
---|
| 4565 | bp_y = - b_y + p_y |
---|
| 4566 | |
---|
| 4567 | IF ( ab_x * ap_x + ab_y * ap_y <= 0. ) THEN |
---|
| 4568 | dist_point_to_edge = SQRT((a_x - p_x)**2 + (a_y - p_y)**2) |
---|
| 4569 | ELSEIF ( ba_x * bp_x + ba_y * bp_y <= 0. ) THEN |
---|
| 4570 | dist_point_to_edge = SQRT((b_x - p_x)**2 + (b_y - p_y)**2) |
---|
| 4571 | ELSE |
---|
| 4572 | ab_d = 1./SQRT((ab_x)**2+(ab_y)**2) |
---|
| 4573 | ab_u_x = ab_x*ab_d |
---|
| 4574 | ab_u_y = ab_y*ab_d |
---|
| 4575 | dist_x = ap_x - (ap_x*ab_u_x+ap_y*ab_u_y)*ab_u_x |
---|
| 4576 | dist_y = ap_y - (ap_x*ab_u_x+ap_y*ab_u_y)*ab_u_y |
---|
| 4577 | dist_point_to_edge = SQRT( dist_x**2 + dist_y**2 ) |
---|
| 4578 | ENDIF |
---|
| 4579 | |
---|
| 4580 | END FUNCTION dist_point_to_edge |
---|
| 4581 | |
---|
| 4582 | ! |
---|
| 4583 | !-- Returns the heuristic between points A and B (currently the straight |
---|
| 4584 | !-- distance) |
---|
| 4585 | FUNCTION heuristic ( ax, ay, bx, by ) |
---|
| 4586 | |
---|
| 4587 | IMPLICIT NONE |
---|
| 4588 | |
---|
| 4589 | REAL(wp) :: ax !< x-coordinate of point A |
---|
| 4590 | REAL(wp) :: ay !< y-coordinate of point A |
---|
| 4591 | REAL(wp) :: bx !< x-coordinate of point B |
---|
| 4592 | REAL(wp) :: by !< y-coordinate of point B |
---|
| 4593 | REAL(wp) :: heuristic !< return value |
---|
| 4594 | |
---|
| 4595 | heuristic = SQRT(( ax - bx )**2 + ( ay - by )**2) |
---|
| 4596 | |
---|
| 4597 | END FUNCTION heuristic |
---|
| 4598 | |
---|
| 4599 | ! |
---|
| 4600 | !-- Calculates if point P is left of the infinite |
---|
| 4601 | !-- line that contains A and B (direction: A to B) |
---|
| 4602 | !-- Concept: 2D rotation of two vectors |
---|
| 4603 | FUNCTION is_left ( ax, ay, bx, by, px, py ) |
---|
| 4604 | |
---|
| 4605 | IMPLICIT NONE |
---|
| 4606 | |
---|
| 4607 | LOGICAL :: is_left !< return value; TRUE if P is left of AB |
---|
| 4608 | |
---|
| 4609 | REAL(wp) :: ax !< x-coordinate of point A |
---|
| 4610 | REAL(wp) :: ay !< y-coordinate of point A |
---|
| 4611 | REAL(wp) :: bx !< x-coordinate of point B |
---|
| 4612 | REAL(wp) :: by !< y-coordinate of point B |
---|
| 4613 | REAL(wp) :: px !< x-coordinate of point P |
---|
| 4614 | REAL(wp) :: py !< y-coordinate of point P |
---|
| 4615 | |
---|
| 4616 | is_left = (bx-ax)*(py-ay)-(px-ax)*(by-ay) > 0 |
---|
| 4617 | IF ( (ABS(ax-px) < .001 .AND. ABS(ay-py) < .001) .OR. & |
---|
| 4618 | (ABS(bx-px) < .001 .AND. ABS(by-py) < .001) ) & |
---|
| 4619 | THEN |
---|
| 4620 | is_left = .FALSE. |
---|
| 4621 | ENDIF |
---|
| 4622 | |
---|
| 4623 | RETURN |
---|
| 4624 | |
---|
| 4625 | END FUNCTION is_left |
---|
| 4626 | |
---|
| 4627 | ! |
---|
| 4628 | !-- Calculates if point P is right of the infinite |
---|
| 4629 | !-- line that contains A and B (direction: A to B) |
---|
| 4630 | !-- Concept: 2D rotation of two vectors |
---|
| 4631 | FUNCTION is_right ( ax, ay, bx, by, px, py ) |
---|
| 4632 | |
---|
| 4633 | IMPLICIT NONE |
---|
| 4634 | |
---|
| 4635 | LOGICAL :: is_right !< return value; TRUE if P is right of AB |
---|
| 4636 | |
---|
| 4637 | REAL(wp), INTENT(IN) :: ax !< x-coordinate of point A |
---|
| 4638 | REAL(wp), INTENT(IN) :: ay !< y-coordinate of point A |
---|
| 4639 | REAL(wp), INTENT(IN) :: bx !< x-coordinate of point B |
---|
| 4640 | REAL(wp), INTENT(IN) :: by !< y-coordinate of point B |
---|
| 4641 | REAL(wp), INTENT(IN) :: px !< x-coordinate of point P |
---|
| 4642 | REAL(wp), INTENT(IN) :: py !< y-coordinate of point P |
---|
| 4643 | |
---|
| 4644 | is_right = (bx-ax)*(py-ay)-(px-ax)*(by-ay) < 0 |
---|
| 4645 | IF ( (ABS(ax-px) < .001 .AND. ABS(ay-py) < .001) .OR. & |
---|
| 4646 | (ABS(bx-px) < .001 .AND. ABS(by-py) < .001) ) & |
---|
| 4647 | THEN |
---|
| 4648 | is_right = .FALSE. |
---|
| 4649 | ENDIF |
---|
| 4650 | |
---|
| 4651 | RETURN |
---|
| 4652 | |
---|
| 4653 | END FUNCTION is_right |
---|
| 4654 | |
---|
| 4655 | ! |
---|
| 4656 | !-- Returns true if the line segments AB and PQ share an intersection |
---|
| 4657 | FUNCTION intersect ( ax, ay, bx, by, px, py, qx, qy ) |
---|
| 4658 | |
---|
| 4659 | IMPLICIT NONE |
---|
| 4660 | |
---|
| 4661 | LOGICAL :: intersect !< return value; TRUE if intersection was found |
---|
| 4662 | LOGICAL :: la !< T if a is left of PQ |
---|
| 4663 | LOGICAL :: lb !< T if b is left of PQ |
---|
| 4664 | LOGICAL :: lp !< T if p is left of AB |
---|
| 4665 | LOGICAL :: lq !< T if q is left of AB |
---|
| 4666 | LOGICAL :: poss !< flag that indicates if an intersection is still possible |
---|
| 4667 | LOGICAL :: ra !< T if a is right of PQ |
---|
| 4668 | LOGICAL :: rb !< T if b is right of PQ |
---|
| 4669 | LOGICAL :: rp !< T if p is right of AB |
---|
| 4670 | LOGICAL :: rq !< T if q is right of AB |
---|
| 4671 | |
---|
| 4672 | REAL(wp) :: ax !< x-coordinate of point A |
---|
| 4673 | REAL(wp) :: ay !< y-coordinate of point A |
---|
| 4674 | REAL(wp) :: bx !< x-coordinate of point B |
---|
| 4675 | REAL(wp) :: by !< y-coordinate of point B |
---|
| 4676 | REAL(wp) :: px !< x-coordinate of point P |
---|
| 4677 | REAL(wp) :: py !< y-coordinate of point P |
---|
| 4678 | REAL(wp) :: qx !< x-coordinate of point Q |
---|
| 4679 | REAL(wp) :: qy !< y-coordinate of point Q |
---|
| 4680 | |
---|
| 4681 | intersect = .FALSE. |
---|
| 4682 | poss = .FALSE. |
---|
| 4683 | ! |
---|
| 4684 | !-- Intersection is possible only if P and Q are on opposing sides of AB |
---|
| 4685 | lp = is_left(ax,ay,bx,by,px,py) |
---|
| 4686 | rq = is_right(ax,ay,bx,by,qx,qy) |
---|
| 4687 | IF ( lp .AND. rq ) poss = .TRUE. |
---|
| 4688 | IF ( .NOT. poss ) THEN |
---|
| 4689 | lq = is_left(ax,ay,bx,by,qx,qy) |
---|
| 4690 | rp = is_right(ax,ay,bx,by,px,py) |
---|
| 4691 | IF ( lq .AND. rp ) poss = .TRUE. |
---|
| 4692 | ENDIF |
---|
| 4693 | ! |
---|
| 4694 | !-- Intersection occurs only if above test (poss) was true AND |
---|
| 4695 | !-- A and B are on opposing sides of PQ |
---|
| 4696 | IF ( poss ) THEN |
---|
| 4697 | la = is_left(px,py,qx,qy,ax,ay) |
---|
| 4698 | rb = is_right(px,py,qx,qy,bx,by) |
---|
| 4699 | IF ( la .AND. rb ) intersect = .TRUE. |
---|
| 4700 | IF ( .NOT. intersect ) THEN |
---|
| 4701 | lb = is_left(px,py,qx,qy,bx,by) |
---|
| 4702 | ra = is_right(px,py,qx,qy,ax,ay) |
---|
| 4703 | IF ( lb .AND. ra ) intersect = .TRUE. |
---|
| 4704 | ENDIF |
---|
| 4705 | ENDIF |
---|
| 4706 | |
---|
| 4707 | RETURN |
---|
| 4708 | |
---|
| 4709 | END FUNCTION intersect |
---|
| 4710 | |
---|
| 4711 | ! |
---|
| 4712 | !-- Gives a nuber randomly distributed around an average |
---|
| 4713 | FUNCTION random_normal ( avg, variation ) |
---|
| 4714 | |
---|
| 4715 | IMPLICIT NONE |
---|
| 4716 | |
---|
| 4717 | REAL(wp) :: avg !< x-coordinate of vector from A to B |
---|
| 4718 | REAL(wp) :: variation !< y-coordinate of vector from A to B |
---|
| 4719 | REAL(wp) :: random_normal !< y-coordinate of vector from A to B |
---|
| 4720 | |
---|
| 4721 | REAL(wp), DIMENSION(12) :: random_arr !< inverse length of vector from A to B |
---|
| 4722 | |
---|
| 4723 | CALL RANDOM_NUMBER(random_arr) |
---|
| 4724 | random_normal = avg + variation*(SUM(random_arr)-6.) |
---|
| 4725 | |
---|
| 4726 | END FUNCTION random_normal |
---|
| 4727 | |
---|
| 4728 | |
---|
| 4729 | END MODULE multi_agent_system_mod |
---|