Ignore:
Timestamp:
Jun 9, 2016 4:25:25 PM (8 years ago)
Author:
suehring
Message:

several bugfixes in particle model and serial mode

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm_exchange_horiz.f90

    r1874 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Bugfixes:
     22! - reallocation of new particles
     23!   ( did not work for small number of min_nr_particle )
     24! - dynamical reallocation of north-south exchange arrays ( particles got lost )
     25! - north-south exchange ( nr_move_north, nr_move_south were overwritten by zero )
     26! - horizontal particle boundary conditions in serial mode
     27!
     28! Remove unused variables
     29! Descriptions in variable declaration blocks added
    2230!
    2331! Former revisions:
     
    121129    INTEGER(iwp)            ::  nr_move_south               !<
    122130
    123     TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_north
    124     TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_south
     131    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_north
     132    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_south
    125133
    126134    SAVE
     
    161169    IMPLICIT NONE
    162170
    163     INTEGER(iwp) ::  i                                       !<
    164     INTEGER(iwp) ::  ip                                      !<
    165     INTEGER(iwp) ::  j                                       !<
    166     INTEGER(iwp) ::  jp                                      !<
    167     INTEGER(iwp) ::  kp                                      !<
    168     INTEGER(iwp) ::  n                                       !<
    169     INTEGER(iwp) ::  trlp_count                              !<
    170     INTEGER(iwp) ::  trlp_count_recv                         !<
    171     INTEGER(iwp) ::  trlpt_count                             !<
    172     INTEGER(iwp) ::  trlpt_count_recv                        !<
    173     INTEGER(iwp) ::  trnp_count                              !<
    174     INTEGER(iwp) ::  trnp_count_recv                         !<
    175     INTEGER(iwp) ::  trnpt_count                             !<
    176     INTEGER(iwp) ::  trnpt_count_recv                        !<
    177     INTEGER(iwp) ::  trrp_count                              !<
    178     INTEGER(iwp) ::  trrp_count_recv                         !<
    179     INTEGER(iwp) ::  trrpt_count                             !<
    180     INTEGER(iwp) ::  trrpt_count_recv                        !<
    181     INTEGER(iwp) ::  trsp_count                              !<
    182     INTEGER(iwp) ::  trsp_count_recv                         !<
    183     INTEGER(iwp) ::  trspt_count                             !<
    184     INTEGER(iwp) ::  trspt_count_recv                        !<
    185 
    186     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !<
    187     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !<
    188     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !<
    189     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !<
    190     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !<
    191     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !<
    192     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !<
    193     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !<
     171    INTEGER(iwp) ::  i                 !< grid index (x) of particle positition
     172    INTEGER(iwp) ::  ip                !< index variable along x
     173    INTEGER(iwp) ::  j                 !< grid index (y) of particle positition
     174    INTEGER(iwp) ::  jp                !< index variable along y
     175    INTEGER(iwp) ::  kp                !< index variable along z
     176    INTEGER(iwp) ::  n                 !< particle index variable
     177    INTEGER(iwp) ::  trlp_count        !< number of particles send to left PE
     178    INTEGER(iwp) ::  trlp_count_recv   !< number of particles receive from right PE
     179    INTEGER(iwp) ::  trnp_count        !< number of particles send to north PE
     180    INTEGER(iwp) ::  trnp_count_recv   !< number of particles receive from south PE
     181    INTEGER(iwp) ::  trrp_count        !< number of particles send to right PE
     182    INTEGER(iwp) ::  trrp_count_recv   !< number of particles receive from left PE
     183    INTEGER(iwp) ::  trsp_count        !< number of particles send to south PE
     184    INTEGER(iwp) ::  trsp_count_recv   !< number of particles receive from north PE
     185
     186    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !< particles received from right PE
     187    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !< particles received from south PE
     188    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !< particles received from left PE
     189    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !< particles received from north PE
     190    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !< particles send to left PE
     191    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !< particles send to north PE
     192    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !< particles send to right PE
     193    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !< particles send to south PE
    194194
    195195    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
     
    209209!-- be adjusted.
    210210    trlp_count  = 0
    211     trlpt_count = 0
    212211    trrp_count  = 0
    213     trrpt_count = 0
    214212
    215213    trlp_count_recv   = 0
    216     trlpt_count_recv  = 0
    217214    trrp_count_recv   = 0
    218     trrpt_count_recv  = 0
    219215
    220216    IF ( pdims(1) /= 1 )  THEN
     
    232228                   IF ( particles(n)%particle_mask )  THEN
    233229                      i = ( particles(n)%x + 0.5_wp * dx ) * ddx
    234    !
    235    !--                Above calculation does not work for indices less than zero
     230!
     231!--                   Above calculation does not work for indices less than zero
    236232                      IF ( particles(n)%x < -0.5_wp * dx )  i = -1
    237233
     
    249245
    250246       IF ( trlp_count  == 0 )  trlp_count  = 1
    251        IF ( trlpt_count == 0 )  trlpt_count = 1
    252247       IF ( trrp_count  == 0 )  trrp_count  = 1
    253        IF ( trrpt_count == 0 )  trrpt_count = 1
    254248
    255249       ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
     
    259253
    260254       trlp_count  = 0
    261        trlpt_count = 0
    262255       trrp_count  = 0
    263        trrpt_count = 0
    264256
    265257    ENDIF
    266258!
    267259!-- Compute only first (nxl) and last (nxr) loop iterration
    268     DO  ip = nxl, nxr,nxr-nxl
     260    DO  ip = nxl, nxr, nxr-nxl
    269261       DO  jp = nys, nyn
    270262          DO  kp = nzb+1, nzt
     
    279271
    280272                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
    281    !
    282    !--             Above calculation does not work for indices less than zero
     273!
     274!--                Above calculation does not work for indices less than zero
    283275                   IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
    284276
    285277                   IF ( i <  nxl )  THEN
    286278                      IF ( i < 0 )  THEN
    287    !
    288    !--                   Apply boundary condition along x
     279!
     280!--                   Apply boundary condition along x
    289281                         IF ( ibc_par_lr == 0 )  THEN
    290    !
    291    !--                      Cyclic condition
     282!
     283!--                         Cyclic condition
    292284                            IF ( pdims(1) == 1 )  THEN
    293285                               particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
     
    312304
    313305                         ELSEIF ( ibc_par_lr == 1 )  THEN
    314    !
    315    !--                      Particle absorption
     306!
     307!--                         Particle absorption
    316308                            particles(n)%particle_mask = .FALSE.
    317309                            deleted_particles = deleted_particles + 1
    318310
    319311                         ELSEIF ( ibc_par_lr == 2 )  THEN
    320    !
    321    !--                      Particle reflection
     312!
     313!--                         Particle reflection
    322314                            particles(n)%x       = -particles(n)%x
    323315                            particles(n)%speed_x = -particles(n)%speed_x
     
    325317                         ENDIF
    326318                      ELSE
    327    !
    328    !--                   Store particle data in the transfer array, which will be
    329    !--                   send to the neighbouring PE
     319!
     320!--                      Store particle data in the transfer array, which will be
     321!--                      send to the neighbouring PE
    330322                         trlp_count = trlp_count + 1
    331323                         trlp(trlp_count) = particles(n)
     
    337329                   ELSEIF ( i > nxr )  THEN
    338330                      IF ( i > nx )  THEN
    339    !
    340    !--                   Apply boundary condition along x
     331!
     332!--                      Apply boundary condition along x
    341333                         IF ( ibc_par_lr == 0 )  THEN
    342    !
    343    !--                      Cyclic condition
     334!
     335!--                         Cyclic condition
    344336                            IF ( pdims(1) == 1 )  THEN
    345337                               particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
     
    358350
    359351                         ELSEIF ( ibc_par_lr == 1 )  THEN
    360    !
    361    !--                      Particle absorption
     352!
     353!--                         Particle absorption
    362354                            particles(n)%particle_mask = .FALSE.
    363355                            deleted_particles = deleted_particles + 1
    364356
    365357                         ELSEIF ( ibc_par_lr == 2 )  THEN
    366    !
    367    !--                      Particle reflection
     358!
     359!--                         Particle reflection
    368360                            particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
    369361                            particles(n)%speed_x = -particles(n)%speed_x
     
    371363                         ENDIF
    372364                      ELSE
    373    !
    374    !--                   Store particle data in the transfer array, which will be send
    375    !--                   to the neighbouring PE
     365!
     366!--                      Store particle data in the transfer array, which will be send
     367!--                      to the neighbouring PE
    376368                         trrp_count = trrp_count + 1
    377369                         trrp(trrp_count) = particles(n)
     
    383375                   ENDIF
    384376                ENDIF
     377
    385378             ENDDO
    386379          ENDDO
     
    389382
    390383!
     384!-- Allocate arrays required for north-south exchange, as these
     385!-- are used directly after particles are exchange along x-direction.
     386    ALLOCATE( move_also_north(1:NR_2_direction_move) )
     387    ALLOCATE( move_also_south(1:NR_2_direction_move) )
     388
     389    nr_move_north = 0
     390    nr_move_south = 0
     391!
    391392!-- Send left boundary, receive right boundary (but first exchange how many
    392393!-- and check, if particle storage must be extended)
     
    427428
    428429    ENDIF
    429 
    430430
    431431!
     
    435435!-- Find out first the number of particles to be transferred and allocate
    436436!-- temporary arrays needed to store them.
    437 !-- For a one-dimensional decomposition along x, no transfer is necessary,
     437!-- For a one-dimensional decomposition along y, no transfer is necessary,
    438438!-- because the particle remains on the PE.
    439439    trsp_count  = nr_move_south
    440     trspt_count = 0
    441440    trnp_count  = nr_move_north
    442     trnpt_count = 0
    443441
    444442    trsp_count_recv   = 0
    445     trspt_count_recv  = 0
    446443    trnp_count_recv   = 0
    447     trnpt_count_recv  = 0
    448444
    449445    IF ( pdims(2) /= 1 )  THEN
     
    452448!--    data
    453449       DO  ip = nxl, nxr
    454           DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
     450          DO  jp = nys, nyn, nyn-nys    !compute only first (nys) and last (nyn) loop iterration
    455451             DO  kp = nzb+1, nzt
    456452                number_of_particles = prt_count(kp,jp,ip)
     
    476472
    477473       IF ( trsp_count  == 0 )  trsp_count  = 1
    478        IF ( trspt_count == 0 )  trspt_count = 1
    479474       IF ( trnp_count  == 0 )  trnp_count  = 1
    480        IF ( trnpt_count == 0 )  trnpt_count = 1
    481475
    482476       ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
     
    486480
    487481       trsp_count  = nr_move_south
    488        trspt_count = 0
    489482       trnp_count  = nr_move_north
    490        trnpt_count = 0
    491 
     483       
    492484       trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
    493485       trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
     
    506498!--             be moved.
    507499                IF ( particles(n)%particle_mask )  THEN
     500
    508501                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
    509502!
     
    521514                               particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
    522515                               particles(n)%origin_y = ( ny + 1 ) * dy + &
    523                                particles(n)%origin_y
     516                                                     particles(n)%origin_y
    524517                            ELSE
    525                                trsp_count = trsp_count + 1
    526                                trsp(trsp_count) = particles(n)
     518                               trsp_count         = trsp_count + 1
     519                               trsp(trsp_count)   = particles(n)
    527520                               trsp(trsp_count)%y = ( ny + 1 ) * dy + &
    528                                trsp(trsp_count)%y
     521                                                 trsp(trsp_count)%y
    529522                               trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
    530                                + ( ny + 1 ) * dy
     523                                                + ( ny + 1 ) * dy
    531524                               particles(n)%particle_mask = .FALSE.
    532525                               deleted_particles = deleted_particles + 1
     
    536529                                  !++ why is 1 subtracted in next statement???
    537530                                  trsp(trsp_count)%origin_y =                        &
    538                                   trsp(trsp_count)%origin_y - 1
     531                                                  trsp(trsp_count)%origin_y - 1
    539532                               ENDIF
    540533
     
    545538!--                         Particle absorption
    546539                            particles(n)%particle_mask = .FALSE.
    547                             deleted_particles = deleted_particles + 1
     540                            deleted_particles          = deleted_particles + 1
    548541
    549542                         ELSEIF ( ibc_par_ns == 2 )  THEN
     
    568561                      IF ( j > ny )  THEN
    569562!
    570 !--                       Apply boundary condition along x
     563!--                       Apply boundary condition along y
    571564                         IF ( ibc_par_ns == 0 )  THEN
    572565!
    573566!--                         Cyclic condition
    574567                            IF ( pdims(2) == 1 )  THEN
    575                                particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
     568                               particles(n)%y        = particles(n)%y - ( ny + 1 ) * dy
    576569                               particles(n)%origin_y =                         &
    577                                   particles(n)%origin_y - ( ny + 1 ) * dy
     570                                          particles(n)%origin_y - ( ny + 1 ) * dy
    578571                            ELSE
    579                                trnp_count = trnp_count + 1
    580                                trnp(trnp_count) = particles(n)
     572                               trnp_count         = trnp_count + 1
     573                               trnp(trnp_count)   = particles(n)
    581574                               trnp(trnp_count)%y =                            &
    582                                   trnp(trnp_count)%y - ( ny + 1 ) * dy
     575                                          trnp(trnp_count)%y - ( ny + 1 ) * dy
    583576                               trnp(trnp_count)%origin_y =                     &
    584                                   trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
     577                                         trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
    585578                               particles(n)%particle_mask = .FALSE.
    586                                deleted_particles = deleted_particles + 1
     579                               deleted_particles          = deleted_particles + 1
    587580                            ENDIF
    588581
     
    628621
    629622       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
    630 
     623 
    631624       CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
    632625                          psouth, 1, rvnp(1)%radius,                             &
     
    661654    ENDIF
    662655
     656    DEALLOCATE( move_also_north )
     657    DEALLOCATE( move_also_south )
     658
    663659#else
    664660
    665 !
    666 !-- Apply boundary conditions
    667     DO  n = 1, number_of_particles
    668 
    669        IF ( particles(n)%x < -0.5_wp * dx )  THEN
    670 
    671           IF ( ibc_par_lr == 0 )  THEN
    672 !
    673 !--          Cyclic boundary. Relevant coordinate has to be changed.
    674              particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
    675 
    676           ELSEIF ( ibc_par_lr == 1 )  THEN
    677 !
    678 !--          Particle absorption
    679              particles(n)%particle_mask = .FALSE.
    680              deleted_particles = deleted_particles + 1
    681 
    682           ELSEIF ( ibc_par_lr == 2 )  THEN
    683 !
    684 !--          Particle reflection
    685              particles(n)%x       = -dx - particles(n)%x
    686              particles(n)%speed_x = -particles(n)%speed_x
    687           ENDIF
    688 
    689        ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
    690 
    691           IF ( ibc_par_lr == 0 )  THEN
    692 !
    693 !--          Cyclic boundary. Relevant coordinate has to be changed.
    694              particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
    695 
    696           ELSEIF ( ibc_par_lr == 1 )  THEN
    697 !
    698 !--          Particle absorption
    699              particles(n)%particle_mask = .FALSE.
    700              deleted_particles = deleted_particles + 1
    701 
    702           ELSEIF ( ibc_par_lr == 2 )  THEN
    703 !
    704 !--          Particle reflection
    705              particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
    706              particles(n)%speed_x = -particles(n)%speed_x
    707           ENDIF
    708 
    709        ENDIF
    710 
    711        IF ( particles(n)%y < -0.5_wp * dy )  THEN
    712 
    713           IF ( ibc_par_ns == 0 )  THEN
    714 !
    715 !--          Cyclic boundary. Relevant coordinate has to be changed.
    716              particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
    717 
    718           ELSEIF ( ibc_par_ns == 1 )  THEN
    719 !
    720 !--          Particle absorption
    721              particles(n)%particle_mask = .FALSE.
    722              deleted_particles = deleted_particles + 1
    723 
    724           ELSEIF ( ibc_par_ns == 2 )  THEN
    725 !
    726 !--          Particle reflection
    727              particles(n)%y       = -dy - particles(n)%y
    728              particles(n)%speed_y = -particles(n)%speed_y
    729           ENDIF
    730 
    731        ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
    732 
    733           IF ( ibc_par_ns == 0 )  THEN
    734 !
    735 !--          Cyclic boundary. Relevant coordinate has to be changed.
    736              particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
    737 
    738           ELSEIF ( ibc_par_ns == 1 )  THEN
    739 !
    740 !--          Particle absorption
    741              particles(n)%particle_mask = .FALSE.
    742              deleted_particles = deleted_particles + 1
    743 
    744           ELSEIF ( ibc_par_ns == 2 )  THEN
    745 !
    746 !--          Particle reflection
    747              particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
    748              particles(n)%speed_y = -particles(n)%speed_y
    749           ENDIF
    750 
    751        ENDIF
     661    DO  ip = nxl, nxr, nxr-nxl
     662       DO  jp = nys, nyn
     663          DO  kp = nzb+1, nzt
     664             number_of_particles = prt_count(kp,jp,ip)
     665             IF ( number_of_particles <= 0 )  CYCLE
     666             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     667             DO  n = 1, number_of_particles
     668!
     669!--             Apply boundary conditions
     670
     671                IF ( particles(n)%x < -0.5_wp * dx )  THEN
     672
     673                   IF ( ibc_par_lr == 0 )  THEN
     674!
     675!--                   Cyclic boundary. Relevant coordinate has to be changed.
     676                      particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
     677
     678                   ELSEIF ( ibc_par_lr == 1 )  THEN
     679!
     680!--                   Particle absorption
     681                      particles(n)%particle_mask = .FALSE.
     682                      deleted_particles = deleted_particles + 1
     683
     684                   ELSEIF ( ibc_par_lr == 2 )  THEN
     685!
     686!--                   Particle reflection
     687                      particles(n)%x       = -dx - particles(n)%x
     688                      particles(n)%speed_x = -particles(n)%speed_x
     689                   ENDIF
     690
     691                ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
     692
     693                   IF ( ibc_par_lr == 0 )  THEN
     694!
     695!--                   Cyclic boundary. Relevant coordinate has to be changed.
     696                      particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
     697
     698                   ELSEIF ( ibc_par_lr == 1 )  THEN
     699!
     700!--                   Particle absorption
     701                      particles(n)%particle_mask = .FALSE.
     702                      deleted_particles = deleted_particles + 1
     703
     704                   ELSEIF ( ibc_par_lr == 2 )  THEN
     705!
     706!--                   Particle reflection
     707                      particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
     708                      particles(n)%speed_x = -particles(n)%speed_x
     709                   ENDIF
     710
     711                ENDIF
     712             ENDDO
     713          ENDDO
     714       ENDDO
    752715    ENDDO
    753716
     717    DO  ip = nxl, nxr
     718       DO  jp = nys, nyn, nyn-nys
     719          DO  kp = nzb+1, nzt
     720             number_of_particles = prt_count(kp,jp,ip)
     721             IF ( number_of_particles <= 0 )  CYCLE
     722             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     723             DO  n = 1, number_of_particles
     724
     725                IF ( particles(n)%y < -0.5_wp * dy )  THEN
     726
     727                   IF ( ibc_par_ns == 0 )  THEN
     728!
     729!--                   Cyclic boundary. Relevant coordinate has to be changed.
     730                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
     731
     732                   ELSEIF ( ibc_par_ns == 1 )  THEN
     733!
     734!--                   Particle absorption
     735                      particles(n)%particle_mask = .FALSE.
     736                      deleted_particles = deleted_particles + 1
     737
     738                   ELSEIF ( ibc_par_ns == 2 )  THEN
     739!
     740!--                   Particle reflection
     741                      particles(n)%y       = -dy - particles(n)%y
     742                      particles(n)%speed_y = -particles(n)%speed_y
     743                   ENDIF
     744
     745                ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
     746
     747                   IF ( ibc_par_ns == 0 )  THEN
     748!
     749!--                   Cyclic boundary. Relevant coordinate has to be changed.
     750                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
     751
     752                   ELSEIF ( ibc_par_ns == 1 )  THEN
     753!
     754!--                   Particle absorption
     755                      particles(n)%particle_mask = .FALSE.
     756                      deleted_particles = deleted_particles + 1
     757
     758                   ELSEIF ( ibc_par_ns == 2 )  THEN
     759!
     760!--                   Particle reflection
     761                      particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
     762                      particles(n)%speed_y = -particles(n)%speed_y
     763                   ENDIF
     764
     765                ENDIF
     766
     767             ENDDO
     768          ENDDO
     769       ENDDO
     770    ENDDO
    754771#endif
    755772
     
    782799    IMPLICIT NONE
    783800
    784     INTEGER(iwp)        ::  ip        !<
    785     INTEGER(iwp)        ::  jp        !<
    786     INTEGER(iwp)        ::  kp        !<
    787     INTEGER(iwp)        ::  n         !<
    788     INTEGER(iwp)        ::  pindex    !<
     801    INTEGER(iwp)        ::  ip        !< grid index (x) of particle
     802    INTEGER(iwp)        ::  jp        !< grid index (x) of particle
     803    INTEGER(iwp)        ::  kp        !< grid index (x) of particle
     804    INTEGER(iwp)        ::  n         !< index variable of particle
     805    INTEGER(iwp)        ::  pindex    !< dummy argument for new number of particles per grid box
    789806
    790807    LOGICAL             ::  pack_done !<
    791808
    792     TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
     809    TYPE(particle_type), DIMENSION(:), INTENT(IN)  ::  particle_array !< new particles in a grid box
     810    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  temp_ns        !< temporary particle array for reallocation
    793811
    794812    pack_done     = .FALSE.
    795813
    796     nr_move_north = 0
    797     nr_move_south = 0
    798 
    799814    DO n = 1, SIZE(particle_array)
     815
     816       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
     817
    800818       ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
    801819       jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
    802        kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
     820       kp =   particle_array(n)%z / dz + 1 + offset_ocean_nzt
    803821
    804822       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
     
    824842          prt_count(kp,jp,ip) = pindex
    825843       ELSE
    826           IF ( jp == nys - 1 )  THEN
     844          IF ( jp <= nys - 1 )  THEN
    827845             nr_move_south = nr_move_south+1
     846!
     847!--          Before particle information is swapped to exchange-array, check
     848!--          if enough memory is allocated. If required, reallocate exchange
     849!--          array.
     850             IF ( nr_move_south > SIZE(move_also_south) )  THEN
     851!
     852!--             At first, allocate further temporary array to swap particle
     853!--             information.
     854                ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move) )
     855                temp_ns(1:nr_move_south-1) = move_also_south(1:nr_move_south-1)
     856                DEALLOCATE( move_also_south )
     857                ALLOCATE( move_also_south(SIZE(temp_ns)) )
     858                move_also_south(1:nr_move_south-1) = temp_ns(1:nr_move_south-1)
     859                DEALLOCATE( temp_ns )
     860
     861             ENDIF
     862
    828863             move_also_south(nr_move_south) = particle_array(n)
     864
    829865             IF ( jp == -1 )  THEN
    830866                move_also_south(nr_move_south)%y =                             &
     
    833869                   move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
    834870             ENDIF
    835           ELSEIF ( jp == nyn+1 )  THEN
     871          ELSEIF ( jp >= nyn+1 )  THEN
    836872             nr_move_north = nr_move_north+1
     873!
     874!--          Before particle information is swapped to exchange-array, check
     875!--          if enough memory is allocated. If required, reallocate exchange
     876!--          array.
     877             IF ( nr_move_north > SIZE(move_also_north) )  THEN
     878!
     879!--             At first, allocate further temporary array to swap particle
     880!--             information.
     881                ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move) )
     882                temp_ns(1:nr_move_north-1) = move_also_south(1:nr_move_north-1)
     883                DEALLOCATE( move_also_north )
     884                ALLOCATE( move_also_north(SIZE(temp_ns)) )
     885                move_also_north(1:nr_move_north-1) = temp_ns(1:nr_move_north-1)
     886                DEALLOCATE( temp_ns )
     887
     888             ENDIF
     889
    837890             move_also_north(nr_move_north) = particle_array(n)
    838891             IF ( jp == ny+1 )  THEN
     
    867920    IMPLICIT NONE
    868921
    869     INTEGER(iwp)        ::  i           !<
    870     INTEGER(iwp)        ::  ip          !<
    871     INTEGER(iwp)        ::  j           !<
    872     INTEGER(iwp)        ::  jp          !<
    873     INTEGER(iwp)        ::  k           !<
    874     INTEGER(iwp)        ::  kp          !<
    875     INTEGER(iwp)        ::  n           !<
    876     INTEGER(iwp)        ::  np_old_cell !<
    877     INTEGER(iwp)        ::  n_start     !<
    878     INTEGER(iwp)        ::  pindex      !<
     922    INTEGER(iwp)        ::  i           !< grid index (x) of particle position
     923    INTEGER(iwp)        ::  ip          !< index variable along x
     924    INTEGER(iwp)        ::  j           !< grid index (y) of particle position
     925    INTEGER(iwp)        ::  jp          !< index variable along y
     926    INTEGER(iwp)        ::  k           !< grid index (z) of particle position
     927    INTEGER(iwp)        ::  kp          !< index variable along z
     928    INTEGER(iwp)        ::  n           !< index variable for particle array
     929    INTEGER(iwp)        ::  np_old_cell !< number of particles per grid box before moving
     930    INTEGER(iwp)        ::  n_start     !< start index
     931    INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
    879932
    880933    LOGICAL             ::  pack_done   !<
    881934
    882     TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !<
     935    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !< particles before moving
    883936
    884937    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
     
    9471000                            IF ( pack_done )  THEN
    9481001                               CALL realloc_particles_array(i,j,k)
    949                                pindex = prt_count(k,j,i)+1
    9501002                            ELSE
    9511003                               CALL lpm_pack_arrays
    9521004                               prt_count(k,j,i) = number_of_particles
    953 
    954                                pindex = prt_count(k,j,i)+1
    9551005!
    9561006!--                            If number of particles in the new grid box
     
    9591009                               IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
    9601010                                  CALL realloc_particles_array(i,j,k)
    961                                   pindex = prt_count(k,j,i)+1
    9621011                               ENDIF
    9631012
     
    9871036! Description:
    9881037! ------------
    989 !> @todo Missing subroutine description.
     1038!> If the allocated memory for the particle array do not suffice to add arriving
     1039!> particles from neighbour grid cells, this subrouting reallocates the
     1040!> particle array to assure enough memory is available.
    9901041!------------------------------------------------------------------------------!
    9911042 SUBROUTINE realloc_particles_array (i,j,k,size_in)
     
    10121063    ENDIF
    10131064
    1014     new_size = MAX( new_size, min_nr_particle )
     1065    new_size = MAX( new_size, min_nr_particle, old_size + 1 )
    10151066
    10161067    IF ( old_size <= 500 )  THEN
Note: See TracChangeset for help on using the changeset viewer.