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