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

Last change on this file since 4879 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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