source: palm/trunk/SOURCE/pmc_particle_interface.f90 @ 4630

Last change on this file since 4630 was 4629, checked in by raasch, 4 years ago

support for MPI Fortran77 interface (mpif.h) removed

  • Property svn:keywords set to Id
File size: 46.1 KB
RevLine 
[2801]1MODULE pmc_particle_interface
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
7! terms of the GNU General Public License as published by the Free Software
8! Foundation, either version 3 of the License, or (at your option) any later
9! version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
[4360]18! Copyright 1997-2020 Leibniz Universitaet Hannover
[2801]19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
[4629]23!
24!
[2801]25! Former revisions:
[4629]26! -----------------
[2809]27! $Id: pmc_particle_interface.f90 4629 2020-07-29 09:37:56Z suehring $
[4629]28! support for MPI Fortran77 interface (mpif.h) removed
29!
30! 4444 2020-03-05 15:59:50Z raasch
[4444]31! bugfix: preprocessor directives for serial mode added
32!
33! 4360 2020-01-07 11:25:50Z suehring
[4182]34! Corrected "Former revisions" section
35!
36! 4043 2019-06-18 16:59:00Z schwenkel
[4043]37! Remove min_nr_particle
38!
39! 4017 2019-06-06 12:16:46Z schwenkel
[3948]40! coarse bound renamed as parent_bound and icl, icr, jcs, jcn as ipl, ipr, jps, jpn.
41!
42! 3883 2019-04-10 12:51:50Z hellstea
[3883]43! Function get_number_of_childs renamed to get_number_of_children and cg
44! renamed to pg according to their definitions in pmc_interface_mod
45!
46! 3655 2019-01-07 16:51:22Z knoop
[3241]47! unused variables removed
[4182]48!
49! Initial Version (by Klaus Ketelsen)
50!
[3241]51!
[2801]52! Description:
53! ------------
54! Introduce particle transfer in nested models.
55! Interface to palm lpm model to handel particle transfer between parent and
56! child model.
57!------------------------------------------------------------------------------!
[4444]58#if defined( __parallel )
[2801]59
60   USE, INTRINSIC ::  ISO_C_BINDING
[2841]61
[2801]62   USE MPI
63
64   USE kinds
65
66   USE pegrid,                                                                 &
67       ONLY: myidx,myidy
68
69   USE indices,                                                                &
70       ONLY: nx, ny, nxl, nxr, nys, nyn, nzb, nzt, nxlg, nxrg, nysg, nyng, nbgp
71
72   USE grid_variables,                                                         &
73       ONLY:  dx, dy
74
75   USE arrays_3d,                                                              &
76        ONLY: zw
77
78   USE control_parameters,                                                     &
79       ONLY:  message_string
80
81   USE particle_attributes,                                                    &
82       ONLY:  prt_count, particles, grid_particles,                            &
83              particle_type, number_of_particles, zero_particle,               &
[4043]84              ibc_par_t, ibc_par_lr, ibc_par_ns, alloc_factor
[2801]85
[4017]86!    USE lpm_pack_and_sort_mod
[2801]87
[2967]88#if defined( __parallel )
[2801]89   USE pmc_general,                                                            &
90       ONLY: pedef
91
92   USE pmc_parent,                                                             &
93       ONLY: children, pmc_s_fillbuffer, pmc_s_getdata_from_buffer,            &
94             pmc_s_get_child_npes
95
96   USE pmc_child,                                                              &
97       ONLY:  me, pmc_c_getbuffer, pmc_c_putbuffer
98
99   USE pmc_interface,                                                          &
[3883]100       ONLY: cpl_id, get_number_of_children, nr_part, part_adr, nested_run,    &
101             get_childid, get_child_edges, nr_partc, part_adrc,                &
[3948]102             parent_bound, coord_x, coord_y, pg, get_child_gridspacing,        &
[2801]103             lower_left_coord_x, lower_left_coord_y
104
105    USE pmc_handle_communicator,                                               &
106        ONLY:  pmc_parent_for_child
107
108   USE pmc_mpi_wrapper,                                                        &
109       ONLY: pmc_send_to_parent, pmc_recv_from_child
110
[2967]111#endif
[2801]112
113   IMPLICIT NONE
[2841]114
115
[2801]116   PRIVATE
117   SAVE
118
119   TYPE  coarse_particle_def
120      INTEGER(iwp) ::  nr_particle
121     
122      TYPE(particle_type),ALLOCATABLE,DIMENSION(:) ::  parent_particles
123   END TYPE  coarse_particle_def
124
125   INTEGER(iwp),PARAMETER ::  Min_particles_per_column = 100      !<
126   INTEGER(iwp),PARAMETER ::  max_nr_particle_in_rma_win = 100000 !<
127
128   INTEGER(iwp) ::  nr_fine_in_coarse  !< Number of fine grid cells in coarse grid (one direction)
129   INTEGER(iwp) ::  particle_win_child !<
130
131   INTEGER(iwp),ALLOCATABLE,DIMENSION(:) ::  particle_win_parent !<
132   
133   TYPE(C_PTR), ALLOCATABLE,DIMENSION(:) ::  buf_ptr !<
134   
135   TYPE(particle_type), DIMENSION(:),POINTER ::  particle_in_win !<
136
137   TYPE(coarse_particle_def),ALLOCATABLE,DIMENSION(:,:) ::  coarse_particles !<
138
139   
140
141   INTERFACE pmcp_g_init
142      MODULE PROCEDURE pmcp_g_init
143   END  INTERFACE pmcp_g_init
144
145   INTERFACE pmcp_g_alloc_win
146      MODULE PROCEDURE pmcp_g_alloc_win
147   END  INTERFACE pmcp_g_alloc_win
148
149   INTERFACE pmcp_c_get_particle_from_parent
150      MODULE PROCEDURE pmcp_c_get_particle_from_parent
151   END  INTERFACE pmcp_c_get_particle_from_parent
152
153   INTERFACE pmcp_c_send_particle_to_parent
154      MODULE PROCEDURE pmcp_c_send_particle_to_parent
155   END  INTERFACE pmcp_c_send_particle_to_parent
156
157   INTERFACE pmcp_p_fill_particle_win
158      MODULE PROCEDURE pmcp_p_fill_particle_win
159   END  INTERFACE pmcp_p_fill_particle_win
160
161   INTERFACE pmcp_p_empty_particle_win
162      MODULE PROCEDURE pmcp_p_empty_particle_win
163   END  INTERFACE pmcp_p_empty_particle_win
164
165   INTERFACE pmcp_g_print_number_of_particles
166      MODULE PROCEDURE pmcp_g_print_number_of_particles
167   END  INTERFACE pmcp_g_print_number_of_particles
168
169   INTERFACE pmcp_p_delete_particles_in_fine_grid_area
170      MODULE PROCEDURE pmcp_p_delete_particles_in_fine_grid_area
171   END  INTERFACE pmcp_p_delete_particles_in_fine_grid_area
172
173   PUBLIC pmcp_g_init, pmcp_g_alloc_win, pmcp_c_get_particle_from_parent
174   PUBLIC pmcp_c_send_particle_to_parent, pmcp_p_fill_particle_win, pmcp_g_print_number_of_particles
175   PUBLIC pmcp_p_empty_particle_win, pmcp_p_delete_particles_in_fine_grid_area
176
177 CONTAINS
178
179!------------------------------------------------------------------------------!
180! Description:
181! ------------
182!> general routine:
183!> Initializing actions of the particle interface
184!> check particle boundary conditions for the child models
185!------------------------------------------------------------------------------!
186 SUBROUTINE pmcp_g_init
187 
188    IMPLICIT NONE
189 
190    INTEGER(iwp) ::  nr_childs !< Number of child models of the current model
191
[2967]192#if defined( __parallel )
193
[3883]194    nr_childs = get_number_of_children()
[2801]195!
196!-- Check if the current model has child models
197    IF ( nr_childs > 0 )   THEN
198       ALLOCATE( nr_part(nysg:nyng, nxlg:nxrg) )
199       ALLOCATE( part_adr(nysg:nyng, nxlg:nxrg) )
200       nr_part = 0
201       part_adr   = 0
202    ENDIF
203
204!
205!-- Set the boundary conditions to nested for all non root (i.e child) models
206    IF ( cpl_id > 1 )  THEN
207   
208       IF ( ibc_par_t /= 3 )   THEN
209          ibc_par_t  = 3
[3046]210          message_string = 'In Child model:  ibc_par_t is automatically set to nested '
211          CALL message( 'pmcp_g_init ', 'PA0477', 0, 1, 0, 6, 0 )
[2801]212       ENDIF
213   
214       IF ( ibc_par_lr /= 3 )   THEN
215          ibc_par_lr = 3
[3046]216          message_string = 'In Child model:  ibc_par_lr is automatically set to nested '
217          CALL message( 'pmcp_g_init ', 'PA0478', 0, 1, 0, 6, 0 )
[2801]218       ENDIF
219       
220       IF ( ibc_par_ns /= 3 )   THEN
221          ibc_par_ns = 3
[3046]222          message_string = 'In Child model:  ibc_par_ns is automatically set to nested '
223          CALL message( 'pmcp_g_init ', 'PA0479', 0, 1, 0, 6, 0 )
[2801]224       ENDIF
225       
226    ENDIF
227
[2967]228#endif
[2801]229 END SUBROUTINE pmcp_g_init
230!------------------------------------------------------------------------------!
231! Description:
232! ------------
233!> general routine:
234!> allocate the MPI windows
235!------------------------------------------------------------------------------!
236 SUBROUTINE pmcp_g_alloc_win
237 
238    IMPLICIT NONE
239   
240    INTEGER(iwp) ::  m         !< loop index
241    INTEGER(iwp) ::  ierr      !< error code
[3948]242    INTEGER(iwp) ::  ipl       !< left boundary in coarse(parent) index space
243    INTEGER(iwp) ::  ipr       !< right boundary in coarse(parent) index space
244    INTEGER(iwp) ::  jps       !< south boundary in coarse(parent) index space
245    INTEGER(iwp) ::  jpn       !< north boundary in coarse(parent) index space
[2801]246    INTEGER(iwp) ::  child_id  !< Id of a child model
247    INTEGER(iwp) ::  nr_childs !< Number of child models of the current model
248   
[2884]249    INTEGER ::  parsize !<
[2801]250    TYPE(C_PTR), SAVE ::  ptr !<
251   
252    TYPE(particle_type),DIMENSION(:),POINTER ::  win_buffer !<
253   
254    INTEGER(iwp),DIMENSION(1) ::  buf_shape !<
255
[2967]256#if defined( __parallel )
257    INTEGER(KIND=MPI_ADDRESS_KIND) ::  parsize_mpi_address_kind !<
258    INTEGER(KIND=MPI_ADDRESS_KIND) ::  winsize !<
259
[2801]260!
261!-- If the model has a parent model prepare the structures for transfer
262    IF ( cpl_id > 1 )  THEN
[2809]263
[2884]264       parsize_mpi_address_kind = STORAGE_SIZE(zero_particle)/8
[2809]265
[2884]266       CALL MPI_ALLOC_MEM( parsize_mpi_address_kind , MPI_INFO_NULL, ptr, ierr )
267       parsize = parsize_mpi_address_kind
[2801]268       buf_shape(1) = 1
269       CALL C_F_POINTER( ptr, win_buffer, buf_shape )
[2884]270       CALL MPI_WIN_CREATE( win_buffer, parsize_mpi_address_kind, parsize,     &
[2801]271                            MPI_INFO_NULL, me%intra_comm, particle_win_child,  &
272                            ierr )
273
274!
275!--    Child domain boundaries in the parent index space
[3948]276       ipl = parent_bound(1)
277       ipr = parent_bound(2)
278       jps = parent_bound(3)
279       jpn = parent_bound(4)
[2801]280
[3948]281       ALLOCATE( coarse_particles(jps:jpn,ipl:ipr) )
[2801]282
283       coarse_particles(:,:)%nr_particle = 0
284    ENDIF
285
286!
287!-- If the model has child models prepare the structures for transfer
[3883]288    nr_childs = get_number_of_children()
[2801]289    IF ( nr_childs > 0 )   THEN
290       ALLOCATE( particle_win_parent(nr_childs) )
291       ALLOCATE( buf_ptr(nr_childs) )
292       DO  m = 1, nr_childs
293          child_id = get_childid(m)
[2884]294          parsize_mpi_address_kind = STORAGE_SIZE(zero_particle)/8
295          parsize = parsize_mpi_address_kind
[2809]296
[2884]297          winsize = max_nr_particle_in_rma_win * parsize_mpi_address_kind
[2801]298          CALL MPI_ALLOC_MEM( winsize , MPI_INFO_NULL, buf_ptr(m), ierr )
299          buf_shape(1) = max_nr_particle_in_rma_win
300          CALL C_F_POINTER( buf_ptr(m), win_buffer, buf_shape )
[2884]301          CALL MPI_WIN_CREATE( win_buffer, winsize, parsize, MPI_INFO_NULL,    &
302                               children(child_id)%intra_comm,                  &
[2801]303                               particle_win_parent(m), ierr )
304          ENDDO
305    ENDIF
306
[2967]307#endif
[2801]308 END SUBROUTINE pmcp_g_alloc_win
[2967]309
310
[2801]311!------------------------------------------------------------------------------!
312! Description:
313! ------------
314!> child routine:
315!> Read/get particles out of the parent MPI window
316!------------------------------------------------------------------------------!
317 SUBROUTINE pmcp_c_get_particle_from_parent
318 
319    IMPLICIT NONE
320
321    INTEGER(iwp) ::  i    !< x grid index
[3948]322    INTEGER(iwp) ::  ipl  !< left boundary in coarse(parent) index space
[2801]323    INTEGER(iwp) ::  ierr !< error code
324    INTEGER(iwp) ::  ij   !< combined xy index for the buffer array
325    INTEGER(iwp) ::  ip   !< loop index (child PEs)
326    INTEGER(iwp) ::  j    !< y grid index
[3948]327    INTEGER(iwp) ::  jps  !< south boundary in coarse(parent) index space
[2801]328    INTEGER(iwp) ::  nr   !< number of particles to receive from a parent box
329   
[2884]330    INTEGER ::  parsize !<
[2967]331
332#if defined( __parallel )
333    TYPE(pedef), POINTER ::  ape !< TO_DO Klaus: give a description and better name of the variable
334
[2884]335    INTEGER(KIND=MPI_ADDRESS_KIND) ::  parsize_mpi_address_kind !<
[2801]336    INTEGER(KIND=MPI_ADDRESS_KIND) ::  target_disp !<
337
338    IF ( cpl_id > 1 )  THEN
339   
340       CALL pmc_c_getbuffer( particle_transfer = .TRUE. ) !Get number of particle/column and offset in RMA window xx
341
342!
343!--    Wait for buffer to fill.
344!
345!--    The parent side (in pmc_s_fillbuffer) is filling the buffer in the MPI RMA window
346!--    When the filling is complete, a MPI_BARRIER is called.
347!--    The child is not allowd to access the parent-buffer before it is completely filled
348!--    Synchronization is done implicitely in pmc_c_getbuffer and pmc_s_fillbuffer on the parent side
349
[3948]350       ipl = parent_bound(1)
351       jps = parent_bound(3)
[2801]352
353       DO  ip = 1, me%inter_npes
354
355          ape => me%pes(ip)
356
357          DO  ij = 1, ape%nrele
[3948]358              j = ape%locind(ij)%j + jps - 1
359              i = ape%locind(ij)%i + ipl - 1
[2801]360              nr = nr_partc(j,i)
361              IF ( nr > 0 )  THEN
362
363                 CALL check_and_alloc_coarse_particle (i, j, nr)
[2884]364                 parsize_mpi_address_kind = STORAGE_SIZE(zero_particle)/8
365                 parsize = parsize_mpi_address_kind
[2801]366                 target_disp = part_adrc(j,i) - 1
[2884]367                 CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , ip - 1, 0,               &
368                                    particle_win_child, ierr )
[2801]369                 CALL MPI_GET( coarse_particles(j,i)%parent_particles,         &
[2884]370                               nr * parsize, MPI_BYTE, ip - 1, target_disp,    &
371                               nr * parsize, MPI_BYTE, particle_win_child, ierr )
[2801]372                 CALL MPI_WIN_UNLOCK( ip - 1, particle_win_child, ierr )
373             ENDIF
374             coarse_particles(j,i)%nr_particle = nr
375          ENDDO
376       ENDDO
377
378       CALL c_copy_particle_to_child_grid
379    ENDIF
380
[2967]381#endif
[2801]382 END SUBROUTINE pmcp_c_get_particle_from_parent
[2967]383
384
[2801]385!------------------------------------------------------------------------------!
386! Description:
387! ------------
388!> child routine:
389!> Write/put particles into the parent MPI window
390!------------------------------------------------------------------------------!
391  SUBROUTINE pmcp_c_send_particle_to_parent
392 
393    IMPLICIT NONE
394   
395    INTEGER(iwp) ::  disp_offset            !<
396    INTEGER(iwp) ::  i                      !< x loop index
[3948]397    INTEGER(iwp) ::  ipl                    !< left boundary in coarse(parent) index space
398    INTEGER(iwp) ::  ipr                    !< right boundary in coarse(parent) index space
[2801]399    INTEGER(iwp) ::  ierr                   !< error code
400    INTEGER(iwp) ::  ij                     !< combined xy index for the buffer array
401    INTEGER(iwp) ::  ip                     !< loop index (child PEs)
402    INTEGER(iwp) ::  j                      !< y loop index
[3948]403    INTEGER(iwp) ::  jps                    !< south boundary in coarse(parent) index space
404    INTEGER(iwp) ::  jpn                    !< north boundary in coarse(parent) index space
[2801]405    INTEGER(iwp) ::  max_nr_particle_per_pe !< maximum number of particles per PE (depending on grid apect ratio)
406    INTEGER(iwp) ::  n                      !< shorter variable name for nr_fine_in_coarse
407    INTEGER(iwp) ::  nr                     !< shorter variabel name for nr_partc
408    INTEGER(iwp) ::  pe_offset              !< offset index of the current PE
409   
[2884]410    INTEGER ::  parsize !<
[2801]411   
412    REAL(wp) ::  eps=0.00001 !< used in calculations to avoid rounding errors
413    REAL(wp) ::  xx          !< number of fine grid cells inside a coarse grid cell in x-direction
414    REAL(wp) ::  yy          !< number of fine grid cells inside a coarse grid cell in y-direction
415
[2967]416 !   TYPE(particle_type) ::  dummy_part !< dummy particle (needed for size calculations)
417
418#if defined( __parallel )
[2801]419    TYPE(pedef), POINTER ::  ape !< TO_DO Klaus: give a description and better name of the variable
[2967]420
421    INTEGER(KIND=MPI_ADDRESS_KIND) ::  parsize_mpi_address_kind !<
422    INTEGER(KIND=MPI_ADDRESS_KIND) ::  target_disp !<
[2801]423   
424   
425    IF ( cpl_id > 1 )  THEN
426       CALL c_copy_particle_to_coarse_grid
427
428!
429!--    Child domain boundaries in the parent index space
430
[3948]431       ipl = parent_bound(1)
432       ipr = parent_bound(2)
433       jps = parent_bound(3)
434       jpn = parent_bound(4)
[2801]435
436       nr_partc = 0
437       
[3948]438       DO i = ipl, ipr
439          DO j = jps, jpn
[2801]440             nr_partc(j,i) = coarse_particles(j,i)%nr_particle
441          ENDDO
442       ENDDO
443       part_adrc = 0
444
445!
446!--    compute number of fine grid cells in coarse grid (one direction)
[3883]447       xx = ( pg%dx + eps ) / dx ! +eps to avoid rounding error
448       yy = ( pg%dy + eps ) / dy
[2801]449       nr_fine_in_coarse = MAX( INT(xx), INT(yy) )
450
[3883]451       IF ( MOD( coord_x(0), pg%dx ) /= 0.0 .OR. MOD( coord_y(0), pg%dy ) /= 0.0 ) THEN
[2801]452          nr_fine_in_coarse = nr_fine_in_coarse + 1
453       ENDIF
454
455!
456!--    Assign a number to my child PE to select different areas in the RMA window on server side
457!--    With this number a square of child PEs is defined which share the same coarse grid cells
458
459       n           = nr_fine_in_coarse ! local variable n to make folloing statements shorter
460       pe_offset   = MOD( myidx, n ) * n + MOD( myidy, n)
461       max_nr_particle_per_pe = max_nr_particle_in_rma_win / ( n * n )
462       disp_offset            = pe_offset * max_nr_particle_per_pe
[2884]463       parsize_mpi_address_kind = STORAGE_SIZE(zero_particle)/8
464       parsize = parsize_mpi_address_kind
[2801]465       DO  ip = 1, me%inter_npes
466
467          ape => me%pes(ip)
468
469          target_disp = disp_offset
470          DO ij = 1, ape%nrele
[3948]471             j  = ape%locind(ij)%j + jps - 1
472             i  = ape%locind(ij)%i + ipl - 1
[2801]473             nr = nr_partc(j,i)
474             IF( nr > 0 ) THEN
475                IF ( target_disp + nr - disp_offset >= max_nr_particle_per_pe ) THEN
476                   WRITE(9,*) 'RMA window too small on child ',                &
477                              target_disp + nr - disp_offset,                  &
478                              max_nr_particle_per_pe, max_nr_particle_in_rma_win
479                   message_string = 'RMA window too small on child'
[3046]480                   CALL message( 'pmci_create_child_arrays', 'PA0480', 3, 2, 0, 6, 0 )
[2801]481                ENDIF
[2884]482                CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , ip - 1, 0,                &
483                                   particle_win_child, ierr )
[2801]484                CALL MPI_PUT( coarse_particles(j,i)%parent_particles,          &
[2884]485                              nr * parsize, MPI_BYTE, ip - 1, target_disp,     &
486                              nr * parsize, MPI_BYTE,  particle_win_child, ierr )
[2801]487                CALL MPI_WIN_UNLOCK( ip - 1, particle_win_child, ierr )
488                part_adrc(j,i) = target_disp + 1
489                target_disp    = target_disp + nr
490             ENDIF
491          ENDDO
492       ENDDO
493
494       CALL pmc_c_putbuffer ( particle_transfer = .TRUE. )   !Send new number of particle/column and offset to parent
495
496    ENDIF
497
[2967]498#endif
[2801]499 END SUBROUTINE pmcp_c_send_particle_to_parent
[2967]500
501
[2801]502!------------------------------------------------------------------------------!
503! Description:
504! ------------
505!> parent routine:
506!> write particles into the MPI window
507!------------------------------------------------------------------------------!
508 SUBROUTINE pmcp_p_fill_particle_win
509 
510    IMPLICIT NONE
511
512    LOGICAL      ::  active_particle !< Particles located in the fine/child grid area are marked as active (to be transferred)
513    LOGICAL,SAVE ::  lfirst = .TRUE. !<
514   
515    INTEGER(iwp) ::  child_id            !< id of the child model
516    INTEGER(iwp) ::  i                   !< x grid box index
517    INTEGER(iwp) ::  ij                  !< combined xy index for the buffer array
518    INTEGER(iwp) ::  ip                  !< loop index (child PEs)
519    INTEGER(iwp) ::  j                   !< y grid box index
520    INTEGER(iwp) ::  k                   !< z grid box index
521    INTEGER(iwp) ::  m                   !< loop index (number of childs)
522    INTEGER(iwp) ::  n                   !< loop index (number of particles)
523    INTEGER(iwp) ::  nr_part_col         !< Number of particles to transfer per column
524    INTEGER(iwp) ::  number_of_particles !<
525    INTEGER(iwp) ::  pindex              !<
526    INTEGER(iwp) ::  tot_particle_count  !< Total number of particles per child
527   
528    REAL(wp) ::  dx_child   !< child grid spacing
529    REAL(wp) ::  dy_child   !< child grid spacing
530    REAL(wp) ::  dz_child   !< child grid spacing
531    REAL(wp) ::  ny_coord   !< north coordinate of child grid
532    REAL(wp) ::  ny_coord_b !< north coordinate of child grid boundary
533    REAL(wp) ::  lx_coord   !< left coordinate of child grid
534    REAL(wp) ::  lx_coord_b !< left coordinate of child grid boundary
535    REAL(wp) ::  rx_coord   !< right coordinate of child grid
536    REAL(wp) ::  rx_coord_b !< right coordinate of child grid boundary
537    REAL(wp) ::  sy_coord   !< south coordinate of child grid
538    REAL(wp) ::  sy_coord_b !< south coordinate of child grid boundary
539    REAL(wp) ::  uz_coord   !< upper coordinate of child grid
540    REAL(wp) ::  uz_coord_b !< upper coordinate of child grid boundary
541    REAL(wp) ::  x          !< particle position
542    REAL(wp) ::  y          !< particle position
543    REAL(wp) ::  z          !< particle position
544   
545    INTEGER(iwp),DIMENSION(1) ::  buf_shape !<
546   
[2967]547#if defined( __parallel )
[2801]548    TYPE(pedef), POINTER ::  ape !< TO_DO Klaus: give a description and better name of the variable
549
[3883]550    DO  m = 1, get_number_of_children()
[2801]551
552       child_id = get_childid(m)
553
554       CALL get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b,    &
555                                sy_coord, sy_coord_b, ny_coord, ny_coord_b,    &
556                                uz_coord, uz_coord_b) 
557
558       CALL get_child_gridspacing( m, dx_child, dy_child, dz_child )
559
560       IF ( lfirst )   THEN
561          WRITE(9,'(a,5f10.2)') 'edges          ',lx_coord,rx_coord,sy_coord,ny_coord,uz_coord
562          WRITE(9,'(a,5f10.2)') 'edges boundary ',lx_coord_b,rx_coord_b,sy_coord_b,ny_coord_b,uz_coord_b
563          WRITE(9,'(a,5f10.2)') 'child spacing  ',dx_child, dy_child, dz_child,lower_left_coord_x,lower_left_coord_y
564       ENDIF
565!
566!--    reset values for every child
567       tot_particle_count = 0 
568       nr_part  = 0           
569       part_adr = 0
570       pindex   = 1
571
572       buf_shape(1) = max_nr_particle_in_rma_win
573       CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
574
575       DO  ip = 1, children(child_id)%inter_npes
576
577          ape => children(child_id)%pes(ip)
578
579          nr_part_col   = 0 
580         
581          DO  ij = 1, ape%nrele
582             
583!
584!--          Inside the PMC adressing of 3d arrays starts with 1
585             i = ape%locind(ij)%i + nxl - nbgp - 1
586             j = ape%locind(ij)%j + nys - nbgp - 1
[3241]587             nr_part_col = 0   ! Number of particles to transfer per column
[2801]588             part_adr(j,i) = pindex
589             
590             DO  k = nzb + 1, nzt
591                number_of_particles = prt_count(k,j,i)
592               
593                IF ( number_of_particles <= 0 )  CYCLE
594               
595                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
596
597                ! Select particles within boundary area
598
599                DO n = 1, number_of_particles
600                   x = particles(n)%x
601                   y = particles(n)%y
602                   z = particles(n)%z
603!
604!--                check if the particle is located in the fine grid area
605                   active_particle = ((x > lx_coord .AND. x < rx_coord) .AND.  & 
606                                      (y > sy_coord .AND. y < ny_coord) .AND.  &
607                                      (z > 0.000000001 .AND. z < uz_coord))
608                   IF ( active_particle .AND. particles(n)%particle_mask )  THEN
609                     
610                      particle_in_win(pindex) = particles(n)
611!
612!--                   Change particle positions and origin relative to global origin
613                      particle_in_win(pindex)%x = particle_in_win(pindex)%x + lower_left_coord_x
614                      particle_in_win(pindex)%y = particle_in_win(pindex)%y + lower_left_coord_y
615                      particle_in_win(pindex)%origin_x = particle_in_win(pindex)%origin_x + lower_left_coord_x
616                      particle_in_win(pindex)%origin_y = particle_in_win(pindex)%origin_y + lower_left_coord_y
617
618                      tot_particle_count = tot_particle_count + 1
619                      nr_part_col        = nr_part_col + 1
620                      pindex             = pindex + 1
621                      IF ( pindex > max_nr_particle_in_rma_win ) THEN
622                         WRITE(9,*) 'RMA window too small on parent ',pindex, max_nr_particle_in_rma_win
623                         message_string = 'RMA window too small on parent'
624                         CALL message( 'pmci_create_child_arrays', 'PA0481', 3, 2, 0, 6, 0 )   ! PA number has to be adjusted
625                     ENDIF
626                   END IF
627                ENDDO
628             ENDDO
629             nr_part(j,i) = nr_part_col
630          ENDDO
631       ENDDO
632
633       CALL pmc_s_fillbuffer( child_id, particle_transfer = .TRUE. )
634    ENDDO
635
636    lfirst = .FALSE.
637
[2967]638#endif
[2801]639 END SUBROUTINE pmcp_p_fill_particle_win
[2967]640
[2801]641 
642!------------------------------------------------------------------------------!
643! Description:
644! ------------
645!> parent routine:
646!> delete particles from the MPI window
647!------------------------------------------------------------------------------!
648 SUBROUTINE pmcp_p_empty_particle_win
[2967]649
[2801]650    IMPLICIT NONE
651
652    INTEGER(iwp) ::  child_id           !< model id of the child
653    INTEGER(iwp) ::  ip                 !< loop index (child PEs)
654    INTEGER(iwp) ::  m                  !< loop index (number of childs)
655
656    INTEGER(iwp),DIMENSION(1) ::  buf_shape !<
657
[2967]658#if defined( __parallel )
[3883]659    DO  m = 1, get_number_of_children()
[2801]660
661       child_id = get_childid(m)
662
663       buf_shape(1) = max_nr_particle_in_rma_win
664       CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
665
666!
[3241]667!--    In some cells of the coarse grid, there are contributions from more than
668!--    one child process. Therefore p_copy_particle_to_org_grid is done for one
669!--    child process per call
670       DO ip = 1, pmc_s_get_child_npes( child_id )
[2801]671
672          nr_part  = 0
673          part_adr = 0
674
[3241]675          CALL pmc_s_getdata_from_buffer( child_id, particle_transfer = .TRUE.,&
676                                          child_process_nr = ip )
677          CALL p_copy_particle_to_org_grid( m )
[2801]678       ENDDO
679
680    ENDDO
681
[2967]682#endif
[2801]683 END SUBROUTINE pmcp_p_empty_particle_win
[2967]684
685
[2801]686!------------------------------------------------------------------------------!
687! Description:
688! ------------
689!> parent routine:
690!> After the transfer mark all parent particles that are still inside on of the
691!> child areas for deletion.
692!------------------------------------------------------------------------------!
693 SUBROUTINE pmcp_p_delete_particles_in_fine_grid_area
694
695    IMPLICIT NONE
696
697    LOGICAL ::  to_delete !< particles outside of model domain are marked as to_delete
698   
699    INTEGER(iwp) ::  i !< loop index (x grid)
700    INTEGER(iwp) ::  j !< loop index (y grid)
701    INTEGER(iwp) ::  k !< loop index (z grid)
702    INTEGER(iwp) ::  m !< loop index (number of particles)
703    INTEGER(iwp) ::  n !< loop index (number of childs)
704   
705    REAL(wp) ::  dx_child   !< child grid spacing
706    REAL(wp) ::  dy_child   !< child grid spacing
707    REAL(wp) ::  dz_child   !< child grid spacing
708    REAL(wp) ::  ny_coord   !< north coordinate of child grid
709    REAL(wp) ::  ny_coord_b !< north coordinate of child grid boundary
710    REAL(wp) ::  lx_coord   !< left coordinate of child grid
711    REAL(wp) ::  lx_coord_b !< left coordinate of child grid boundary
712    REAL(wp) ::  rx_coord   !< right coordinate of child grid
713    REAL(wp) ::  rx_coord_b !< right coordinate of child grid boundary
714    REAL(wp) ::  sy_coord   !< south coordinate of child grid
715    REAL(wp) ::  sy_coord_b !< south coordinate of child grid boundary
716    REAL(wp) ::  uz_coord   !< upper coordinate of child grid
717    REAL(wp) ::  uz_coord_b !< upper coordinate of child grid boundary
718    REAL(wp) ::  x          !< particle position
719    REAL(wp) ::  y          !< particle position
720    REAL(wp) ::  z          !< particle position
721   
[2967]722#if defined( __parallel )
[3883]723    DO  m = 1, get_number_of_children()
[2801]724       CALL get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b,    &
725                                sy_coord, sy_coord_b, ny_coord, ny_coord_b,    &
726                                uz_coord, uz_coord_b )
727
728
729       CALL get_child_gridspacing( m, dx_child, dy_child, dz_child ) 
730
731       DO i = nxl, nxr
732          DO j = nys, nyn
733             DO k = nzb, nzt
734                number_of_particles = prt_count(k,j,i)
735
736                IF ( number_of_particles == 0 ) CYCLE
737
738                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
739
740                DO n = 1, number_of_particles
741                   x = particles(n)%x
742                   y = particles(n)%y
743                   z = particles(n)%z
744
745                   to_delete = ((x > lx_coord .AND. x < rx_coord) .AND.        &
746                                (y > sy_coord .AND. y < ny_coord) .AND.        &
747                                (z > 0.000000001 .AND. z < uz_coord))
748
749                   IF ( to_delete ) THEN
750                      particles(n)%particle_mask = .FALSE.
751                   ENDIF
752                ENDDO
753             ENDDO
754          ENDDO
755       ENDDO
756
757    ENDDO
758
[2967]759#endif
[2801]760 END SUBROUTINE pmcp_p_delete_particles_in_fine_grid_area
[2967]761
762
[2801]763!------------------------------------------------------------------------------!
764! Description:
765! ------------
766!> general routine:
767!> print the total number of of the current model and its child models into a file
768!------------------------------------------------------------------------------!
769 SUBROUTINE pmcp_g_print_number_of_particles (my_time,local_nr_particles)
770 
771    USE pegrid,                                                                &
772        ONLY: myid
773 
774    IMPLICIT NONE
775
776    INTEGER(iwp),INTENT(IN) ::  local_nr_particles !<
777   
778    REAL(wp),INTENT(IN) ::  my_time !<
779     
780    LOGICAL, SAVE ::  is_file_open=.FALSE. !<
781   
782    INTEGER(iwp) ::  child_id           !<
783    INTEGER(iwp) ::  child_nr_particles !< total number of particles in all child models
784    INTEGER(iwp) ::  ierr               !< error code
785    INTEGER(iwp) ::  m                  !< loop index (number of childs)
786   
787    CHARACTER(LEN=16) ::  fname !< filename
788   
789    INTEGER(iwp),DIMENSION(2) ::  ivalr !< integer value to be received
790    INTEGER(iwp),DIMENSION(2) ::  ivals !< integer value to be send
791   
[2967]792#if defined( __parallel )
[2801]793    child_nr_particles = 0
794    IF ( myid == 0 )  THEN
795       IF ( cpl_id > 1 )  THEN
796          ivals(1) = local_nr_particles
797          CALL pmc_send_to_parent( ivals, 1, 0, 400, ierr )
798       ENDIF
799       DO  m = 1, SIZE( pmc_parent_for_child ) - 1
800
801          child_id = pmc_parent_for_child(m)
802          CALL pmc_recv_from_child( child_id, ivalr, 1, 0, 400, ierr )
803          child_nr_particles = child_nr_particles + ivalr(1)
804       ENDDO
805
806       IF ( SIZE( pmc_parent_for_child ) > 1 ) THEN
807          IF ( .NOT. is_file_open )  THEN !kk muss noch auf file_open umgestellt werden
808             WRITE(fname,'(a,i2.2)') 'nr_particles_',cpl_id
809             OPEN (333,file = fname)
810             is_file_open = .true.
811          ENDIF
812          WRITE(333,'(f12.2,3i10)') my_time,local_nr_particles + child_nr_particles,local_nr_particles,child_nr_particles
813       ENDIF
814    ENDIF
815
[2967]816#endif
[2801]817 END SUBROUTINE pmcp_g_print_number_of_particles
818
[2967]819
[2801]820!------------------------------------------------------------------------------!
821!------------------------------------------------------------------------------!
822! Private subroutines
823!------------------------------------------------------------------------------!
824!------------------------------------------------------------------------------!
825
826!------------------------------------------------------------------------------!
827! Description:
828! ------------
829!> child routine
830!> update the size of the local buffer (coarse_particles)
831!------------------------------------------------------------------------------!
832 SUBROUTINE check_and_alloc_coarse_particle (ic,jc,nr,with_copy)
833   
834    IMPLICIT NONE
835   
836    INTEGER(iwp),INTENT(IN) ::  ic !< coarse x grid index
837    INTEGER(iwp),INTENT(IN) ::  jc !< coarse y grid index
838    INTEGER(iwp),INTENT(IN) ::  nr !< number of particles to be transferred/stored in local buffer
839   
840    LOGICAL,INTENT(IN),OPTIONAL ::  with_copy !< copy particles in buffer? or reallocate empty buffer
841
842    LOGICAL :: with_copy_lo !< local variable of with copy
843   
844    INTEGER(iwp) ::  new_size !< new size of the local buffer
845    INTEGER(iwp) ::  old_size !< old size of the local buffer
846   
847    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles_d !<
848
[2967]849#if defined( __parallel )
[2801]850    with_copy_lo = .FALSE.
851    IF ( PRESENT( with_copy ) ) with_copy_lo = with_copy
852
853    IF ( .NOT. ALLOCATED( coarse_particles(jc,ic)%parent_particles ) ) THEN
854       new_size = MAX( nr, Min_particles_per_column )
855       ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
856    ELSEIF ( nr > SIZE( coarse_particles(jc,ic)%parent_particles ) ) THEN
857
858       old_size = SIZE( coarse_particles(jc,ic)%parent_particles )
859       new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
860       new_size = MAX( nr, new_size, old_size + Min_particles_per_column )
861
862!
863!--    Copy existing particles to new particle buffer
864       IF ( with_copy_lo ) THEN
865          ALLOCATE( tmp_particles_d(old_size) )
866          tmp_particles_d(1:old_size) = coarse_particles(jc,ic)%parent_particles
867
868          DEALLOCATE( coarse_particles(jc,ic)%parent_particles )
869          ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
870
871          coarse_particles(jc,ic)%parent_particles(1:old_size)          = tmp_particles_d(1:old_size)
872          coarse_particles(jc,ic)%parent_particles(old_size+1:new_size) = zero_particle
873
874          DEALLOCATE( tmp_particles_d )
875!
876!--    allocate or reallocate an empty buffer
877       ELSE
878          DEALLOCATE( coarse_particles(jc,ic)%parent_particles )
879          ALLOCATE( coarse_particles(jc,ic)%parent_particles(new_size) )
880       ENDIF
881    ENDIF
882
[2967]883#endif
[2801]884 END SUBROUTINE check_and_alloc_coarse_particle
885
[2967]886
[2801]887!------------------------------------------------------------------------------!
888! Description:
889! ------------
890!> child routine:
891!> copy/sort particles out of the local buffer into the respective grid boxes
892!------------------------------------------------------------------------------!
893 SUBROUTINE c_copy_particle_to_child_grid
894 
895    IMPLICIT NONE
896 
897    INTEGER(iwp) ::  ic  !< coarse x grid index
[3948]898    INTEGER(iwp) ::  ipl !< left boundary in coarse(parent) index space
899    INTEGER(iwp) ::  ipr !< right boundary in coarse(parent) index space
[2801]900    INTEGER(iwp) ::  ip  !< x grid index
901    INTEGER(iwp) ::  jc  !< coarse y grid index
[3948]902    INTEGER(iwp) ::  jpn !< north boundary in coarse(parent) index space
903    INTEGER(iwp) ::  jps !< south boundary in coarse(parent) index space
[2801]904    INTEGER(iwp) ::  jp  !< y grid index
905    INTEGER(iwp) ::  kp  !< z grid index
906    INTEGER(iwp) ::  n   !< loop index (number of particles)
907    INTEGER(iwp) ::  nr  !< short variable for name number or particles
908   
909    REAL(wp) ::  xc  !< child x coordinate
910    REAL(wp) ::  xoc !< child x origin
911    REAL(wp) ::  yc  !< child y coordinate
912    REAL(wp) ::  yoc !< child y origin
913    REAL(wp) ::  zc  !< child z coordinate
914
[2967]915#if defined( __parallel )
[2801]916!
917!-- Child domain boundaries in the parent index space
[3948]918    ipl = parent_bound(1)
919    ipr = parent_bound(2)
920    jps = parent_bound(3)
921    jpn = parent_bound(4)
[2801]922
[3948]923    DO ic = ipl, ipr
924       DO jc = jps, jpn
[2801]925          nr = coarse_particles(jc,ic)%nr_particle
926
927          IF ( nr > 0 ) THEN
928             DO n = 1, nr
929                xc =  coarse_particles(jc,ic)%parent_particles(n)%x-lower_left_coord_x
930                yc =  coarse_particles(jc,ic)%parent_particles(n)%y-lower_left_coord_y
931                zc =  coarse_particles(jc,ic)%parent_particles(n)%z
932                xoc = coarse_particles(jc,ic)%parent_particles(n)%origin_x-lower_left_coord_x
933                yoc = coarse_particles(jc,ic)%parent_particles(n)%origin_y-lower_left_coord_y
934                ip = xc / dx
935                jp = yc / dy
936                kp = nzt
937                DO WHILE ( zw(kp-1) > zc .AND. kp > nzb + 1 )         ! kk search loop has to be optimzed !!!
938                   kp = kp - 1
939                ENDDO 
940
941                IF ( ip >= nxl .AND. ip <= nxr .AND. jp >= nys .AND. jp <= nyn ) THEN
942                   prt_count(kp,jp,ip) = prt_count(kp,jp,ip) + 1
943                   IF ( prt_count(kp,jp,ip) > SIZE( grid_particles(kp,jp,ip)%particles ) ) THEN
[4017]944                      CALL pmc_realloc_particles_array( ip, jp, kp )
[2801]945                   ENDIF
946                   coarse_particles(jc,ic)%parent_particles(n)%x = xc                   ! Adjust coordinates to child grid
[3241]947                   coarse_particles(jc,ic)%parent_particles(n)%y = yc
[2801]948                   coarse_particles(jc,ic)%parent_particles(n)%origin_x = xoc           ! Adjust origins to child grid
949                   coarse_particles(jc,ic)%parent_particles(n)%origin_y = yoc
950                   grid_particles(kp,jp,ip)%particles(prt_count(kp,jp,ip)) = coarse_particles(jc,ic)%parent_particles(n)
951                ENDIF
952             ENDDO
953          ENDIF
954       ENDDO
955    ENDDO
956
[2967]957#endif
[2801]958 END SUBROUTINE c_copy_particle_to_child_grid
[2967]959
960
[2801]961!------------------------------------------------------------------------------!
962! Description:
963! ------------
964!> child routine:
965!> copy particles which left the model area into the local buffer
966!------------------------------------------------------------------------------!
967 SUBROUTINE c_copy_particle_to_coarse_grid
968 
969    IMPLICIT NONE
970   
971    LOGICAL ::  boundary_particle !<
972   
973    INTEGER(iwp) ::  i    !< loop index (x grid)
974    INTEGER(iwp) ::  ic   !< loop index (coarse x grid)
[3948]975    INTEGER(iwp) ::  ipl  !< left boundary in coarse(parent) index space
976    INTEGER(iwp) ::  ipr  !< left boundary in coarse(parent) index space
[2801]977    INTEGER(iwp) ::  ierr !< error code
978    INTEGER(iwp) ::  j    !< loop index (y grid)
979    INTEGER(iwp) ::  jc   !< loop index (coarse y grid)
[3948]980    INTEGER(iwp) ::  jps  !< south boundary in coarse(parent) index space
981    INTEGER(iwp) ::  jpn  !< north boundary in coarse(parent) index space
[2801]982    INTEGER(iwp) ::  k    !< loop index (z grid)
983    INTEGER(iwp) ::  n    !< loop index (number of particles)
984   
[3123]985    REAL(wp) ::  x       !< x coordinate
986    REAL(wp) ::  xo      !< x origin
987    REAL(wp) ::  x_left  !< absolute left boundary
988    REAL(wp) ::  x_right !< absolute right boundary
989    REAL(wp) ::  y       !< left boundary in coarse(parent) index space
990    REAL(wp) ::  yo      !< y origin
991    REAL(wp) ::  y_north !< absolute north boundary
992    REAL(wp) ::  y_south !< absoulte south boundary
993    REAL(wp) ::  z       !< z coordinate
[2801]994
[2967]995#if defined( __parallel )
[2801]996!
997!-- Child domain boundaries in the parent index space
998
[3948]999    ipl = parent_bound(1)
1000    ipr = parent_bound(2)
1001    jps = parent_bound(3)
1002    jpn = parent_bound(4)
[2801]1003
1004!
1005!-- absolute coordinates
1006    x_left  = coord_x(0)
1007    x_right = coord_x(nx) + dx
1008    y_south = coord_y(0)
1009    y_north = coord_y(ny) + dy
1010
1011!   Clear Particle Buffer
1012
[3948]1013    DO ic = ipl, ipr
1014       DO jc = jps, jpn
[2801]1015          coarse_particles(jc,ic)%nr_particle = 0
1016       ENDDO
1017    ENDDO
1018
1019!
1020!-- Determine particles which leave the inner area in east or west dirextion
1021!-- Compute only first (nxl) and last (nxr) loop iterration.
1022    DO  i = nxl, nxr
1023       DO  j = nys, nyn
1024          DO  k = nzb + 1, nzt
1025             number_of_particles = prt_count(k,j,i)
1026             IF ( number_of_particles <= 0 )  CYCLE
1027             particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1028             DO  n = 1, number_of_particles
1029                IF ( particles(n)%particle_mask )  THEN
1030                   x =  particles(n)%x+ lower_left_coord_x 
1031                   y =  particles(n)%y+ lower_left_coord_y 
1032                   xo = particles(n)%origin_x + lower_left_coord_x
1033                   yo = particles(n)%origin_y + lower_left_coord_y
1034                   z =  particles(n)%z
1035                   
1036                   boundary_particle = .FALSE.
1037                   boundary_particle = boundary_particle .OR. (x < x_left  .OR. x > x_right)
1038                   boundary_particle = boundary_particle .OR. (y < y_south .OR. y > y_north)
1039                   boundary_particle = boundary_particle .OR. (z > zw(nzt))
1040
1041                   IF ( boundary_particle ) THEN                     
[3883]1042                      ic = x / pg%dx                     !TODO anpassen auf Mehrfachnesting
1043                      jc = y / pg%dy
[2801]1044
[3948]1045                      IF ( ic >= ipl .AND. ic <= ipr .AND. jc >= jps .AND. jc <= jpn ) THEN
[2801]1046                         coarse_particles(jc,ic)%nr_particle = coarse_particles(jc,ic)%nr_particle + 1
1047                         CALL check_and_alloc_coarse_particle( ic, jc, coarse_particles(jc,ic)%nr_particle, with_copy=.TRUE. )
1048
1049                         coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)   = particles(n)
1050                         coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%x = x   !adapt to global coordinates
1051                         coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%y = y
1052                         coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%origin_x = xo
1053                         coarse_particles(jc,ic)%parent_particles(coarse_particles(jc,ic)%nr_particle)%origin_y = yo
1054!--                      Mark particle as deleted after copying it to the transfer buffer
1055                         grid_particles(k,j,i)%particles(n)%particle_mask = .FALSE.
1056                      ELSE
[3948]1057                         WRITE(9,'(a,10i6)') 'This should not happen ',i,j,k,ic,jc,ipl,ipr,jps,jpn
[2801]1058                         CALL MPI_Abort( MPI_COMM_WORLD, 9999, ierr )
1059                      ENDIF
1060                   ENDIF
1061                ENDIF
1062             ENDDO
1063          ENDDO
1064       ENDDO
1065    ENDDO
1066
1067!
1068!- Pack particles (eliminate those marked for deletion),
1069!- determine new number of particles
[4017]1070!    CALL lpm_sort_in_subboxes
[2801]1071
[2967]1072#endif
[2801]1073 END SUBROUTINE c_copy_particle_to_coarse_grid
[2967]1074
1075
[2801]1076!------------------------------------------------------------------------------!
1077! Description:
1078! ------------
1079!> parent routine:
1080!> copy/sort particles from the MPI window into the respective grid boxes
1081!------------------------------------------------------------------------------!
[3241]1082 SUBROUTINE p_copy_particle_to_org_grid( m )
[2801]1083
1084    IMPLICIT NONE
1085
1086    INTEGER(iwp),INTENT(IN) ::  m        !<
1087
1088    INTEGER(iwp) ::  i      !< loop index (x grid)
1089    INTEGER(iwp) ::  j      !< loop index (y grid)
1090    INTEGER(iwp) ::  k      !< loop index (z grid)
1091    INTEGER(iwp) ::  n      !< loop index (nr part)
1092    INTEGER(iwp) ::  nr     !< short variable name for nr_part
1093    INTEGER(iwp) ::  pindex !< short variable name part_adr
1094   
[3123]1095    REAL(wp) ::  x  !< x coordinate
1096    REAL(wp) ::  xo !< x origin
1097    REAL(wp) ::  y  !< y coordinate
1098    REAL(wp) ::  yo !< y origin
1099    REAL(wp) ::  z  !< z coordinate
[2801]1100   
1101    INTEGER(iwp),DIMENSION(1) ::  buf_shape !<
1102
[2967]1103#if defined( __parallel )
[2801]1104    buf_shape(1) = max_nr_particle_in_rma_win
1105    CALL C_F_POINTER( buf_ptr(m), particle_in_win , buf_shape )
1106
1107    DO i = nxl, nxr
1108       DO j = nys, nyn
1109          nr = nr_part(j,i)
1110          IF ( nr > 0 ) THEN
1111             pindex = part_adr(j,i)
1112             DO n = 1, nr
1113                x = particle_in_win(pindex)%x-lower_left_coord_x
1114                y = particle_in_win(pindex)%y-lower_left_coord_y
1115                z = particle_in_win(pindex)%z
1116                xo = particle_in_win(pindex)%origin_x-lower_left_coord_x
1117                yo = particle_in_win(pindex)%origin_y-lower_left_coord_y
1118                k = nzt + 1
[3948]1119                DO WHILE ( zw(k-1) > z .AND. k > nzb + 1 )           ! kk search loop has to be optimzed !!!
[2801]1120                   k = k - 1
1121                END DO
1122
1123                prt_count(k,j,i) = prt_count(k,j,i) + 1
1124                IF ( prt_count(k,j,i) > SIZE( grid_particles(k,j,i)%particles ) ) THEN
[4017]1125                   CALL pmc_realloc_particles_array( i, j, k )
[2801]1126                ENDIF
1127                grid_particles(k,j,i)%particles(prt_count(k,j,i)) = particle_in_win(pindex)
1128               
1129!
1130!--             Update particle positions and origins relative to parent domain
1131                grid_particles(k,j,i)%particles(prt_count(k,j,i))%x = x
1132                grid_particles(k,j,i)%particles(prt_count(k,j,i))%y = y
1133                grid_particles(k,j,i)%particles(prt_count(k,j,i))%origin_x = xo
1134                grid_particles(k,j,i)%particles(prt_count(k,j,i))%origin_y = yo
1135                pindex = pindex + 1
1136             ENDDO
1137          ENDIF
1138       ENDDO
1139    ENDDO
[2967]1140
1141#endif
[2801]1142 END SUBROUTINE p_copy_particle_to_org_grid
[4017]1143 
1144!------------------------------------------------------------------------------!
1145! Description:
1146! ------------
1147!> If the allocated memory for the particle array do not suffice to add arriving
1148!> particles from neighbour grid cells, this subrouting reallocates the
1149!> particle array to assure enough memory is available.
1150!------------------------------------------------------------------------------!
1151 SUBROUTINE pmc_realloc_particles_array ( i, j, k, size_in )
[2801]1152
[4017]1153    INTEGER(iwp), INTENT(IN)                       ::  i              !<
1154    INTEGER(iwp), INTENT(IN)                       ::  j              !<
1155    INTEGER(iwp), INTENT(IN)                       ::  k              !<
1156    INTEGER(iwp), INTENT(IN), OPTIONAL             ::  size_in        !<
[2801]1157
[4017]1158    INTEGER(iwp)                                   ::  old_size        !<
1159    INTEGER(iwp)                                   ::  new_size        !<
1160    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles_d !<
1161    TYPE(particle_type), DIMENSION(500)            ::  tmp_particles_s !<
1162
1163    old_size = SIZE(grid_particles(k,j,i)%particles)
1164
1165    IF ( PRESENT(size_in) )   THEN
1166       new_size = size_in
1167    ELSE
1168       new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
1169    ENDIF
1170
[4043]1171    new_size = MAX( new_size, 1, old_size + 1 )
[4017]1172
1173    IF ( old_size <= 500 )  THEN
1174
1175       tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
1176
1177       DEALLOCATE(grid_particles(k,j,i)%particles)
1178       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1179
1180       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_s(1:old_size)
1181       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1182
1183    ELSE
1184
1185       ALLOCATE(tmp_particles_d(new_size))
1186       tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
1187
1188       DEALLOCATE(grid_particles(k,j,i)%particles)
1189       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1190
1191       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_d(1:old_size)
1192       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1193
1194       DEALLOCATE(tmp_particles_d)
1195
1196    ENDIF
1197    particles => grid_particles(k,j,i)%particles(1:new_size)
1198
1199    RETURN
1200   
1201 END SUBROUTINE pmc_realloc_particles_array
1202
[4444]1203#endif
[2801]1204END MODULE pmc_particle_interface
Note: See TracBrowser for help on using the repository browser.