Changeset 2606


Ignore:
Timestamp:
Nov 10, 2017 10:36:31 AM (6 years ago)
Author:
schwenkel
Message:

Modified particle box location and further changes in particle model

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

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

    r2418 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2418 2017-09-06 15:24:24Z suehring
    2734! Major bugfixes in modeling SGS particle speeds (since revision 1359).
    2835! Particle sorting added to distinguish between already completed and
     
    144151        ONLY: lpm_create_particle, PHASE_RELEASE
    145152
    146     USE lpm_pack_arrays_mod,                                                   &
    147         ONLY:  lpm_pack_all_arrays, lpm_sort
     153    USE lpm_pack_and_sort_mod,                                                   &
     154        ONLY:  lpm_sort_in_subboxes, lpm_sort_timeloop_done
    148155
    149156    USE particle_attributes,                                                   &
     
    281288!--    according to their sub-boxes.
    282289       IF ( .NOT. first_loop_stride  .AND.  use_sgs_for_particles )            &
    283           CALL lpm_sort
     290          CALL lpm_sort_timeloop_done
    284291
    285292       DO  i = nxl, nxr
     
    408415!--    Pack particles (eliminate those marked for deletion),
    409416!--    determine new number of particles
    410        CALL lpm_pack_all_arrays
     417       CALL lpm_sort_in_subboxes
    411418!
    412419!--    Initialize variables for the next (sub-) timestep, i.e., for marking
  • palm/trunk/SOURCE/lpm_advec.f90

    r2417 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2417 2017-09-06 15:22:27Z suehring
    2734! Particle loops adapted for sub-box structure, i.e. for each sub-box the
    2835! particle loop runs from start_index up to end_index instead from 1 to
     
    273280
    274281    DO  nb = 0, 7
    275        
     282!
     283!--    Interpolate u velocity-component       
    276284       i = ip
    277285       j = jp + block_offset(nb)%j_off
    278286       k = kp + block_offset(nb)%k_off
    279 !
    280 !--    Interpolate u velocity-component
     287
    281288       DO  n = start_index(nb), end_index(nb)
    282289!
     
    289296!--       First, check if particle is located below first vertical grid level
    290297!--       above topography (Prandtl-layer height)
    291           ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
    292           jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
     298          ilog = particles(n)%x * ddx
     299          jlog = particles(n)%y * ddy
    293300!
    294301!--       Determine vertical index of topography top
     
    378385
    379386       ENDDO
    380 
     387!
     388!--    Same procedure for interpolation of the v velocity-component
    381389       i = ip + block_offset(nb)%i_off
    382390       j = jp
    383391       k = kp + block_offset(nb)%k_off
    384 !
    385 !--    Same procedure for interpolation of the v velocity-component
     392
    386393       DO  n = start_index(nb), end_index(nb)
    387394
    388           ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
    389           jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
     395          ilog = particles(n)%x * ddx
     396          jlog = particles(n)%y * ddy
    390397!
    391398!--       Determine vertical index of topography top
     
    474481
    475482       ENDDO
    476 
     483!
     484!--    Same procedure for interpolation of the w velocity-component
    477485       i = ip + block_offset(nb)%i_off
    478486       j = jp + block_offset(nb)%j_off
    479487       k = kp - 1
    480 !
    481 !--    Same procedure for interpolation of the w velocity-component
     488
    482489       DO  n = start_index(nb), end_index(nb)
    483490
     
    660667                i = particles(n)%x * ddx
    661668                j = particles(n)%y * ddy
    662                 k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
     669                k = ( zv(n) + dz * atmos_ocean_sign ) / dz  &
    663670                    + offset_ocean_nzt                      ! only exact if eq.dist
    664671!
  • palm/trunk/SOURCE/lpm_boundary_conds.f90

    r2318 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2318 2017-07-20 17:27:44Z suehring
    2734! Get topography top index via Function call
    2835!
     
    254261!
    255262!--       Obtain x/y indices for current particle position
    256           i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
    257           j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
     263          i2 = particles(n)%x * ddx
     264          j2 = particles(n)%y * ddy
    258265!
    259266!--       Save current particle positions
     
    268275!
    269276!--       Obtain x/y indices for old particle positions
    270           i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
    271           j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
     277          i1 = pos_x_old * ddx
     278          j1 = pos_y_old * ddy
    272279!
    273280!--       Determine horizontal as well as vertical walls at which particle can
  • palm/trunk/SOURCE/lpm_exchange_horiz.f90

    r2330 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! lpm_move_particle now moves inner particles that would leave the core
     30! to the respective boundary gridbox. Afterwards they get transferred by
     31! lpm_exchange_horiz. This allows particles to move more than one gridbox
     32! independent of domain composition. 
     33! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     34! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     35! lpm_sort -> lpm_sort_timeloop_done
     36!
     37! 2330 2017-08-03 14:26:02Z schwenkel
    2738! Bugfix: Also for gfortran compilable, function c_sizeof is not used.
    2839!
     
    128139    USE kinds
    129140
    130     USE lpm_pack_arrays_mod,                                                   &
    131         ONLY:  lpm_pack_arrays
     141    USE lpm_pack_and_sort_mod,                                                   &
     142        ONLY:  lpm_pack
    132143
    133144    USE netcdf_interface,                                                      &
     
    255266                DO  n = 1, number_of_particles
    256267                   IF ( particles(n)%particle_mask )  THEN
    257                       i = ( particles(n)%x + 0.5_wp * dx ) * ddx
     268                      i = particles(n)%x * ddx
    258269!
    259270!--                   Above calculation does not work for indices less than zero
    260                       IF ( particles(n)%x < -0.5_wp * dx )  i = -1
     271                      IF ( particles(n)%x < 0.0_wp)  i = -1
    261272
    262273                      IF ( i < nxl )  THEN
     
    298309                IF ( particles(n)%particle_mask )  THEN
    299310
    300                    i = ( particles(n)%x + 0.5_wp * dx ) * ddx
     311                   i = particles(n)%x * ddx
    301312!
    302313!--                Above calculation does not work for indices less than zero
    303                    IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
     314                   IF ( particles(n)%x < 0.0_wp )  i = -1
    304315
    305316                   IF ( i <  nxl )  THEN
     
    323334                               deleted_particles = deleted_particles + 1
    324335
    325                                IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp )  THEN
     336                               IF ( trlp(trlp_count)%x >= (nx + 1)* dx - 1.0E-12_wp )  THEN
    326337                                  trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
    327338                                  !++ why is 1 subtracted in next statement???
     
    501512                DO  n = 1, number_of_particles
    502513                   IF ( particles(n)%particle_mask )  THEN
    503                       j = ( particles(n)%y + 0.5_wp * dy ) * ddy
     514                      j = particles(n)%y * ddy
    504515!
    505516!--                   Above calculation does not work for indices less than zero
    506                       IF ( particles(n)%y < -0.5_wp * dy )  j = -1
     517                      IF ( particles(n)%y < 0.0_wp)  j = -1
    507518
    508519                      IF ( j < nys )  THEN
     
    545556                IF ( particles(n)%particle_mask )  THEN
    546557
    547                    j = ( particles(n)%y + 0.5_wp * dy ) * ddy
     558                   j = particles(n)%y * ddy
    548559!
    549560!--                Above calculation does not work for indices less than zero
    550                    IF ( particles(n)%y < -0.5_wp * dy )  j = -1
     561                   IF ( particles(n)%y < 0.0_wp )  j = -1
    551562
    552563                   IF ( j < nys )  THEN
     
    571582                               deleted_particles = deleted_particles + 1
    572583
    573                                IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp )  THEN
     584                               IF ( trsp(trsp_count)%y >= (ny+1)* dy - 1.0E-12_wp )  THEN
    574585                                  trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
    575586                                  !++ why is 1 subtracted in next statement???
     
    719730!--             Apply boundary conditions
    720731
    721                 IF ( particles(n)%x < -0.5_wp * dx )  THEN
     732                IF ( particles(n)%x < 0.0_wp )  THEN
    722733
    723734                   IF ( ibc_par_lr == 0 )  THEN
     
    725736!--                   Cyclic boundary. Relevant coordinate has to be changed.
    726737                      particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
    727 
     738                      particles(n)%origin_x = ( nx + 1 ) * dx + &
     739                               particles(n)%origin_x
    728740                   ELSEIF ( ibc_par_lr == 1 )  THEN
    729741!
     
    739751                   ENDIF
    740752
    741                 ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
     753                ELSEIF ( particles(n)%x >= ( nx + 1) * dx )  THEN
    742754
    743755                   IF ( ibc_par_lr == 0 )  THEN
     
    745757!--                   Cyclic boundary. Relevant coordinate has to be changed.
    746758                      particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
     759                      particles(n)%origin_x = particles(n)%origin_x - &
     760                               ( nx + 1 ) * dx
    747761
    748762                   ELSEIF ( ibc_par_lr == 1 )  THEN
     
    773787             DO  n = 1, number_of_particles
    774788
    775                 IF ( particles(n)%y < -0.5_wp * dy )  THEN
     789                IF ( particles(n)%y < 0.0_wp)  THEN
    776790
    777791                   IF ( ibc_par_ns == 0 )  THEN
     
    779793!--                   Cyclic boundary. Relevant coordinate has to be changed.
    780794                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
     795                      particles(n)%origin_y = ( ny + 1 ) * dy + &
     796                           particles(n)%origin_y
    781797
    782798                   ELSEIF ( ibc_par_ns == 1 )  THEN
     
    793809                   ENDIF
    794810
    795                 ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
     811                ELSEIF ( particles(n)%y >= ( ny + 1) * dy )  THEN
    796812
    797813                   IF ( ibc_par_ns == 0 )  THEN
     
    799815!--                   Cyclic boundary. Relevant coordinate has to be changed.
    800816                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
     817                      particles(n)%origin_y = particles(n)%origin_y - &
     818                                ( ny + 1 ) * dy
    801819
    802820                   ELSEIF ( ibc_par_ns == 1 )  THEN
     
    866884       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
    867885
    868        ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
    869        jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
    870        kp =   particle_array(n)%z / dz + 1 + offset_ocean_nzt
     886       ip = particle_array(n)%x * ddx
     887       jp = particle_array(n)%y * ddy
     888       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
    871889
    872890       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
     
    880898                CALL realloc_particles_array (ip,jp,kp)
    881899             ELSE
    882                 CALL lpm_pack_arrays
     900                CALL lpm_pack
    883901                prt_count(kp,jp,ip) = number_of_particles
    884902                pindex = prt_count(kp,jp,ip)+1
     
    10141032    INTEGER(iwp)        ::  kp          !< index variable along z
    10151033    INTEGER(iwp)        ::  n           !< index variable for particle array
    1016     INTEGER(iwp)        ::  np_old_cell !< number of particles per grid box before moving
    1017     INTEGER(iwp)        ::  n_start     !< start index
     1034    INTEGER(iwp)        ::  np_before_move !< number of particles per grid box before moving
    10181035    INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
    10191036
    1020     LOGICAL             ::  pack_done   !<
    1021 
    1022     TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !< particles before moving
     1037    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_before_move !< particles before moving
    10231038
    10241039    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
     
    10281043          DO  kp = nzb+1, nzt
    10291044
    1030              np_old_cell = prt_count(kp,jp,ip)
    1031              IF ( np_old_cell <= 0 )  CYCLE
    1032              particles_old_cell => grid_particles(kp,jp,ip)%particles(1:np_old_cell)
    1033              n_start = -1
     1045             np_before_move = prt_count(kp,jp,ip)
     1046             IF ( np_before_move <= 0 )  CYCLE
     1047             particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move)
    10341048             
    1035              DO  n = 1, np_old_cell
    1036                 i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
    1037                 j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
    1038                 k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
     1049             DO  n = 1, np_before_move
     1050                i = particles_before_move(n)%x * ddx
     1051                j = particles_before_move(n)%y * ddy
     1052                k = particles_before_move(n)%z / dz + 1 + offset_ocean_nzt
     1053
     1054!--             For lpm_exchange_horiz to work properly particles need to be moved to the outermost gridboxes
     1055!--             of the respective processor. If the particle index is inside the processor the following lines
     1056!--             will not change the index
     1057                i = MIN ( i , nxr )
     1058                i = MAX ( i , nxl )
     1059                j = MIN ( j , nyn )
     1060                j = MAX ( j , nys )
     1061                k = MIN ( k , nzt )
     1062                k = MAX ( k , nzb+1 )
    10391063!
    10401064!--             Check, if particle has moved to another grid cell.
    10411065                IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
    1042 !
    1043 !--                The particle has moved to another grid cell. Now check, if
    1044 !--                particle stays on the same processor.
    1045                    IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.      &
    1046                         j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
    1047 !
    1048 !--                   If the particle stays on the same processor, the particle
    1049 !--                   will be added to the particle array of the new processor.
    1050                       number_of_particles = prt_count(k,j,i)
    1051                       particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    1052 
    1053                       pindex = prt_count(k,j,i)+1
    1054                       IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )  &
    1055                       THEN
    1056                          n_start = n
    1057                          EXIT
    1058                       ENDIF
    1059 
    1060                       grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
    1061                       prt_count(k,j,i) = pindex
    1062 
    1063                       particles_old_cell(n)%particle_mask = .FALSE.
     1066!!
     1067!--                If the particle stays on the same processor, the particle
     1068!--                will be added to the particle array of the new processor.
     1069                   number_of_particles = prt_count(k,j,i)
     1070                   particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     1071
     1072                   pindex = prt_count(k,j,i)+1
     1073                   IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )     &
     1074                   THEN
     1075                      CALL realloc_particles_array(i,j,k)
    10641076                   ENDIF
     1077
     1078                   grid_particles(k,j,i)%particles(pindex) = particles_before_move(n)
     1079                   prt_count(k,j,i) = pindex
     1080
     1081                   particles_before_move(n)%particle_mask = .FALSE.
    10651082                ENDIF
    10661083             ENDDO
    10671084
    1068              IF ( n_start >= 0 )  THEN
    1069                 pack_done = .FALSE.
    1070                 DO  n = n_start, np_old_cell
    1071                    i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
    1072                    j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
    1073                    k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
    1074                    IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
    1075 !
    1076 !--                   Particle is in different box
    1077                       IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.   &
    1078                            j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
    1079 !
    1080 !--                      Particle stays on processor
    1081                          number_of_particles = prt_count(k,j,i)
    1082                          particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    1083 
    1084                          pindex = prt_count(k,j,i)+1
    1085                          IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
    1086                          THEN
    1087                             IF ( pack_done )  THEN
    1088                                CALL realloc_particles_array(i,j,k)
    1089                             ELSE
    1090                                CALL lpm_pack_arrays
    1091                                prt_count(k,j,i) = number_of_particles
    1092 !
    1093 !--                            Set pindex to next index of the modified
    1094 !--                            particle-array (due to lpm_pack_arrays)
    1095                                pindex = prt_count(k,j,i)+1
    1096 !
    1097 !--                            If number of particles in the new grid box
    1098 !--                            exceeds its allocated memory, the particle array
    1099 !--                            will be reallocated
    1100                                IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
    1101                                   CALL realloc_particles_array(i,j,k)
    1102                                ENDIF
    1103 
    1104                                pack_done = .TRUE.
    1105                             ENDIF
    1106                          ENDIF
    1107 
    1108                          grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
    1109                          prt_count(k,j,i) = pindex
    1110 
    1111                          particles_old_cell(n)%particle_mask = .FALSE.
    1112                       ENDIF
    1113                    ENDIF
    1114                 ENDDO
    1115              ENDIF
    11161085          ENDDO
    11171086       ENDDO
  • palm/trunk/SOURCE/lpm_init.f90

    r2375 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2375 2017-08-29 14:10:28Z schwenkel
    2734! Initialization of chemical aerosol composition
    2835!
     
    277284!
    278285!-- Define block offsets for dividing a gridcell in 8 sub cells
    279 
    280     block_offset(0) = block_offset_def (-1,-1,-1)
    281     block_offset(1) = block_offset_def (-1,-1, 0)
    282     block_offset(2) = block_offset_def (-1, 0,-1)
    283     block_offset(3) = block_offset_def (-1, 0, 0)
    284     block_offset(4) = block_offset_def ( 0,-1,-1)
    285     block_offset(5) = block_offset_def ( 0,-1, 0)
    286     block_offset(6) = block_offset_def ( 0, 0,-1)
    287     block_offset(7) = block_offset_def ( 0, 0, 0)
     286!-- See documentation for List of subgrid boxes
     287!-- See pack_and_sort in lpm_pack_arrays.f90 for assignment of the subgrid boxes
     288    block_offset(0) = block_offset_def ( 0, 0, 0)
     289    block_offset(1) = block_offset_def ( 0, 0,-1)
     290    block_offset(2) = block_offset_def ( 0,-1, 0)
     291    block_offset(3) = block_offset_def ( 0,-1,-1)
     292    block_offset(4) = block_offset_def (-1, 0, 0)
     293    block_offset(5) = block_offset_def (-1, 0,-1)
     294    block_offset(6) = block_offset_def (-1,-1, 0)
     295    block_offset(7) = block_offset_def (-1,-1,-1)
    288296!
    289297!-- Check the number of particle groups.
     
    309317!
    310318!-- Set default start positions, if necessary
    311     IF ( psl(1) == 9999999.9_wp )  psl(1) = -0.5_wp * dx
    312     IF ( psr(1) == 9999999.9_wp )  psr(1) = ( nx + 0.5_wp ) * dx
    313     IF ( pss(1) == 9999999.9_wp )  pss(1) = -0.5_wp * dy
    314     IF ( psn(1) == 9999999.9_wp )  psn(1) = ( ny + 0.5_wp ) * dy
     319    IF ( psl(1) == 9999999.9_wp )  psl(1) = 0.0_wp
     320    IF ( psr(1) == 9999999.9_wp )  psr(1) = ( nx +1 ) * dx
     321    IF ( pss(1) == 9999999.9_wp )  pss(1) = 0.0_wp
     322    IF ( psn(1) == 9999999.9_wp )  psn(1) = ( ny +1 ) * dy
    315323    IF ( psb(1) == 9999999.9_wp )  psb(1) = zu(nz/2)
    316324    IF ( pst(1) == 9999999.9_wp )  pst(1) = psb(1)
     
    617625        ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
    618626
    619     USE lpm_pack_arrays_mod,                                                   &
    620         ONLY: lpm_pack_all_arrays
     627    USE lpm_pack_and_sort_mod,                                                   &
     628        ONLY: lpm_sort_in_subboxes
    621629
    622630    USE particle_attributes,                                                   &
     
    683691                DO WHILE ( pos_y <= psn(i) )
    684692
    685                    IF ( pos_y >= ( nys - 0.5_wp ) * dy  .AND.                  &
    686                         pos_y <  ( nyn + 0.5_wp ) * dy ) THEN
     693                   IF ( pos_y >= nys * dy  .AND.                  &
     694                        pos_y <  ( nyn + 1 ) * dy  ) THEN
    687695
    688696                      pos_x = psl(i)
     
    690698               xloop: DO WHILE ( pos_x <= psr(i) )
    691699
    692                          IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            &
    693                               pos_x <  ( nxr + 0.5_wp ) * dx ) THEN
     700                         IF ( pos_x >= nxl * dx  .AND.            &
     701                              pos_x <  ( nxr + 1) * dx ) THEN
    694702
    695703                            DO  j = 1, particles_per_point
     
    729737!
    730738!--                            Determine the grid indices of the particle position
    731                                ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
    732                                jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
     739                               ip = tmp_particle%x * ddx
     740                               jp = tmp_particle%y * ddy
    733741                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
    734742!
     
    916924                       grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    917925                DO  n = local_start(kp,jp,ip), number_of_particles
    918                    i = ( particles(n)%x + 0.5_wp * dx ) * ddx
    919                    j = ( particles(n)%y + 0.5_wp * dy ) * ddy
     926                   i = particles(n)%x * ddx
     927                   j = particles(n)%y * ddy
    920928                   k =   particles(n)%z / dz + 1 + offset_ocean_nzt
    921929!
     
    941949!-- which is needed for a fast interpolation of the LES fields on the particle
    942950!-- position.
    943     CALL lpm_pack_all_arrays
     951    CALL lpm_sort_in_subboxes
    944952
    945953!
  • palm/trunk/SOURCE/lpm_pack_arrays.f90

    r2417 r2606  
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2417 2017-09-06 15:22:27Z suehring
    2734! New routine which sorts particles into particles that completed and not
    2835! completed the LES timestep.
     
    6875!> its timestep.
    6976!------------------------------------------------------------------------------!
    70  MODULE lpm_pack_arrays_mod
     77 MODULE lpm_pack_and_sort_mod
    7178 
    7279
     
    7683
    7784    PRIVATE
    78     PUBLIC lpm_pack_all_arrays, lpm_pack_arrays, lpm_sort
    79 
    80     INTERFACE lpm_pack_all_arrays
    81        MODULE PROCEDURE lpm_pack_all_arrays
    82     END INTERFACE lpm_pack_all_arrays
    83 
    84     INTERFACE lpm_pack_arrays
    85        MODULE PROCEDURE lpm_pack_arrays
    86     END INTERFACE lpm_pack_arrays
    87 
    88     INTERFACE lpm_sort
    89        MODULE PROCEDURE lpm_sort
    90     END INTERFACE lpm_sort
     85    PUBLIC lpm_sort_in_subboxes, lpm_pack, lpm_sort_timeloop_done
     86
     87    INTERFACE lpm_sort_in_subboxes
     88       MODULE PROCEDURE lpm_sort_in_subboxes
     89    END INTERFACE lpm_sort_in_subboxes
     90
     91    INTERFACE lpm_pack
     92       MODULE PROCEDURE lpm_pack
     93    END INTERFACE lpm_pack
     94
     95    INTERFACE lpm_sort_timeloop_done
     96       MODULE PROCEDURE lpm_sort_timeloop_done
     97    END INTERFACE lpm_sort_timeloop_done
    9198
    9299
    93100 CONTAINS
     101
     102!------------------------------------------------------------------------------!
     103! Description:
     104! -----------
     105!> Routine for the whole processor
     106!> Sort all particles into the 8 respective subgrid boxes
     107!------------------------------------------------------------------------------!
     108    SUBROUTINE lpm_sort_in_subboxes
     109
     110       USE cpulog,                                                              &
     111           ONLY:  cpu_log, log_point_s
     112
     113       USE indices,                                                             &
     114           ONLY:  nxl, nxr, nys, nyn, nzb, nzt
     115
     116       USE kinds
     117
     118       USE control_parameters,                                                  &
     119           ONLY: dz
     120       
     121       USE grid_variables,                                                      &
     122           ONLY: dx,dy,ddx, ddy
     123         
     124       IMPLICIT NONE
     125
     126       INTEGER(iwp) ::  ip !<
     127       INTEGER(iwp) ::  jp !<
     128       INTEGER(iwp) ::  kp !<
     129
     130       CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'start' )
     131       DO  ip = nxl, nxr
     132          DO  jp = nys, nyn
     133             DO  kp = nzb+1, nzt
     134                number_of_particles = prt_count(kp,jp,ip)
     135                IF ( number_of_particles <= 0 )  CYCLE
     136                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     137               
     138                CALL lpm_pack_and_sort(ip,jp,kp)
     139               
     140                prt_count(kp,jp,ip) = number_of_particles
     141             ENDDO
     142          ENDDO
     143       ENDDO
     144       CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'stop' )
     145       RETURN
     146
     147    END SUBROUTINE lpm_sort_in_subboxes
    94148
    95149!------------------------------------------------------------------------------!
    96150! Description:
    97151! ------------
    98 !> @todo Missing subroutine description.
    99 !------------------------------------------------------------------------------!
    100     SUBROUTINE lpm_pack_all_arrays
    101 
    102        USE cpulog,                                                             &
    103            ONLY:  cpu_log, log_point_s
    104 
    105        USE indices,                                                            &
    106            ONLY:  nxl, nxr, nys, nyn, nzb, nzt
    107 
    108        USE kinds
    109 
    110        IMPLICIT NONE
    111 
    112        INTEGER(iwp) ::  i !<
    113        INTEGER(iwp) ::  j !<
    114        INTEGER(iwp) ::  k !<
    115 
    116        CALL cpu_log( log_point_s(51), 'lpm_pack_all_arrays', 'start' )
    117        DO  i = nxl, nxr
    118           DO  j = nys, nyn
    119              DO  k = nzb+1, nzt
    120                 number_of_particles = prt_count(k,j,i)
    121                 IF ( number_of_particles <= 0 )  CYCLE
    122                 particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    123                 CALL lpm_pack_and_sort(i,j,k)
    124                 prt_count(k,j,i) = number_of_particles
    125              ENDDO
    126           ENDDO
    127        ENDDO
    128        CALL cpu_log( log_point_s(51), 'lpm_pack_all_arrays', 'stop' )
    129        RETURN
    130 
    131     END SUBROUTINE lpm_pack_all_arrays
    132 
    133 !------------------------------------------------------------------------------!
    134 ! Description:
    135 ! ------------
    136 !> @todo Missing subroutine description.
    137 !------------------------------------------------------------------------------!
    138     SUBROUTINE lpm_pack_arrays
     152!> Move all particles not marked for deletion to lowest indices (packing)
     153!------------------------------------------------------------------------------!
     154    SUBROUTINE lpm_pack
    139155
    140156       USE kinds
     
    173189       number_of_particles = nn
    174190
    175     END SUBROUTINE lpm_pack_arrays
    176 
    177 !------------------------------------------------------------------------------!
    178 ! Description:
    179 ! ------------
    180 !> @todo Missing subroutine description.
    181 !------------------------------------------------------------------------------!
     191    END SUBROUTINE lpm_pack
     192
    182193    SUBROUTINE lpm_pack_and_sort (ip,jp,kp)
    183194
    184       USE control_parameters,                                                  &
    185           ONLY: dz
    186 
    187       USE kinds
    188 
    189       USE grid_variables,                                                      &
    190           ONLY: ddx, ddy
    191 
    192       IMPLICIT NONE
    193 
    194       INTEGER(iwp), INTENT(IN) :: ip
    195       INTEGER(iwp), INTENT(IN) :: jp
    196       INTEGER(iwp), INTENT(IN) :: kp
    197 
    198       INTEGER(iwp)             :: i
    199       INTEGER(iwp)             :: is
    200       INTEGER(iwp)             :: j
    201       INTEGER(iwp)             :: k
    202       INTEGER(iwp)             :: n
    203       INTEGER(iwp)             :: nn
    204       INTEGER(iwp)             :: m
    205       INTEGER(iwp)             :: sort_index
    206 
    207       INTEGER(iwp), DIMENSION(0:7) :: sort_count
    208 
    209       TYPE(particle_type), DIMENSION(number_of_particles,0:7) :: sort_particles
    210 
     195       USE control_parameters,                                                 &
     196           ONLY: dz
     197
     198       USE kinds
     199
     200       USE grid_variables,                                                     &
     201           ONLY: dx,dy,ddx, ddy
     202
     203       IMPLICIT NONE
     204
     205       INTEGER(iwp), INTENT(IN) :: ip
     206       INTEGER(iwp), INTENT(IN) :: jp
     207       INTEGER(iwp), INTENT(IN) :: kp
     208
     209       INTEGER(iwp)             :: i
     210       INTEGER(iwp)             :: is
     211       INTEGER(iwp)             :: j
     212       INTEGER(iwp)             :: k
     213       INTEGER(iwp)             :: n
     214       INTEGER(iwp)             :: nn
     215       INTEGER(iwp)             :: m
     216       INTEGER(iwp)             :: sort_index
     217
     218       INTEGER(iwp), DIMENSION(0:7) :: sort_count
     219
     220       TYPE(particle_type), DIMENSION(number_of_particles,0:7) :: sort_particles
     221     
    211222       nn = 0
    212223       sort_count = 0
    213 
     224               
    214225       DO  n = 1, number_of_particles
    215226          sort_index = 0
     
    217228          IF ( particles(n)%particle_mask )  THEN
    218229             nn = nn + 1
    219              i = particles(n)%x * ddx
    220              j = particles(n)%y * ddy
    221              k = ( particles(n)%z + 0.5_wp * dz ) / dz + offset_ocean_nzt
    222              IF ( i == ip )  sort_index = sort_index+4
    223              IF ( j == jp )  sort_index = sort_index+2
    224              IF ( k == kp )  sort_index = sort_index+1
    225              sort_count(sort_index) = sort_count(sort_index)+1
     230!
     231!--          Sorting particles with a binary scheme
     232!--          sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox
     233!--          sort_index=000_2=0_10 -> particle at the right,north,top subgridbox
     234!--          For this the center of the gridbox is calculated
     235             i = (particles(n)%x + 0.5_wp * dx) * ddx
     236             j = (particles(n)%y + 0.5_wp * dy) * ddy
     237             k = ( particles(n)%z+ 0.5_wp *dz ) / dz+1 + offset_ocean_nzt
     238                     
     239             IF ( i == ip )  sort_index = sort_index + 4
     240             IF ( j == jp )  sort_index = sort_index + 2
     241             IF ( k == kp )  sort_index = sort_index + 1
     242                     
     243             sort_count(sort_index) = sort_count(sort_index) + 1
    226244             m = sort_count(sort_index)
    227245             sort_particles(m,sort_index) = particles(n)
     
    236254          grid_particles(kp,jp,ip)%start_index(is) = nn + 1
    237255          DO n = 1,sort_count(is)
    238              nn = nn+1
     256             nn = nn + 1
    239257             particles(nn) = sort_particles(n,is)
    240           ENDDO
    241           grid_particles(kp,jp,ip)%end_index(is) = nn
    242        ENDDO
    243 
    244        number_of_particles = nn
    245        RETURN
    246 
     258             ENDDO
     259             grid_particles(kp,jp,ip)%end_index(is) = nn
     260       ENDDO
     261
     262       number_of_particles = nn   
     263               
    247264    END SUBROUTINE lpm_pack_and_sort
    248 
     265               
    249266!------------------------------------------------------------------------------!
    250267! Description:
     
    254271!> complete the LES timestep.
    255272!------------------------------------------------------------------------------!
    256     SUBROUTINE lpm_sort
     273    SUBROUTINE lpm_sort_timeloop_done
    257274
    258275       USE control_parameters,                                                 &
     
    342359       ENDDO
    343360
    344     END SUBROUTINE lpm_sort
    345 
    346 
    347  END module lpm_pack_arrays_mod
     361    END SUBROUTINE lpm_sort_timeloop_done
     362
     363
     364 END module lpm_pack_and_sort_mod
  • palm/trunk/SOURCE/lpm_read_restart_file.f90

    r2312 r2606  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 !
     22! 
     23! 
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Changed particle box locations: center of particle box now coincides
     28! with scalar grid point of same index.
     29! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
     30! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     31! lpm_sort -> lpm_sort_timeloop_done
     32!
     33! 2312 2017-07-14 20:26:51Z hoffmann
    2734! Extended particle data type.
    2835!
     
    7481    USE kinds
    7582
    76     USE lpm_pack_arrays_mod,                                                   &
    77         ONLY:  lpm_pack_all_arrays
     83    USE lpm_pack_and_sort_mod,                                                   &
     84        ONLY:  lpm_sort_in_subboxes
    7885
    7986    USE particle_attributes,                                                   &
     
    178185!-- Must be called to sort particles into blocks, which is needed for a fast
    179186!-- interpolation of the LES fields on the particle position.
    180     CALL lpm_pack_all_arrays
     187    CALL lpm_sort_in_subboxes
    181188
    182189
Note: See TracChangeset for help on using the changeset viewer.