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

Last change on this file since 4883 was 4883, checked in by hellstea, 3 years ago

user switch for particle coupling added

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