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

Last change on this file since 4595 was 4444, checked in by raasch, 5 years ago

bugfix: cpp-directives for serial mode added

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