source: palm/trunk/SOURCE/lpm_exchange_horiz.f90 @ 1359

Last change on this file since 1359 was 1359, checked in by hoffmann, 8 years ago

new Lagrangian particle structure integrated

  • Property svn:keywords set to Id
File size: 53.3 KB
Line 
1 MODULE lpm_exchange_horiz_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! New particle structure integrated.
23! Kind definition added to all floating point numbers.
24!
25! Former revisions:
26! -----------------
27! $Id: lpm_exchange_horiz.f90 1359 2014-04-11 17:15:14Z hoffmann $
28!
29! 1327 2014-03-21 11:00:16Z raasch
30! -netcdf output queries
31!
32! 1320 2014-03-20 08:40:49Z raasch
33! ONLY-attribute added to USE-statements,
34! kind-parameters added to all INTEGER and REAL declaration statements,
35! kinds are defined in new module kinds,
36! comment fields (!:) to be used for variable explanations added to
37! all variable declaration statements
38!
39! 1318 2014-03-17 13:35:16Z raasch
40! module interfaces removed
41!
42! 1036 2012-10-22 13:43:42Z raasch
43! code put under GPL (PALM 3.9)
44!
45! 851 2012-03-15 14:32:58Z raasch
46! Bugfix: resetting of particle_mask and tail mask moved from end of this
47! routine to lpm
48!
49! 849 2012-03-15 10:35:09Z raasch
50! initial revision (former part of advec_particles)
51!
52!
53! Description:
54! ------------
55! Exchange of particles (and tails) between the subdomains.
56!------------------------------------------------------------------------------!
57
58    USE control_parameters,                                                    &
59        ONLY:  dz, message_string, netcdf_data_format, simulated_time
60
61    USE cpulog,                                                                &
62        ONLY:  cpu_log, log_point_s
63
64    USE grid_variables,                                                        &
65        ONLY:  ddx, ddy, dx, dy
66
67    USE indices,                                                               &
68        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzt
69
70    USE kinds
71
72    USE lpm_pack_arrays_mod,                                                   &
73        ONLY:  lpm_pack_arrays
74
75    USE particle_attributes,                                                   &
76        ONLY:  alloc_factor, deleted_particles, deleted_tails, grid_particles, &
77               ibc_par_lr, ibc_par_ns, maximum_number_of_tails,                &
78               maximum_number_of_tailpoints, min_nr_particle,                  &
79               mpi_particle_type, number_of_tails, number_of_particles,        &
80               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
81               particle_tail_coordinates, particle_type, prt_count,            &
82               tail_mask, trlp_count_sum,                                      &
83               trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
84               trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
85               trsp_count_recv_sum, use_particle_tails, zero_particle
86
87    USE pegrid
88
89    IMPLICIT NONE
90
91    INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !:
92    INTEGER(iwp)            ::  nr_move_north               !:
93    INTEGER(iwp)            ::  nr_move_south               !:
94
95    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_north
96    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_south
97
98    SAVE
99
100    PRIVATE
101    PUBLIC lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
102
103    INTERFACE lpm_exchange_horiz
104       MODULE PROCEDURE lpm_exchange_horiz
105    END INTERFACE lpm_exchange_horiz
106
107    INTERFACE lpm_move_particle
108       MODULE PROCEDURE lpm_move_particle
109    END INTERFACE lpm_move_particle
110
111    INTERFACE realloc_particles_array
112       MODULE PROCEDURE realloc_particles_array
113    END INTERFACE realloc_particles_array
114
115CONTAINS
116
117 SUBROUTINE lpm_exchange_horiz
118
119    IMPLICIT NONE
120
121    INTEGER(iwp) ::  i                                       !:
122    INTEGER(iwp) ::  ip                                      !:
123    INTEGER(iwp) ::  j                                       !:
124    INTEGER(iwp) ::  jp                                      !:
125    INTEGER(iwp) ::  k                                       !:
126    INTEGER(iwp) ::  kp                                      !:
127    INTEGER(iwp) ::  n                                       !:
128    INTEGER(iwp) ::  nn                                      !:
129    INTEGER(iwp) ::  tlength                                 !:
130    INTEGER(iwp) ::  trlp_count                              !:
131    INTEGER(iwp) ::  trlp_count_recv                         !:
132    INTEGER(iwp) ::  trlpt_count                             !:
133    INTEGER(iwp) ::  trlpt_count_recv                        !:
134    INTEGER(iwp) ::  trnp_count                              !:
135    INTEGER(iwp) ::  trnp_count_recv                         !:
136    INTEGER(iwp) ::  trnpt_count                             !:
137    INTEGER(iwp) ::  trnpt_count_recv                        !:
138    INTEGER(iwp) ::  trrp_count                              !:
139    INTEGER(iwp) ::  trrp_count_recv                         !:
140    INTEGER(iwp) ::  trrpt_count                             !:
141    INTEGER(iwp) ::  trrpt_count_recv                        !:
142    INTEGER(iwp) ::  trsp_count                              !:
143    INTEGER(iwp) ::  trsp_count_recv                         !:
144    INTEGER(iwp) ::  trspt_count                             !:
145    INTEGER(iwp) ::  trspt_count_recv                        !:
146
147    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trlpt        !:
148    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trnpt        !:
149    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trrpt        !:
150    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trspt        !:
151
152    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !:
153    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !:
154    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !:
155    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !:
156    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !:
157    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !:
158    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !:
159    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !:
160
161    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
162
163#if defined( __parallel )
164
165!
166!-- Exchange between subdomains.
167!-- As soon as one particle has moved beyond the boundary of the domain, it
168!-- is included in the relevant transfer arrays and marked for subsequent
169!-- deletion on this PE.
170!-- First sweep for crossings in x direction. Find out first the number of
171!-- particles to be transferred and allocate temporary arrays needed to store
172!-- them.
173!-- For a one-dimensional decomposition along y, no transfer is necessary,
174!-- because the particle remains on the PE, but the particle coordinate has to
175!-- be adjusted.
176    trlp_count  = 0
177    trlpt_count = 0
178    trrp_count  = 0
179    trrpt_count = 0
180
181    trlp_count_recv   = 0
182    trlpt_count_recv  = 0
183    trrp_count_recv   = 0
184    trrpt_count_recv  = 0
185
186    IF ( pdims(1) /= 1 )  THEN
187!
188!--    First calculate the storage necessary for sending and receiving the data.
189!--    Compute only first (nxl) and last (nxr) loop iterration.
190       DO  ip = nxl, nxr, nxr - nxl
191          DO  jp = nys, nyn
192             DO  kp = nzb+1, nzt
193
194                number_of_particles = prt_count(kp,jp,ip)
195                IF ( number_of_particles <= 0 )  CYCLE
196                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
197                DO  n = 1, number_of_particles
198                   IF ( particles(n)%particle_mask )  THEN
199                      i = ( particles(n)%x + 0.5_wp * dx ) * ddx
200   !
201   !--                Above calculation does not work for indices less than zero
202                      IF ( particles(n)%x < -0.5_wp * dx )  i = -1
203
204                      IF ( i < nxl )  THEN
205                         trlp_count = trlp_count + 1
206                         IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
207                      ELSEIF ( i > nxr )  THEN
208                         trrp_count = trrp_count + 1
209                         IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
210                      ENDIF
211                   ENDIF
212                ENDDO
213
214             ENDDO
215          ENDDO
216       ENDDO
217
218       IF ( trlp_count  == 0 )  trlp_count  = 1
219       IF ( trlpt_count == 0 )  trlpt_count = 1
220       IF ( trrp_count  == 0 )  trrp_count  = 1
221       IF ( trrpt_count == 0 )  trrpt_count = 1
222
223       ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
224
225       trlp = zero_particle
226       trrp = zero_particle
227
228       IF ( use_particle_tails )  THEN
229          ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
230                    trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
231          tlength = maximum_number_of_tailpoints * 5
232       ENDIF
233
234       trlp_count  = 0
235       trlpt_count = 0
236       trrp_count  = 0
237       trrpt_count = 0
238
239    ENDIF
240!
241!-- Compute only first (nxl) and last (nxr) loop iterration
242    DO  ip = nxl, nxr,nxr-nxl
243       DO  jp = nys, nyn
244          DO  kp = nzb+1, nzt
245             number_of_particles = prt_count(kp,jp,ip)
246             IF ( number_of_particles <= 0 ) CYCLE
247             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
248             DO  n = 1, number_of_particles
249
250                nn = particles(n)%tail_id
251!
252!--             Only those particles that have not been marked as 'deleted' may
253!--             be moved.
254                IF ( particles(n)%particle_mask )  THEN
255
256                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
257   !
258   !--             Above calculation does not work for indices less than zero
259                   IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
260
261                   IF ( i <  nxl )  THEN
262                      IF ( i < 0 )  THEN
263   !
264   !--                   Apply boundary condition along x
265                         IF ( ibc_par_lr == 0 )  THEN
266   !
267   !--                      Cyclic condition
268                            IF ( pdims(1) == 1 )  THEN
269                               particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
270                               particles(n)%origin_x = ( nx + 1 ) * dx + &
271                               particles(n)%origin_x
272                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
273                                  = particles(n)%tailpoints
274                                  particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
275                                  + particle_tail_coordinates(1:i,1,nn)
276                               ENDIF
277                            ELSE
278                               trlp_count = trlp_count + 1
279                               trlp(trlp_count)   = particles(n)
280                               trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
281                               trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
282                               ( nx + 1 ) * dx
283                               particles(n)%particle_mask  = .FALSE.
284                               deleted_particles = deleted_particles + 1
285
286                               IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp )  THEN
287                                  trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
288                                  !++ why is 1 subtracted in next statement???
289                                  trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
290                               ENDIF
291
292                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
293                                  trlpt_count = trlpt_count + 1
294                                  trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
295                                  trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
296                                  trlpt(:,1,trlpt_count)
297                                  tail_mask(nn) = .FALSE.
298                                  deleted_tails = deleted_tails + 1
299                               ENDIF
300                            ENDIF
301
302                         ELSEIF ( ibc_par_lr == 1 )  THEN
303   !
304   !--                      Particle absorption
305                            particles(n)%particle_mask = .FALSE.
306                            deleted_particles = deleted_particles + 1
307                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
308                               tail_mask(nn) = .FALSE.
309                               deleted_tails = deleted_tails + 1
310                            ENDIF
311
312                         ELSEIF ( ibc_par_lr == 2 )  THEN
313   !
314   !--                      Particle reflection
315                            particles(n)%x       = -particles(n)%x
316                            particles(n)%speed_x = -particles(n)%speed_x
317
318                         ENDIF
319                      ELSE
320   !
321   !--                   Store particle data in the transfer array, which will be
322   !--                   send to the neighbouring PE
323                         trlp_count = trlp_count + 1
324                         trlp(trlp_count) = particles(n)
325                         particles(n)%particle_mask = .FALSE.
326                         deleted_particles = deleted_particles + 1
327
328                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
329                            trlpt_count = trlpt_count + 1
330                            trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
331                            tail_mask(nn) = .FALSE.
332                            deleted_tails = deleted_tails + 1
333                         ENDIF
334                      ENDIF
335
336                   ELSEIF ( i > nxr )  THEN
337                      IF ( i > nx )  THEN
338   !
339   !--                   Apply boundary condition along x
340                         IF ( ibc_par_lr == 0 )  THEN
341   !
342   !--                      Cyclic condition
343                            IF ( pdims(1) == 1 )  THEN
344                               particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
345                               particles(n)%origin_x = particles(n)%origin_x - &
346                               ( nx + 1 ) * dx
347                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
348                                  i = particles(n)%tailpoints
349                                  particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
350                                  + particle_tail_coordinates(1:i,1,nn)
351                               ENDIF
352                            ELSE
353                               trrp_count = trrp_count + 1
354                               trrp(trrp_count) = particles(n)
355                               trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
356                               trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
357                               ( nx + 1 ) * dx
358                               particles(n)%particle_mask = .FALSE.
359                               deleted_particles = deleted_particles + 1
360
361                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
362                                  trrpt_count = trrpt_count + 1
363                                  trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
364                                  trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
365                                  ( nx + 1 ) * dx
366                                  tail_mask(nn) = .FALSE.
367                                  deleted_tails = deleted_tails + 1
368                               ENDIF
369                            ENDIF
370
371                         ELSEIF ( ibc_par_lr == 1 )  THEN
372   !
373   !--                      Particle absorption
374                            particles(n)%particle_mask = .FALSE.
375                            deleted_particles = deleted_particles + 1
376                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
377                               tail_mask(nn) = .FALSE.
378                               deleted_tails = deleted_tails + 1
379                            ENDIF
380
381                         ELSEIF ( ibc_par_lr == 2 )  THEN
382   !
383   !--                      Particle reflection
384                            particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
385                            particles(n)%speed_x = -particles(n)%speed_x
386
387                         ENDIF
388                      ELSE
389   !
390   !--                   Store particle data in the transfer array, which will be send
391   !--                   to the neighbouring PE
392                         trrp_count = trrp_count + 1
393                         trrp(trrp_count) = particles(n)
394                         particles(n)%particle_mask = .FALSE.
395                         deleted_particles = deleted_particles + 1
396
397                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
398                            trrpt_count = trrpt_count + 1
399                            trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
400                            tail_mask(nn) = .FALSE.
401                            deleted_tails = deleted_tails + 1
402                         ENDIF
403                      ENDIF
404
405                   ENDIF
406                ENDIF
407             ENDDO
408          ENDDO
409       ENDDO
410    ENDDO
411
412!
413!-- Send left boundary, receive right boundary (but first exchange how many
414!-- and check, if particle storage must be extended)
415    IF ( pdims(1) /= 1 )  THEN
416
417       CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
418                          trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
419                          comm2d, status, ierr )
420
421       ALLOCATE(rvrp(MAX(1,trrp_count_recv)))
422
423       CALL MPI_SENDRECV( trlp(1)%radius, max(1,trlp_count), mpi_particle_type,&
424                          pleft, 1, rvrp(1)%radius,                            &
425                          max(1,trrp_count_recv), mpi_particle_type, pright, 1,&
426                          comm2d, status, ierr )
427
428       IF ( trrp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvrp(1:trrp_count_recv))
429
430       DEALLOCATE(rvrp)
431
432       IF ( use_particle_tails )  THEN
433
434          CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
435                             trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
436                             comm2d, status, ierr )
437
438          IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
439          THEN
440             IF ( netcdf_data_format < 3 )  THEN
441                message_string = 'maximum_number_of_tails ' //   &
442                                 'needs to be increased ' //     &
443                                 '&but this is not allowed wi'// &
444                                 'th netcdf_data_format < 3'
445                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
446             ELSE
447                CALL lpm_extend_tail_array( trrpt_count_recv )
448             ENDIF
449          ENDIF
450
451          CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,      &
452                             pleft, 1,                                         &
453                             particle_tail_coordinates(1,1,number_of_tails+1), &
454                             trrpt_count_recv*tlength, MPI_REAL, pright, 1,    &
455                             comm2d, status, ierr )
456!
457!--       Update the tail ids for the transferred particles
458          nn = number_of_tails
459          DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
460             IF ( particles(n)%tail_id /= 0 )  THEN
461                nn = nn + 1
462                particles(n)%tail_id = nn
463             ENDIF
464          ENDDO
465
466       ENDIF
467
468!
469!--    Send right boundary, receive left boundary
470       CALL MPI_SENDRECV( trrp_count,      1, MPI_INTEGER, pright, 0, &
471                          trlp_count_recv, 1, MPI_INTEGER, pleft,  0, &
472                          comm2d, status, ierr )
473
474       ALLOCATE(rvlp(MAX(1,trlp_count_recv)))
475
476       CALL MPI_SENDRECV( trrp(1)%radius, max(1,trrp_count), mpi_particle_type,&
477                          pright, 1, rvlp(1)%radius,                           &
478                          max(1,trlp_count_recv), mpi_particle_type, pleft, 1, &
479                          comm2d, status, ierr )
480
481       IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
482
483       DEALLOCATE(rvlp)
484
485       IF ( use_particle_tails )  THEN
486
487          CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
488                             trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
489                             comm2d, status, ierr )
490
491          IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
492          THEN
493             IF ( netcdf_data_format < 3 )  THEN
494                message_string = 'maximum_number_of_tails ' //   &
495                                 'needs to be increased ' //     &
496                                 '&but this is not allowed wi'// &
497                                 'th netcdf_data_format < 3'
498                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
499             ELSE
500                CALL lpm_extend_tail_array( trlpt_count_recv )
501             ENDIF
502          ENDIF
503
504          CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,      &
505                             pright, 1,                                        &
506                             particle_tail_coordinates(1,1,number_of_tails+1), &
507                             trlpt_count_recv*tlength, MPI_REAL, pleft, 1,     &
508                             comm2d, status, ierr )
509!
510!--       Update the tail ids for the transferred particles
511          nn = number_of_tails
512          DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
513             IF ( particles(n)%tail_id /= 0 )  THEN
514                nn = nn + 1
515                particles(n)%tail_id = nn
516             ENDIF
517          ENDDO
518
519       ENDIF
520
521!       number_of_particles = number_of_particles + trlp_count_recv
522!       number_of_tails     = number_of_tails     + trlpt_count_recv
523
524       IF ( use_particle_tails )  THEN
525          DEALLOCATE( trlpt, trrpt )
526       ENDIF
527       DEALLOCATE( trlp, trrp )
528
529    ENDIF
530
531
532!
533!-- Check whether particles have crossed the boundaries in y direction. Note
534!-- that this case can also apply to particles that have just been received
535!-- from the adjacent right or left PE.
536!-- Find out first the number of particles to be transferred and allocate
537!-- temporary arrays needed to store them.
538!-- For a one-dimensional decomposition along x, no transfer is necessary,
539!-- because the particle remains on the PE.
540    trsp_count  = nr_move_south
541    trspt_count = 0
542    trnp_count  = nr_move_north
543    trnpt_count = 0
544
545    trsp_count_recv   = 0
546    trspt_count_recv  = 0
547    trnp_count_recv   = 0
548    trnpt_count_recv  = 0
549
550    IF ( pdims(2) /= 1 )  THEN
551!
552!--    First calculate the storage necessary for sending and receiving the
553!--    data
554       DO  ip = nxl, nxr
555          DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
556             DO  kp = nzb+1, nzt
557                number_of_particles = prt_count(kp,jp,ip)
558                IF ( number_of_particles <= 0 )  CYCLE
559                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
560                DO  n = 1, number_of_particles
561                   IF ( particles(n)%particle_mask )  THEN
562                      j = ( particles(n)%y + 0.5_wp * dy ) * ddy
563!
564!--                   Above calculation does not work for indices less than zero
565                      IF ( particles(n)%y < -0.5_wp * dy )  j = -1
566
567                      IF ( j < nys )  THEN
568                         trsp_count = trsp_count + 1
569                         IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count + 1
570                      ELSEIF ( j > nyn )  THEN
571                         trnp_count = trnp_count + 1
572                         IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count + 1
573                      ENDIF
574                   ENDIF
575                ENDDO
576             ENDDO
577          ENDDO
578       ENDDO
579
580       IF ( trsp_count  == 0 )  trsp_count  = 1
581       IF ( trspt_count == 0 )  trspt_count = 1
582       IF ( trnp_count  == 0 )  trnp_count  = 1
583       IF ( trnpt_count == 0 )  trnpt_count = 1
584
585       ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
586
587       trsp = zero_particle
588       trnp = zero_particle
589
590       IF ( use_particle_tails )  THEN
591          ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
592                    trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
593          tlength = maximum_number_of_tailpoints * 5
594       ENDIF
595
596       trsp_count  = nr_move_south
597       trspt_count = 0
598       trnp_count  = nr_move_north
599       trnpt_count = 0
600
601       trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
602       trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
603
604    ENDIF
605
606    DO  ip = nxl, nxr
607       DO  jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration
608          DO  kp = nzb+1, nzt
609             number_of_particles = prt_count(kp,jp,ip)
610             IF ( number_of_particles <= 0 )  CYCLE
611             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
612             DO  n = 1, number_of_particles
613
614                nn = particles(n)%tail_id
615!
616!--             Only those particles that have not been marked as 'deleted' may
617!--             be moved.
618                IF ( particles(n)%particle_mask )  THEN
619                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
620!
621!--                Above calculation does not work for indices less than zero
622                   IF ( particles(n)%y < -0.5_wp * dy )  j = -1
623
624                   IF ( j < nys )  THEN
625                      IF ( j < 0 )  THEN
626!
627!--                      Apply boundary condition along y
628                         IF ( ibc_par_ns == 0 )  THEN
629!
630!--                         Cyclic condition
631                            IF ( pdims(2) == 1 )  THEN
632                               particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
633                               particles(n)%origin_y = ( ny + 1 ) * dy + &
634                               particles(n)%origin_y
635                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
636                                  i = particles(n)%tailpoints
637                                  particle_tail_coordinates(1:i,2,nn) =        &
638                                     ( ny+1 ) * dy + particle_tail_coordinates(1:i,2,nn)
639                               ENDIF
640                            ELSE
641                               trsp_count = trsp_count + 1
642                               trsp(trsp_count) = particles(n)
643                               trsp(trsp_count)%y = ( ny + 1 ) * dy + &
644                               trsp(trsp_count)%y
645                               trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
646                               + ( ny + 1 ) * dy
647                               particles(n)%particle_mask = .FALSE.
648                               deleted_particles = deleted_particles + 1
649
650                               IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp )  THEN
651                                  trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
652                                  !++ why is 1 subtracted in next statement???
653                                  trsp(trsp_count)%origin_y =                        &
654                                  trsp(trsp_count)%origin_y - 1
655                               ENDIF
656
657                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
658                                  trspt_count = trspt_count + 1
659                                  trspt(:,:,trspt_count) = &
660                                  particle_tail_coordinates(:,:,nn)
661                                  trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
662                                  trspt(:,2,trspt_count)
663                                  tail_mask(nn) = .FALSE.
664                                  deleted_tails = deleted_tails + 1
665                               ENDIF
666                            ENDIF
667
668                         ELSEIF ( ibc_par_ns == 1 )  THEN
669!
670!--                         Particle absorption
671                            particles(n)%particle_mask = .FALSE.
672                            deleted_particles = deleted_particles + 1
673                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
674                               tail_mask(nn) = .FALSE.
675                               deleted_tails = deleted_tails + 1
676                            ENDIF
677
678                         ELSEIF ( ibc_par_ns == 2 )  THEN
679!
680!--                         Particle reflection
681                            particles(n)%y       = -particles(n)%y
682                            particles(n)%speed_y = -particles(n)%speed_y
683
684                         ENDIF
685                      ELSE
686!
687!--                      Store particle data in the transfer array, which will
688!--                      be send to the neighbouring PE
689                         trsp_count = trsp_count + 1
690                         trsp(trsp_count) = particles(n)
691                         particles(n)%particle_mask = .FALSE.
692                         deleted_particles = deleted_particles + 1
693
694                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
695                            trspt_count = trspt_count + 1
696                            trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
697                            tail_mask(nn) = .FALSE.
698                            deleted_tails = deleted_tails + 1
699                         ENDIF
700                      ENDIF
701
702                   ELSEIF ( j > nyn )  THEN
703                      IF ( j > ny )  THEN
704!
705!--                       Apply boundary condition along x
706                         IF ( ibc_par_ns == 0 )  THEN
707!
708!--                         Cyclic condition
709                            IF ( pdims(2) == 1 )  THEN
710                               particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
711                               particles(n)%origin_y =                         &
712                                  particles(n)%origin_y - ( ny + 1 ) * dy
713                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
714                                  i = particles(n)%tailpoints
715                                  particle_tail_coordinates(1:i,2,nn) =        &
716                                     - (ny+1) * dy + particle_tail_coordinates(1:i,2,nn)
717                               ENDIF
718                            ELSE
719                               trnp_count = trnp_count + 1
720                               trnp(trnp_count) = particles(n)
721                               trnp(trnp_count)%y =                            &
722                                  trnp(trnp_count)%y - ( ny + 1 ) * dy
723                               trnp(trnp_count)%origin_y =                     &
724                                  trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
725                               particles(n)%particle_mask = .FALSE.
726                               deleted_particles = deleted_particles + 1
727                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
728                                  trnpt_count = trnpt_count + 1
729                                  trnpt(:,:,trnpt_count) =                     &
730                                     particle_tail_coordinates(:,:,nn)
731                                  trnpt(:,2,trnpt_count) =                     &
732                                     trnpt(:,2,trnpt_count) - ( ny + 1 ) * dy
733                                  tail_mask(nn) = .FALSE.
734                                  deleted_tails = deleted_tails + 1
735                               ENDIF
736                            ENDIF
737
738                         ELSEIF ( ibc_par_ns == 1 )  THEN
739!
740!--                         Particle absorption
741                            particles(n)%particle_mask = .FALSE.
742                            deleted_particles = deleted_particles + 1
743                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
744                               tail_mask(nn) = .FALSE.
745                               deleted_tails = deleted_tails + 1
746                            ENDIF
747
748                         ELSEIF ( ibc_par_ns == 2 )  THEN
749!
750!--                         Particle reflection
751                            particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
752                            particles(n)%speed_y = -particles(n)%speed_y
753
754                         ENDIF
755                      ELSE
756!
757!--                      Store particle data in the transfer array, which will
758!--                      be send to the neighbouring PE
759                         trnp_count = trnp_count + 1
760                         trnp(trnp_count) = particles(n)
761                         particles(n)%particle_mask = .FALSE.
762                         deleted_particles = deleted_particles + 1
763
764                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
765                            trnpt_count = trnpt_count + 1
766                            trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
767                            tail_mask(nn) = .FALSE.
768                            deleted_tails = deleted_tails + 1
769                         ENDIF
770                      ENDIF
771
772                   ENDIF
773                ENDIF
774             ENDDO
775          ENDDO
776       ENDDO
777    ENDDO
778
779!
780!-- Send front boundary, receive back boundary (but first exchange how many
781!-- and check, if particle storage must be extended)
782    IF ( pdims(2) /= 1 )  THEN
783
784       CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
785                          trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
786                          comm2d, status, ierr )
787
788       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
789
790       CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
791                          psouth, 1, rvnp(1)%radius,                             &
792                          trnp_count_recv, mpi_particle_type, pnorth, 1,   &
793                          comm2d, status, ierr )
794
795       IF ( trnp_count_recv  > 0 )  CALL Add_particles_to_gridcell(rvnp(1:trnp_count_recv))
796
797       DEALLOCATE(rvnp)
798
799       IF ( use_particle_tails )  THEN
800
801          CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
802                             trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
803                             comm2d, status, ierr )
804
805          IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
806          THEN
807             IF ( netcdf_data_format < 3 )  THEN
808                message_string = 'maximum_number_of_tails ' //    &
809                                 'needs to be increased ' //      &
810                                 '&but this is not allowed wi' // &
811                                 'th netcdf_data_format < 3'
812                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
813             ELSE
814                CALL lpm_extend_tail_array( trnpt_count_recv )
815             ENDIF
816          ENDIF
817
818          CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,      &
819                             psouth, 1,                                        &
820                             particle_tail_coordinates(1,1,number_of_tails+1), &
821                             trnpt_count_recv*tlength, MPI_REAL, pnorth, 1,    &
822                             comm2d, status, ierr )
823
824!
825!--       Update the tail ids for the transferred particles
826          nn = number_of_tails
827          DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
828             IF ( particles(n)%tail_id /= 0 )  THEN
829                nn = nn + 1
830                particles(n)%tail_id = nn
831             ENDIF
832          ENDDO
833
834       ENDIF
835
836!       number_of_particles = number_of_particles + trnp_count_recv
837!       number_of_tails     = number_of_tails     + trnpt_count_recv
838
839!
840!--    Send back boundary, receive front boundary
841       CALL MPI_SENDRECV( trnp_count,      1, MPI_INTEGER, pnorth, 0, &
842                          trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
843                          comm2d, status, ierr )
844
845       ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
846
847       CALL MPI_SENDRECV( trnp(1)%radius, trnp_count, mpi_particle_type,      &
848                          pnorth, 1, rvsp(1)%radius,                          &
849                          trsp_count_recv, mpi_particle_type, psouth, 1,   &
850                          comm2d, status, ierr )
851
852       IF ( trsp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvsp(1:trsp_count_recv))
853
854       DEALLOCATE(rvsp)
855
856       IF ( use_particle_tails )  THEN
857
858          CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
859                             trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
860                             comm2d, status, ierr )
861
862          IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
863          THEN
864             IF ( netcdf_data_format < 3 )  THEN
865                message_string = 'maximum_number_of_tails ' //   &
866                                 'needs to be increased ' //     &
867                                 '&but this is not allowed wi'// &
868                                 'th NetCDF output switched on'
869                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
870             ELSE
871                CALL lpm_extend_tail_array( trspt_count_recv )
872             ENDIF
873          ENDIF
874
875          CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,      &
876                             pnorth, 1,                                        &
877                             particle_tail_coordinates(1,1,number_of_tails+1), &
878                             trspt_count_recv*tlength, MPI_REAL, psouth, 1,    &
879                             comm2d, status, ierr )
880!
881!--       Update the tail ids for the transferred particles
882          nn = number_of_tails
883          DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
884             IF ( particles(n)%tail_id /= 0 )  THEN
885                nn = nn + 1
886                particles(n)%tail_id = nn
887             ENDIF
888          ENDDO
889
890       ENDIF
891
892       number_of_particles = number_of_particles + trsp_count_recv
893       number_of_tails     = number_of_tails     + trspt_count_recv
894
895       IF ( use_particle_tails )  THEN
896          DEALLOCATE( trspt, trnpt )
897       ENDIF
898       DEALLOCATE( trsp, trnp )
899
900    ENDIF
901
902#else
903
904!
905!-- Apply boundary conditions
906    DO  n = 1, number_of_particles
907
908       nn = particles(n)%tail_id
909
910       IF ( particles(n)%x < -0.5_wp * dx )  THEN
911
912          IF ( ibc_par_lr == 0 )  THEN
913!
914!--          Cyclic boundary. Relevant coordinate has to be changed.
915             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
916             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
917                i = particles(n)%tailpoints
918                particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
919                                             particle_tail_coordinates(1:i,1,nn)
920             ENDIF
921          ELSEIF ( ibc_par_lr == 1 )  THEN
922!
923!--          Particle absorption
924             particles(n)%particle_mask = .FALSE.
925             deleted_particles = deleted_particles + 1
926             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
927                tail_mask(nn) = .FALSE.
928                deleted_tails = deleted_tails + 1
929             ENDIF
930          ELSEIF ( ibc_par_lr == 2 )  THEN
931!
932!--          Particle reflection
933             particles(n)%x       = -dx - particles(n)%x
934             particles(n)%speed_x = -particles(n)%speed_x
935          ENDIF
936
937       ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
938
939          IF ( ibc_par_lr == 0 )  THEN
940!
941!--          Cyclic boundary. Relevant coordinate has to be changed.
942             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
943             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
944                i = particles(n)%tailpoints
945                particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
946                                             particle_tail_coordinates(1:i,1,nn)
947             ENDIF
948          ELSEIF ( ibc_par_lr == 1 )  THEN
949!
950!--          Particle absorption
951             particles(n)%particle_mask = .FALSE.
952             deleted_particles = deleted_particles + 1
953             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
954                tail_mask(nn) = .FALSE.
955                deleted_tails = deleted_tails + 1
956             ENDIF
957          ELSEIF ( ibc_par_lr == 2 )  THEN
958!
959!--          Particle reflection
960             particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
961             particles(n)%speed_x = -particles(n)%speed_x
962          ENDIF
963
964       ENDIF
965
966       IF ( particles(n)%y < -0.5_wp * dy )  THEN
967
968          IF ( ibc_par_ns == 0 )  THEN
969!
970!--          Cyclic boundary. Relevant coordinate has to be changed.
971             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
972             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
973                i = particles(n)%tailpoints
974                particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
975                                             particle_tail_coordinates(1:i,2,nn)
976             ENDIF
977          ELSEIF ( ibc_par_ns == 1 )  THEN
978!
979!--          Particle absorption
980             particles(n)%particle_mask = .FALSE.
981             deleted_particles = deleted_particles + 1
982             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
983                tail_mask(nn) = .FALSE.
984                deleted_tails = deleted_tails + 1
985             ENDIF
986          ELSEIF ( ibc_par_ns == 2 )  THEN
987!
988!--          Particle reflection
989             particles(n)%y       = -dy - particles(n)%y
990             particles(n)%speed_y = -particles(n)%speed_y
991          ENDIF
992
993       ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
994
995          IF ( ibc_par_ns == 0 )  THEN
996!
997!--          Cyclic boundary. Relevant coordinate has to be changed.
998             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
999             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
1000                i = particles(n)%tailpoints
1001                particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
1002                                             particle_tail_coordinates(1:i,2,nn)
1003             ENDIF
1004          ELSEIF ( ibc_par_ns == 1 )  THEN
1005!
1006!--          Particle absorption
1007             particles(n)%particle_mask = .FALSE.
1008             deleted_particles = deleted_particles + 1
1009             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
1010                tail_mask(nn) = .FALSE.
1011                deleted_tails = deleted_tails + 1
1012             ENDIF
1013          ELSEIF ( ibc_par_ns == 2 )  THEN
1014!
1015!--          Particle reflection
1016             particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
1017             particles(n)%speed_y = -particles(n)%speed_y
1018          ENDIF
1019
1020       ENDIF
1021    ENDDO
1022
1023#endif
1024
1025!
1026!-- Accumulate the number of particles transferred between the subdomains
1027#if defined( __parallel )
1028    trlp_count_sum      = trlp_count_sum      + trlp_count
1029    trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
1030    trrp_count_sum      = trrp_count_sum      + trrp_count
1031    trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
1032    trsp_count_sum      = trsp_count_sum      + trsp_count
1033    trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
1034    trnp_count_sum      = trnp_count_sum      + trnp_count
1035    trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
1036#endif
1037
1038    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'stop' )
1039
1040 END SUBROUTINE lpm_exchange_horiz
1041
1042 SUBROUTINE Add_particles_to_gridcell (particle_array)
1043
1044!
1045!-- If a particle moves from one processor to another, this subroutine moves
1046!-- the corresponding elements from the particle arrays of the old grid cells
1047!-- to the particle arrays of the new grid cells.
1048
1049    IMPLICIT NONE
1050
1051    INTEGER(iwp)        ::  ip        !:
1052    INTEGER(iwp)        ::  jp        !:
1053    INTEGER(iwp)        ::  kp        !:
1054    INTEGER(iwp)        ::  n         !:
1055    INTEGER(iwp)        ::  pindex    !:
1056
1057    LOGICAL             ::  pack_done !:
1058
1059    TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
1060
1061    pack_done     = .FALSE.
1062
1063    nr_move_north = 0
1064    nr_move_south = 0
1065
1066    DO n = 1, SIZE(particle_array)
1067       ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
1068       jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
1069       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt_m1
1070
1071       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
1072            .AND.  kp >= nzb+1  .AND.  kp <= nzt)  THEN ! particle stays on processor
1073          number_of_particles = prt_count(kp,jp,ip)
1074          particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1075
1076          pindex = prt_count(kp,jp,ip)+1
1077          IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
1078             IF ( pack_done )  THEN
1079                CALL realloc_particles_array (ip,jp,kp)
1080             ELSE
1081                CALL lpm_pack_arrays
1082                prt_count(kp,jp,ip) = number_of_particles
1083                pindex = prt_count(kp,jp,ip)+1
1084                IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
1085                   CALL realloc_particles_array (ip,jp,kp)
1086                ENDIF
1087                pack_done = .TRUE.
1088             ENDIF
1089          ENDIF
1090          grid_particles(kp,jp,ip)%particles(pindex) = particle_array(n)
1091          prt_count(kp,jp,ip) = pindex
1092       ELSE
1093          IF ( jp == nys - 1 )  THEN
1094             nr_move_south = nr_move_south+1
1095             move_also_south(nr_move_south) = particle_array(n)
1096             IF ( jp == -1 )  THEN
1097                move_also_south(nr_move_south)%y =                             &
1098                   move_also_south(nr_move_south)%y + ( ny + 1 ) * dy
1099                move_also_south(nr_move_south)%origin_y =                      &
1100                   move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
1101             ENDIF
1102          ELSEIF ( jp == nyn+1 )  THEN
1103             nr_move_north = nr_move_north+1
1104             move_also_north(nr_move_north) = particle_array(n)
1105             IF ( jp == ny+1 )  THEN
1106                move_also_north(nr_move_north)%y =                             &
1107                   move_also_north(nr_move_north)%y - ( ny + 1 ) * dy
1108                move_also_north(nr_move_north)%origin_y =                      &
1109                   move_also_north(nr_move_north)%origin_y - ( ny + 1 ) * dy
1110             ENDIF
1111          ELSE
1112             WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
1113          ENDIF
1114       ENDIF
1115    ENDDO
1116
1117    RETURN
1118
1119 END SUBROUTINE Add_particles_to_gridcell
1120
1121
1122
1123
1124 SUBROUTINE lpm_move_particle
1125
1126!
1127!-- If a particle moves from one grid cell to another (on the current
1128!-- processor!), this subroutine moves the corresponding element from the
1129!-- particle array of the old grid cell to the particle array of the new grid
1130!-- cell.
1131
1132    IMPLICIT NONE
1133
1134    INTEGER(iwp)        ::  i           !:
1135    INTEGER(iwp)        ::  ip          !:
1136    INTEGER(iwp)        ::  j           !:
1137    INTEGER(iwp)        ::  jp          !:
1138    INTEGER(iwp)        ::  k           !:
1139    INTEGER(iwp)        ::  kp          !:
1140    INTEGER(iwp)        ::  n           !:
1141    INTEGER(iwp)        ::  np_old_cell !:
1142    INTEGER(iwp)        ::  n_start     !:
1143    INTEGER(iwp)        ::  pindex      !:
1144
1145    LOGICAL             ::  pack_done   !:
1146
1147    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !:
1148
1149    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
1150
1151    DO  ip = nxl, nxr
1152       DO  jp = nys, nyn
1153          DO  kp = nzb+1, nzt
1154
1155             np_old_cell = prt_count(kp,jp,ip)
1156             IF ( np_old_cell <= 0 )  CYCLE
1157             particles_old_cell => grid_particles(kp,jp,ip)%particles(1:np_old_cell)
1158             n_start = -1
1159             
1160             DO  n = 1, np_old_cell
1161                i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
1162                j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
1163                k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
1164!
1165!--             Check, if particle has moved to another grid cell.
1166                IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
1167!
1168!--                The particle has moved to another grid cell. Now check, if
1169!--                particle stays on the same processor.
1170                   IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.      &
1171                        j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
1172!
1173!--                   If the particle stays on the same processor, the particle
1174!--                   will be added to the particle array of the new processor.
1175                      number_of_particles = prt_count(k,j,i)
1176                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1177
1178                      pindex = prt_count(k,j,i)+1
1179                      IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )  &
1180                      THEN
1181                         n_start = n
1182                         EXIT
1183                      ENDIF
1184
1185                      grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
1186                      prt_count(k,j,i) = pindex
1187
1188                      particles_old_cell(n)%particle_mask = .FALSE.
1189                   ENDIF
1190                ENDIF
1191             ENDDO
1192
1193             IF ( n_start .GE. 0 )  THEN
1194                pack_done = .FALSE.
1195                DO  n = n_start, np_old_cell
1196                   i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
1197                   j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
1198                   k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
1199                   IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
1200!
1201!--                   Particle is in different box
1202                      IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.   &
1203                           j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
1204!
1205!--                      Particle stays on processor
1206                         number_of_particles = prt_count(k,j,i)
1207                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1208
1209                         pindex = prt_count(k,j,i)+1
1210                         IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
1211                         THEN
1212                            IF ( pack_done )  THEN
1213                               CALL realloc_particles_array(i,j,k)
1214                               pindex = prt_count(k,j,i)+1
1215                            ELSE
1216                               CALL lpm_pack_arrays
1217                               prt_count(k,j,i) = number_of_particles
1218
1219                               pindex = prt_count(k,j,i)+1
1220!
1221!--                            If number of particles in the new grid box
1222!--                            exceeds its allocated memory, the particle array
1223!--                            will be reallocated
1224                               IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
1225                                  CALL realloc_particles_array(i,j,k)
1226                                  pindex = prt_count(k,j,i)+1
1227                               ENDIF
1228
1229                               pack_done = .TRUE.
1230                            ENDIF
1231                         ENDIF
1232
1233                         grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
1234                         prt_count(k,j,i) = pindex
1235
1236                         particles_old_cell(n)%particle_mask = .FALSE.
1237                      ENDIF
1238                   ENDIF
1239                ENDDO
1240             ENDIF
1241          ENDDO
1242       ENDDO
1243    ENDDO
1244
1245    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'stop' )
1246
1247    RETURN
1248
1249 END SUBROUTINE lpm_move_particle
1250
1251 SUBROUTINE realloc_particles_array (i,j,k,size_in)
1252
1253    IMPLICIT NONE
1254
1255    INTEGER(iwp), INTENT(in)                       ::  i              !:
1256    INTEGER(iwp), INTENT(in)                       ::  j              !:
1257    INTEGER(iwp), INTENT(in)                       ::  k              !:
1258    INTEGER(iwp), INTENT(in), optional             ::  size_in        !:
1259
1260    INTEGER(iwp)                                   :: old_size        !:
1261    INTEGER(iwp)                                   :: new_size        !:
1262    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !:
1263    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !:
1264
1265
1266    old_size = SIZE(grid_particles(k,j,i)%particles)
1267
1268    IF ( PRESENT(size_in) )   THEN
1269       new_size = size_in
1270    ELSE
1271       new_size = old_size * ( 1.0 + alloc_factor / 100.0 )
1272    ENDIF
1273
1274    new_size = MAX( new_size, min_nr_particle )
1275
1276    IF ( old_size <= 500 )  THEN
1277
1278       tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
1279
1280       DEALLOCATE(grid_particles(k,j,i)%particles)
1281       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1282
1283       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_s(1:old_size)
1284       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1285
1286    ELSE
1287
1288       ALLOCATE(tmp_particles_d(new_size))
1289       tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
1290
1291       DEALLOCATE(grid_particles(k,j,i)%particles)
1292       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1293
1294       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_d(1:old_size)
1295       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1296
1297       DEALLOCATE(tmp_particles_d)
1298
1299    ENDIF
1300    particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1301
1302    RETURN
1303 END SUBROUTINE realloc_particles_array
1304
1305END MODULE lpm_exchange_horiz_mod
Note: See TracBrowser for help on using the repository browser.