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

Last change on this file since 4017 was 4017, checked in by schwenkel, 5 years ago

Modularization of all lagrangian particle model code components

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