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

Modified particle box location and further changes in particle model

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.