source: palm/trunk/SOURCE/multi_agent_system_mod.f90

Last change on this file was 4843, checked in by raasch, 3 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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