Home

Context Navigation

source: palm/trunk/SOURCE/multi_agent_system_mod.f90 @ 3746

Last change on this file since 3746 was 3665, checked in by raasch, 6 years ago
dummy statements added to avoid compiler warnings about unused variables, unused variables removed, ssh-call for submitting batch jobs on remote systems modified again to avoid output of login messages on specific systems
Property svn:keywords set to `Id`
File size: 198.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |