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

Last change on this file since 4043 was 4043, checked in by schwenkel, 2 years ago

further modularization of lpm and delete min_nr_particle

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