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

Last change on this file since 2801 was 2801, checked in by thiele, 3 years ago

Introduce particle transfer in nested models

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