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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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