Changeset 1929


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

several bugfixes in particle model and serial mode

Location:
palm/trunk/SOURCE
Files:
11 edited

Legend:

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

    r1919 r1929  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix in check for use_upstream_for_tke
    2222!
    2323! Former revisions:
     
    894894    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.               &
    895895         .NOT. use_upstream_for_tke  .AND.                                       &
    896          ( scalar_advec /= 'ws-scheme'  .OR.  scalar_advec /= 'ws-scheme-mono' ) &
     896         ( scalar_advec /= 'ws-scheme'  .AND.  scalar_advec /= 'ws-scheme-mono' )&
    897897       )  THEN
    898898       use_upstream_for_tke = .TRUE.
  • palm/trunk/SOURCE/lpm.f90

    r1823 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Call wall boundary conditions only if particles are in the vertical range of
     22! topography.
    2223!
    2324! Former revisions:
     
    111112
    112113    USE indices,                                                               &
    113         ONLY: nxl, nxr, nys, nyn, nzb, nzt
     114        ONLY: nxl, nxr, nys, nyn, nzb, nzb_max, nzt, nzb_w_inner
    114115
    115116    USE kinds
     
    300301                CALL lpm_advec(i,j,k)
    301302!
    302 !--             Particle reflection from walls
    303                 IF ( topography /= 'flat' )  THEN
     303!--             Particle reflection from walls. Only applied if the particles
     304!--             are in the vertical range of the topography. (Here, some
     305!--             optimization is still possible.)
     306                IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
    304307                   CALL lpm_boundary_conds( 'walls' )
    305308                ENDIF
     
    334337       steps = steps + 1
    335338       dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
    336 
    337339!
    338340!--    Find out, if all particles on every PE have completed the LES timestep
     
    359361!--    Horizontal boundary conditions including exchange between subdmains
    360362       CALL lpm_exchange_horiz
    361 
    362363!
    363364!--    Pack particles (eliminate those marked for deletion),
  • palm/trunk/SOURCE/lpm_advec.f90

    r1889 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Put stochastic equation in an extra subroutine.
     22! Set flag for stochastic equation to communicate whether a particle is near
     23! topography. This case, memory and drift term are disabled in the Weil equation.
     24!
     25! Enable vertical logarithmic interpolation also above topography. This case,
     26! set a lower limit for the friction velocity, as it can become very small
     27! in narrow street canyons, leading to too large particle speeds.
    2228!
    2329! Former revisions:
     
    110116        ONLY:  block_offset, c_0, dt_min_part, grid_particles,                 &
    111117               iran_part, log_z_z0, number_of_particles, number_of_sublayers,  &
    112                particles, particle_groups, offset_ocean_nzt, sgs_wfu_part,     &
    113                sgs_wfv_part, sgs_wfw_part, use_sgs_for_particles,              &
    114                vertical_particle_advection, z0_av_global
     118               particles, particle_groups, offset_ocean_nzt, sgs_wf_part,      &
     119               use_sgs_for_particles, vertical_particle_advection, z0_av_global
    115120       
    116121    USE statistics,                                                            &
     
    119124    IMPLICIT NONE
    120125
    121     INTEGER(iwp) ::  agp                         !<
    122     INTEGER(iwp) ::  gp_outside_of_building(1:8) !<
    123     INTEGER(iwp) ::  i                           !<
    124     INTEGER(iwp) ::  ip                          !<
    125     INTEGER(iwp) ::  j                           !<
    126     INTEGER(iwp) ::  jp                          !<
    127     INTEGER(iwp) ::  k                           !<
    128     INTEGER(iwp) ::  kp                          !<
    129     INTEGER(iwp) ::  kw                          !<
    130     INTEGER(iwp) ::  n                           !<
    131     INTEGER(iwp) ::  nb                          !<
    132     INTEGER(iwp) ::  num_gp                      !<
    133 
    134     INTEGER(iwp), DIMENSION(0:7) ::  start_index !<
    135     INTEGER(iwp), DIMENSION(0:7) ::  end_index   !<
    136 
    137     REAL(wp) ::  aa                 !<
    138     REAL(wp) ::  bb                 !<
    139     REAL(wp) ::  cc                 !<
    140     REAL(wp) ::  d_sum              !<
    141     REAL(wp) ::  d_z_p_z0           !<
    142     REAL(wp) ::  dd                 !<
    143     REAL(wp) ::  de_dx_int_l        !<
    144     REAL(wp) ::  de_dx_int_u        !<
    145     REAL(wp) ::  de_dy_int_l        !<
    146     REAL(wp) ::  de_dy_int_u        !<
    147     REAL(wp) ::  de_dt              !<
    148     REAL(wp) ::  de_dt_min          !<
    149     REAL(wp) ::  de_dz_int_l        !<
    150     REAL(wp) ::  de_dz_int_u        !<
     126    INTEGER(iwp) ::  agp                         !< loop variable
     127    INTEGER(iwp) ::  gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
     128    INTEGER(iwp) ::  i                           !< index variable along x
     129    INTEGER(iwp) ::  ip                          !< index variable along x
     130    INTEGER(iwp) ::  ilog                        !< index variable along x
     131    INTEGER(iwp) ::  j                           !< index variable along y
     132    INTEGER(iwp) ::  jp                          !< index variable along y
     133    INTEGER(iwp) ::  jlog                        !< index variable along y
     134    INTEGER(iwp) ::  k                           !< index variable along z
     135    INTEGER(iwp) ::  kp                          !< index variable along z
     136    INTEGER(iwp) ::  kw                          !< index variable along z
     137    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
     138    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
     139    INTEGER(iwp) ::  num_gp                      !< number of adjacent grid points inside topography
     140
     141    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
     142    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
     143
     144    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
     145    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
     146    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
     147    REAL(wp) ::  d_sum              !< dummy argument for horizontal particle interpolation in case of topography
     148    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation
     149    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation
     150    REAL(wp) ::  de_dx_int_l        !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
     151    REAL(wp) ::  de_dx_int_u        !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
     152    REAL(wp) ::  de_dy_int_l        !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
     153    REAL(wp) ::  de_dy_int_u        !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
     154    REAL(wp) ::  de_dt              !< temporal derivative of TKE experienced by the particle
     155    REAL(wp) ::  de_dt_min          !< lower level for temporal TKE derivative
     156    REAL(wp) ::  de_dz_int_l        !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
     157    REAL(wp) ::  de_dz_int_u        !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
    151158    REAL(wp) ::  diameter           !< diamter of droplet
    152     REAL(wp) ::  diss_int_l         !<
    153     REAL(wp) ::  diss_int_u         !<
    154     REAL(wp) ::  dt_gap             !<
    155     REAL(wp) ::  dt_particle_m      !<
    156     REAL(wp) ::  e_int_l            !<
    157     REAL(wp) ::  e_int_u            !<
    158     REAL(wp) ::  e_mean_int         !<
     159    REAL(wp) ::  diss_int_l         !< x/y-interpolated dissipation at particle position at lower vertical level
     160    REAL(wp) ::  diss_int_u         !< x/y-interpolated dissipation at particle position at upper vertical level
     161    REAL(wp) ::  dt_gap             !< remaining time until particle time integration reaches LES time
     162    REAL(wp) ::  dt_particle_m      !< previous particle time step
     163    REAL(wp) ::  e_int_l            !< x/y-interpolated TKE at particle position at lower vertical level
     164    REAL(wp) ::  e_int_u            !< x/y-interpolated TKE at particle position at upper vertical level
     165    REAL(wp) ::  e_mean_int         !< horizontal mean TKE at particle height
    159166    REAL(wp) ::  exp_arg            !<
    160167    REAL(wp) ::  exp_term           !<
    161     REAL(wp) ::  gg                 !<
    162     REAL(wp) ::  height_p           !<
     168    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
     169    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
    163170    REAL(wp) ::  lagr_timescale     !< Lagrangian timescale
    164     REAL(wp) ::  location(1:30,1:3) !<
     171    REAL(wp) ::  location(1:30,1:3) !< wall locations
     172    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
    165173    REAL(wp) ::  random_gauss       !<
    166174    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
     
    169177    REAL(wp) ::  rg3                !< Gaussian distributed random number
    170178    REAL(wp) ::  sigma              !< velocity standard deviation
    171     REAL(wp) ::  u_int_l            !<
    172     REAL(wp) ::  u_int_u            !<
    173     REAL(wp) ::  us_int             !<
    174     REAL(wp) ::  v_int_l            !<
    175     REAL(wp) ::  v_int_u            !<
     179    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
     180    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
     181    REAL(wp) ::  us_int             !< friction velocity at particle grid box
     182    REAL(wp) ::  v_int_l            !< x/y-interpolated v-component at particle position at lower vertical level
     183    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
    176184    REAL(wp) ::  vv_int             !<
    177     REAL(wp) ::  w_int_l            !<
    178     REAL(wp) ::  w_int_u            !<
     185    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
     186    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
    179187    REAL(wp) ::  w_s                !< terminal velocity of droplets
    180     REAL(wp) ::  x                  !<
    181     REAL(wp) ::  y                  !<
    182     REAL(wp) ::  z_p                !<
     188    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation
     189    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
     190    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
    183191
    184192    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
     
    189197    REAL(wp), PARAMETER ::  d0_rog = 0.745_wp    !< separation diameter
    190198
    191     REAL(wp), DIMENSION(1:30) ::  d_gp_pl !<
    192     REAL(wp), DIMENSION(1:30) ::  de_dxi  !<
    193     REAL(wp), DIMENSION(1:30) ::  de_dyi  !<
    194     REAL(wp), DIMENSION(1:30) ::  de_dzi  !<
    195     REAL(wp), DIMENSION(1:30) ::  dissi   !<
    196     REAL(wp), DIMENSION(1:30) ::  ei      !<
    197 
     199    REAL(wp), DIMENSION(1:30) ::  d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
     200    REAL(wp), DIMENSION(1:30) ::  de_dxi  !< horizontal TKE gradient along x at adjacent wall
     201    REAL(wp), DIMENSION(1:30) ::  de_dyi  !< horizontal TKE gradient along y at adjacent wall
     202    REAL(wp), DIMENSION(1:30) ::  de_dzi  !< horizontal TKE gradient along z at adjacent wall
     203    REAL(wp), DIMENSION(1:30) ::  dissi   !< dissipation at adjacent wall
     204    REAL(wp), DIMENSION(1:30) ::  ei      !< TKE at adjacent wall
     205
     206    REAL(wp), DIMENSION(number_of_particles) ::  term_1_2     !< flag to communicate whether a particle is near topography or not
    198207    REAL(wp), DIMENSION(number_of_particles) ::  dens_ratio   !<
    199     REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !<
    200     REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !<
    201     REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !<
    202     REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !<
    203     REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !<
    204     REAL(wp), DIMENSION(number_of_particles) ::  e_int        !<
    205     REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !<
    206     REAL(wp), DIMENSION(number_of_particles) ::  log_z_z0_int !<
    207     REAL(wp), DIMENSION(number_of_particles) ::  u_int        !<
    208     REAL(wp), DIMENSION(number_of_particles) ::  v_int        !<
    209     REAL(wp), DIMENSION(number_of_particles) ::  w_int        !<
    210     REAL(wp), DIMENSION(number_of_particles) ::  xv           !<
    211     REAL(wp), DIMENSION(number_of_particles) ::  yv           !<
    212     REAL(wp), DIMENSION(number_of_particles) ::  zv           !<
    213 
    214     REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !<
     208    REAL(wp), DIMENSION(number_of_particles) ::  de_dx_int    !< horizontal TKE gradient along x at particle position
     209    REAL(wp), DIMENSION(number_of_particles) ::  de_dy_int    !< horizontal TKE gradient along y at particle position
     210    REAL(wp), DIMENSION(number_of_particles) ::  de_dz_int    !< horizontal TKE gradient along z at particle position
     211    REAL(wp), DIMENSION(number_of_particles) ::  diss_int     !< dissipation at particle position
     212    REAL(wp), DIMENSION(number_of_particles) ::  dt_particle  !< particle time step
     213    REAL(wp), DIMENSION(number_of_particles) ::  e_int        !< TKE at particle position
     214    REAL(wp), DIMENSION(number_of_particles) ::  fs_int       !< weighting factor for subgrid-scale particle speed
     215    REAL(wp), DIMENSION(number_of_particles) ::  u_int        !< u-component of particle speed
     216    REAL(wp), DIMENSION(number_of_particles) ::  v_int        !< v-component of particle speed
     217    REAL(wp), DIMENSION(number_of_particles) ::  w_int        !< w-component of particle speed
     218    REAL(wp), DIMENSION(number_of_particles) ::  xv           !< x-position
     219    REAL(wp), DIMENSION(number_of_particles) ::  yv           !< y-position
     220    REAL(wp), DIMENSION(number_of_particles) ::  zv           !< z-position
     221
     222    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
    215223
    216224    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
     
    236244       j = jp + block_offset(nb)%j_off
    237245       k = kp + block_offset(nb)%k_off
     246
    238247
    239248!
     
    249258!--       First, check if particle is located below first vertical grid level
    250259!--       (Prandtl-layer height)
    251           IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
     260          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
     261          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
     262
     263          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
    252264!
    253265!--          Resolved-scale horizontal particle velocity is zero below z0.
    254              IF ( particles(n)%z < z0_av_global )  THEN
     266             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
    255267                u_int(n) = 0.0_wp
    256268             ELSE
    257269!
    258270!--             Determine the sublayer. Further used as index.
    259                 height_p     = ( particles(n)%z - z0_av_global )            &
     271                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global )            &
    260272                                     * REAL( number_of_sublayers, KIND=wp ) &
    261273                                     * d_z_p_z0
     
    263275!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
    264276!--             interpolate linearly between precalculated logarithm.
    265                 log_z_z0_int(n) = log_z_z0(INT(height_p))                      &
     277                log_z_z0_int = log_z_z0(INT(height_p))                         &
    266278                                 + ( height_p - INT(height_p) )                &
    267279                                 * ( log_z_z0(INT(height_p)+1)                 &
    268280                                      - log_z_z0(INT(height_p))                &
    269281                                   )
     282!
     283!--             Limit friction velocity. In narrow canyons or holes the
     284!--             friction velocity can become very small, resulting in a too
     285!--             large particle speed.
     286                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ),  &
     287                                0.01_wp ) 
    270288!
    271289!--             Neutral solution is applied for all situations, e.g. also for
     
    275293!--             as sensitivity studies revealed no significant effect of
    276294!--             using the neutral solution also for un/stable situations.
    277 !--             Calculated left and bottom index on u grid.
    278                 us_int   = 0.5_wp * ( us(j,i) + us(j,i-1) ) 
    279 
    280                 u_int(n) = -usws(j,i) / ( us_int * kappa + 1E-10_wp )          &
    281                          * log_z_z0_int(n) - u_gtrans
    282 
     295                u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          &
     296                            * log_z_z0_int - u_gtrans
     297               
    283298             ENDIF
    284299!
     
    308323                           ( u_int_u - u_int_l )
    309324             ENDIF
     325
    310326          ENDIF
    311327
     
    319335       DO  n = start_index(nb), end_index(nb)
    320336
    321           IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
    322 
    323              IF ( particles(n)%z < z0_av_global )  THEN
     337          ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
     338          jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
     339          IF ( constant_flux_layer .AND. zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p )  THEN
     340
     341             IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global )  THEN
    324342!
    325343!--             Resolved-scale horizontal particle velocity is zero below z0.
    326344                v_int(n) = 0.0_wp
    327345             ELSE       
     346!
     347!--             Determine the sublayer. Further used as index. Please note,
     348!--             logarithmus can not be reused from above, as in in case of
     349!--             topography particle on u-grid can be above surface-layer height,
     350!--             whereas it can be below on v-grid.
     351                height_p     = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) &
     352                                  * REAL( number_of_sublayers, KIND=wp )    &
     353                                  * d_z_p_z0
     354!
     355!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
     356!--             interpolate linearly between precalculated logarithm.
     357                log_z_z0_int = log_z_z0(INT(height_p))                         &
     358                                 + ( height_p - INT(height_p) )                &
     359                                 * ( log_z_z0(INT(height_p)+1)                 &
     360                                      - log_z_z0(INT(height_p))                &
     361                                   )
     362!
     363!--             Limit friction velocity. In narrow canyons or holes the
     364!--             friction velocity can become very small, resulting in a too
     365!--             large particle speed.
     366                us_int   = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), &
     367                                0.01_wp )   
    328368!
    329369!--             Neutral solution is applied for all situations, e.g. also for
     
    333373!--             as sensitivity studies revealed no significant effect of
    334374!--             using the neutral solution also for un/stable situations.
    335 !--             Calculated left and bottom index on v grid.
    336                 us_int   = 0.5_wp * ( us(j,i) + us(j-1,i) )   
    337 
    338                 v_int(n) = -vsws(j,i) / ( us_int * kappa + 1E-10_wp )          &
    339                          * log_z_z0_int(n) - v_gtrans
    340              ENDIF
     375                v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp )          &
     376                         * log_z_z0_int - v_gtrans
     377
     378             ENDIF
     379
    341380          ELSE
    342381             x  = xv(n) - i * dx
     
    361400                                  ( v_int_u - v_int_l )
    362401             ENDIF
     402
    363403          ENDIF
    364404
     
    367407       i = ip + block_offset(nb)%i_off
    368408       j = jp + block_offset(nb)%j_off
    369        k = kp-1
     409       k = kp - 1
    370410!
    371411!--    Same procedure for interpolation of the w velocity-component
     
    535575                ENDIF
    536576
     577!
     578!--             Set flag for stochastic equation.
     579                term_1_2(n) = 1.0_wp
     580
    537581             ENDDO
    538582          ENDDO
     
    541585
    542586          DO  n = 1, number_of_particles
    543 
    544587             i = particles(n)%x * ddx
    545588             j = particles(n)%y * ddy
     
    575618                de_dzi(num_gp) = de_dz(k,j,i)
    576619             ENDIF
    577 
    578              IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
    579              THEN
     620             IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
    580621                num_gp = num_gp + 1
    581622                gp_outside_of_building(2) = 1
     
    603644             ENDIF
    604645
    605              IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
    606              THEN
     646             IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 )  THEN
    607647                num_gp = num_gp + 1
    608648                gp_outside_of_building(4) = 1
     
    617657             ENDIF
    618658
    619              IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
    620              THEN
     659             IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
    621660                num_gp = num_gp + 1
    622661                gp_outside_of_building(5) = 1
     
    631670             ENDIF
    632671
    633              IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 ) &
    634              THEN
     672             IF ( k > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0 )  THEN
    635673                num_gp = num_gp + 1
    636674                gp_outside_of_building(6) = 1
     
    645683             ENDIF
    646684
    647              IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
    648              THEN
     685             IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 )  THEN
    649686                num_gp = num_gp + 1
    650687                gp_outside_of_building(7) = 1
     
    659696             ENDIF
    660697
    661              IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)&
    662              THEN
     698             IF ( k+1 > nzb_s_inner(j+1,i+1)  .OR.  nzb_s_inner(j+1,i+1) == 0)  THEN
    663699                num_gp = num_gp + 1
    664700                gp_outside_of_building(8) = 1
     
    672708                de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    673709             ENDIF
    674 
    675710!
    676711!--          If all gridpoints are situated outside of a building, then the
     
    685720                dd = ( dx - x )**2 + ( dy - y )**2
    686721                gg = aa + bb + cc + dd
    687 
     722 
    688723                e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
    689724                          + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
    690725                          ) / ( 3.0_wp * gg )
    691 
     726   
    692727                IF ( k == nzt )  THEN
    693728                   e_int(n) = e_int_l
     
    699734                             ) / ( 3.0_wp * gg )
    700735                   e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
    701                                      ( e_int_u - e_int_l )
    702                 ENDIF
    703 !
     736                                       ( e_int_u - e_int_l )
     737                ENDIF
     738! 
    704739!--             Needed to avoid NaN particle velocities (this might not be
    705740!--             required any more)
     
    757792                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
    758793                                 ) / ( 3.0_wp * gg )
    759 
     794 
    760795                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    761796                      de_dz_int(n) = de_dz_int_l
     
    778813                               ( gg - dd ) * diss(k,j+1,i+1) &
    779814                             ) / ( 3.0_wp * gg )
    780 
     815 
    781816                IF ( k == nzt )  THEN
    782817                   diss_int(n) = diss_int_l
     
    790825                                           ( diss_int_u - diss_int_l )
    791826                ENDIF
    792 
     827!
     828!--             Set flag for stochastic equation.
     829                term_1_2(n) = 1.0_wp
     830 
    793831             ELSE
    794 
     832 
    795833!
    796834!--             If wall between gridpoint 1 and gridpoint 5, then
     
    810848
    811849                IF ( gp_outside_of_building(5) == 1  .AND. &
    812                    gp_outside_of_building(1) == 0 )  THEN
     850                     gp_outside_of_building(1) == 0 )  THEN
    813851                   num_gp = num_gp + 1
    814852                   location(num_gp,1) = i * dx + 0.5_wp * dx
     
    893931                   de_dxi(num_gp) = de_dx(k,j,i)
    894932                   de_dyi(num_gp) = 0.0_wp
    895                    de_dzi(num_gp) = de_dz(k,j,i)
     933                   de_dzi(num_gp) = de_dz(k,j,i) 
    896934                ENDIF
    897935
     
    10761114                ENDIF
    10771115
    1078 !
     1116! 
    10791117!--             If wall between gridpoint 6 and gridpoint 8, then
    10801118!--             Neumann boundary condition has to be applied
     
    10921130                   de_dzi(num_gp) = 0.0_wp
    10931131                ENDIF
    1094 
     1132 
    10951133!
    10961134!--             Carry out the interpolation
    10971135                IF ( num_gp == 1 )  THEN
    1098 !
     1136! 
    10991137!--                If only one of the gridpoints is situated outside of the
    11001138!--                building, it follows that the values at the particle
    11011139!--                location are the same as the gridpoint values
    1102                    e_int(n)     = ei(num_gp)
    1103                    diss_int(n)  = dissi(num_gp)
     1140                   e_int(n)     = ei(num_gp)   
     1141                   diss_int(n)  = dissi(num_gp) 
    11041142                   de_dx_int(n) = de_dxi(num_gp)
    11051143                   de_dy_int(n) = de_dyi(num_gp)
    11061144                   de_dz_int(n) = de_dzi(num_gp)
     1145!
     1146!--                Set flag for stochastic equation which disables calculation
     1147!--                of drift and memory term near topography.
     1148                   term_1_2(n) = 0.0_wp
    11071149                ELSE IF ( num_gp > 1 )  THEN
    1108 
     1150 
    11091151                   d_sum = 0.0_wp
    1110 !
     1152! 
    11111153!--                Evaluation of the distances between the gridpoints
    11121154!--                contributing to the interpolated values, and the particle
     
    11181160                      d_sum        = d_sum + d_gp_pl(agp)
    11191161                   ENDDO
    1120 
     1162 
    11211163!
    11221164!--                Finally the interpolation can be carried out
    11231165                   e_int(n)     = 0.0_wp
    11241166                   diss_int(n)  = 0.0_wp
    1125                    de_dx_int(n) = 0.0_wp
    1126                    de_dy_int(n) = 0.0_wp
    1127                    de_dz_int(n) = 0.0_wp
     1167                   de_dx_int(n) = 0.0_wp  
     1168                   de_dy_int(n) = 0.0_wp  
     1169                   de_dz_int(n) = 0.0_wp  
    11281170                   DO  agp = 1, num_gp
    11291171                      e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
     
    11381180                                         de_dzi(agp) / ( (num_gp-1) * d_sum )
    11391181                   ENDDO
    1140 
    1141                 ENDIF
    1142 
     1182 
     1183                ENDIF
     1184                e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
     1185                diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
     1186                de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
     1187                de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) )
     1188                de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) )
     1189!
     1190!--             Set flag for stochastic equation which disables calculation
     1191!--             of drift and memory term near topography.
     1192                term_1_2(n) = 0.0_wp
    11431193             ENDIF
    11441194          ENDDO
     
    12091259!
    12101260!--       Calculate the Lagrangian timescale according to Weil et al. (2004).
    1211           lagr_timescale = ( 4.0_wp * e_int(n) ) / &
    1212                            ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )
     1261          lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
     1262                           ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
    12131263
    12141264!
     
    12331283!--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
    12341284!--          from becoming unrealistically large.
    1235              particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wfu_part * e_int(n) ) *   &
     1285             particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    12361286                                  ( rg(n,1) - 1.0_wp )
    1237              particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wfv_part * e_int(n) ) *   &
     1287             particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    12381288                                  ( rg(n,2) - 1.0_wp )
    1239              particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wfw_part * e_int(n) ) *   &
     1289             particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    12401290                                  ( rg(n,3) - 1.0_wp )
    12411291
     
    12671317             ENDIF
    12681318
    1269              particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 *       &
    1270                            diss_int(n) * particles(n)%rvar1 * dt_particle(n) / &
    1271                            ( 4.0_wp * sgs_wfu_part * e_int(n) ) +              &
    1272                            ( 2.0_wp * sgs_wfu_part * de_dt *                   &
    1273                              particles(n)%rvar1 /                              &
    1274                              ( 2.0_wp * sgs_wfu_part * e_int(n) ) +            &
    1275                              de_dx_int(n)                                      &
    1276                            ) * dt_particle(n) / 2.0_wp          +              &
    1277                            SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
    1278                            ( rg(n,1) - 1.0_wp ) *                              &
    1279                            SQRT( dt_particle(n) )
    1280 
    1281              particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 *       &
    1282                            diss_int(n) * particles(n)%rvar2 * dt_particle(n) / &
    1283                            ( 4.0_wp * sgs_wfv_part * e_int(n) ) +              &
    1284                            ( 2.0_wp * sgs_wfv_part * de_dt *                   &
    1285                              particles(n)%rvar2 /                              &
    1286                              ( 2.0_wp * sgs_wfv_part * e_int(n) ) +            &
    1287                              de_dy_int(n)                                      &
    1288                            ) * dt_particle(n) / 2.0_wp          +              &
    1289                            SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
    1290                            ( rg(n,2) - 1.0_wp ) *                              &
    1291                            SQRT( dt_particle(n) )
    1292 
    1293              particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 *       &
    1294                            diss_int(n) * particles(n)%rvar3 * dt_particle(n) / &
    1295                            ( 4.0_wp * sgs_wfw_part * e_int(n) ) +              &
    1296                            ( 2.0_wp * sgs_wfw_part * de_dt *                   &
    1297                              particles(n)%rvar3 /                              &
    1298                              ( 2.0_wp * sgs_wfw_part * e_int(n) ) +            &
    1299                              de_dz_int(n)                                      &
    1300                            ) * dt_particle(n) / 2.0_wp          +              &
    1301                            SQRT( fs_int(n) * c_0 * diss_int(n) ) *             &
    1302                            ( rg(n,3) - 1.0_wp ) *                              &
    1303                            SQRT( dt_particle(n) )
     1319             CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),  &
     1320                                     de_dx_int(n), de_dt, diss_int(n),         &
     1321                                     dt_particle(n), rg(n,1), term_1_2(n) )
     1322
     1323             CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),  &
     1324                                     de_dy_int(n), de_dt, diss_int(n),         &
     1325                                     dt_particle(n), rg(n,2), term_1_2(n) )
     1326
     1327             CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),  &
     1328                                     de_dz_int(n), de_dt, diss_int(n),         &
     1329                                     dt_particle(n), rg(n,3), term_1_2(n) )
    13041330
    13051331          ENDIF
     1332
    13061333          u_int(n) = u_int(n) + particles(n)%rvar1
    13071334          v_int(n) = v_int(n) + particles(n)%rvar2
    13081335          w_int(n) = w_int(n) + particles(n)%rvar3
    1309 
    13101336!
    13111337!--       Store the SGS TKE of the current timelevel which is needed for
     
    15071533    CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
    15081534
     1535
    15091536 END SUBROUTINE lpm_advec
     1537
     1538! Description:
     1539! ------------
     1540!> Calculation of subgrid-scale particle speed using the stochastic model
     1541!> of Weil et al. (2004, JAS, 61, 2877-2887).
     1542!------------------------------------------------------------------------------!
     1543 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
     1544                                dt_n, rg_n, fac )
     1545
     1546    USE kinds
     1547
     1548    USE particle_attributes,                                                   &
     1549        ONLY:  c_0, sgs_wf_part
     1550
     1551    IMPLICIT NONE
     1552
     1553    REAL(wp) ::  a1      !< dummy argument
     1554    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position
     1555    REAL(wp) ::  dedxi_n !< horizontal derivative of TKE at particle position
     1556    REAL(wp) ::  diss_n  !< dissipation at particle position
     1557    REAL(wp) ::  dt_n    !< particle timestep
     1558    REAL(wp) ::  e_n     !< TKE at particle position
     1559    REAL(wp) ::  fac     !< flag to identify adjacent topography
     1560    REAL(wp) ::  fs_n    !< weighting factor to prevent that subgrid-scale particle speed becomes too large
     1561    REAL(wp) ::  sgs_w   !< constant (1/3)
     1562    REAL(wp) ::  rg_n    !< random number
     1563    REAL(wp) ::  term1   !< memory term
     1564    REAL(wp) ::  term2   !< drift correction term
     1565    REAL(wp) ::  term3   !< random term
     1566    REAL(wp) ::  v_sgs   !< subgrid-scale velocity component
     1567
     1568!
     1569!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
     1570!-- multiplied by a flag to switch of both terms near topography.
     1571!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
     1572!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
     1573!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
     1574!-- disabled if one of the adjacent grid points belongs to topography.
     1575!-- Moreover, in this case, the  previous subgrid-scale component is also set
     1576!-- to zero.
     1577
     1578    a1 = fs_n * c_0 * diss_n
     1579!
     1580!-- Memory term
     1581    term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp )    &
     1582                 * fac
     1583!
     1584!-- Drift correction term
     1585    term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n              &
     1586                 * fac
     1587!
     1588!-- Random term
     1589    term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
     1590!
     1591!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
     1592!-- subgrid-scale velocity component is set to zero, in order to prevent a
     1593!-- velocity build-up.
     1594
     1595!-- This case, set also previous subgrid-scale component to zero.
     1596    v_sgs = v_sgs * fac + term1 + term2 + term3
     1597
     1598 END SUBROUTINE weil_stochastic_eq
  • palm/trunk/SOURCE/lpm_boundary_conds.f90

    r1823 r1929  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Rewritten wall reflection
    2222!
    2323! Former revisions:
     
    102102    CHARACTER (LEN=*) ::  range     !<
    103103   
    104     INTEGER(iwp) ::  i              !<
    105     INTEGER(iwp) ::  inc            !<
    106     INTEGER(iwp) ::  ir             !<
    107     INTEGER(iwp) ::  i1             !<
    108     INTEGER(iwp) ::  i2             !<
    109     INTEGER(iwp) ::  i3             !<
    110     INTEGER(iwp) ::  i5             !<
    111     INTEGER(iwp) ::  j              !<
    112     INTEGER(iwp) ::  jr             !<
    113     INTEGER(iwp) ::  j1             !<
    114     INTEGER(iwp) ::  j2             !<
    115     INTEGER(iwp) ::  j3             !<
    116     INTEGER(iwp) ::  j5             !<
    117     INTEGER(iwp) ::  k              !<
    118     INTEGER(iwp) ::  k1             !<
    119     INTEGER(iwp) ::  k2             !<
    120     INTEGER(iwp) ::  k3             !<
    121     INTEGER(iwp) ::  k5             !<
    122     INTEGER(iwp) ::  n              !<
    123     INTEGER(iwp) ::  t_index        !<
    124     INTEGER(iwp) ::  t_index_number !<
     104    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
     105    INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
     106    INTEGER(iwp) ::  i1             !< grid index (x) of old particle position
     107    INTEGER(iwp) ::  i2             !< grid index (x) of current particle position
     108    INTEGER(iwp) ::  i3             !< grid index (x) of intermediate particle position
     109    INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
     110    INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
     111    INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
     112    INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
     113    INTEGER(iwp) ::  n              !< particle number
     114    INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
     115    INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
     116    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
     117    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
     118
     119    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
     120    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
    125121   
    126     LOGICAL  ::  reflect_x   !<
    127     LOGICAL  ::  reflect_y   !<
    128     LOGICAL  ::  reflect_z   !<
    129 
    130     REAL(wp) ::  dt_particle !<
    131     REAL(wp) ::  pos_x       !<
    132     REAL(wp) ::  pos_x_old   !<
    133     REAL(wp) ::  pos_y       !<
    134     REAL(wp) ::  pos_y_old   !<
    135     REAL(wp) ::  pos_z       !<
    136     REAL(wp) ::  pos_z_old   !<
    137     REAL(wp) ::  prt_x       !<
    138     REAL(wp) ::  prt_y       !<
    139     REAL(wp) ::  prt_z       !<
    140     REAL(wp) ::  t(1:200)    !<
    141     REAL(wp) ::  tmp_t       !<
    142     REAL(wp) ::  xline       !<
    143     REAL(wp) ::  yline       !<
    144     REAL(wp) ::  zline       !<
     122    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
     123    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
     124    LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
     125    LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
     126    LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
     127    LOGICAL  ::  reflect_z       !< flag to check if particle is already reflected along z
     128    LOGICAL  ::  tmp_reach_x     !< dummy for sorting algorithmus
     129    LOGICAL  ::  tmp_reach_y     !< dummy for sorting algorithmus
     130    LOGICAL  ::  tmp_reach_z     !< dummy for sorting algorithmus
     131    LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
     132    LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
     133
     134    LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
     135    LOGICAL, DIMENSION(0:10) ::  reach_y  !< flag to check if particle is at a xz-wall
     136    LOGICAL, DIMENSION(0:10) ::  reach_z  !< flag to check if particle is at a xy-wall
     137
     138    REAL(wp) ::  dt_particle    !< particle timestep
     139    REAL(wp) ::  dum            !< dummy argument
     140    REAL(wp) ::  eps = 1E-10_wp !< security number to check if particle has reached a wall
     141    REAL(wp) ::  pos_x          !< intermediate particle position (x)
     142    REAL(wp) ::  pos_x_old      !< particle position (x) at previous particle timestep
     143    REAL(wp) ::  pos_y          !< intermediate particle position (y)
     144    REAL(wp) ::  pos_y_old      !< particle position (y) at previous particle timestep
     145    REAL(wp) ::  pos_z          !< intermediate particle position (z)
     146    REAL(wp) ::  pos_z_old      !< particle position (z) at previous particle timestep
     147    REAL(wp) ::  prt_x          !< current particle position (x)
     148    REAL(wp) ::  prt_y          !< current particle position (y)
     149    REAL(wp) ::  prt_z          !< current particle position (z)
     150    REAL(wp) ::  t_old          !< previous reflection time
     151    REAL(wp) ::  tmp_t          !< dummy for sorting algorithmus
     152    REAL(wp) ::  xwall          !< location of wall in x
     153    REAL(wp) ::  ywall          !< location of wall in y
     154    REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
     155    REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
     156    REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
     157    REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
     158
     159    REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
     160
    145161
    146162    IF ( range == 'bottom/top' )  THEN
     
    211227    ELSEIF ( range == 'walls' )  THEN
    212228
     229
    213230       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
    214231
    215        reflect_x = .FALSE.
    216        reflect_y = .FALSE.
    217        reflect_z = .FALSE.
    218 
    219232       DO  n = 1, number_of_particles
    220 
     233!
     234!--       Recalculate particle timestep
    221235          dt_particle = particles(n)%age - particles(n)%age_m
    222 
     236!
     237!--       Obtain x/y indices for current particle position
    223238          i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
    224239          j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
    225           k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
    226 
     240!
     241!--       Save current particle positions
    227242          prt_x = particles(n)%x
    228243          prt_y = particles(n)%y
    229244          prt_z = particles(n)%z
    230 
    231 !
    232 !--       If the particle position is below the surface, it has to be reflected
    233           IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
    234 
    235              pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
    236              pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
    237              pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
    238              i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
    239              j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
    240              k1 = pos_z_old / dz + offset_ocean_nzt_m1
    241 
    242 !
    243 !--          Case 1
    244              IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old )&
    245              THEN
    246                 t_index = 1
    247 
    248                 DO  i = i1, i2
    249                    xline      = i * dx + 0.5_wp * dx
    250                    t(t_index) = ( xline - pos_x_old ) / &
    251                                 ( particles(n)%x - pos_x_old )
    252                    t_index    = t_index + 1
     245!
     246!--       Recalculate old particle positions
     247          pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
     248          pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
     249          pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
     250!
     251!--       Obtain x/y indices for old particle positions
     252          i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
     253          j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
     254!
     255!--       Determine horizontal as well as vertical walls at which particle can
     256!--       be potentially reflected.
     257!--       Start with walls aligned in yz layer.
     258!--       Wall to the right
     259          IF ( prt_x > pos_x_old )  THEN
     260             xwall = ( i1 + 0.5_wp ) * dx
     261!
     262!--       Wall to the left
     263          ELSE
     264             xwall = ( i1 - 0.5_wp ) * dx
     265          ENDIF
     266!
     267!--       Walls aligned in xz layer
     268!--       Wall to the north
     269          IF ( prt_y > pos_y_old )  THEN
     270             ywall = ( j1 + 0.5_wp ) * dy
     271!--       Wall to the south
     272          ELSE
     273             ywall = ( j1 - 0.5_wp ) * dy
     274          ENDIF
     275!
     276!--       Walls aligned in xy layer at which particle can be possiblly reflected
     277          zwall1 = zw(nzb_s_inner(j2,i2))
     278          zwall2 = zw(nzb_s_inner(j1,i1))
     279          zwall3 = zw(nzb_s_inner(j1,i2))
     280          zwall4 = zw(nzb_s_inner(j2,i1))
     281!
     282!--       Initialize flags to check if particle reflection is necessary
     283          downwards    = .FALSE.
     284          cross_wall_x = .FALSE.
     285          cross_wall_y = .FALSE.
     286!
     287!--       Initialize flags to check if a wall is reached
     288          reach_x      = .FALSE.
     289          reach_y      = .FALSE.
     290          reach_z      = .FALSE.
     291!
     292!--       Initialize flags to check if a particle was already reflected
     293          reflect_x = .FALSE.
     294          reflect_y = .FALSE.
     295          reflect_z = .FALSE.
     296!
     297!--       Initialize flags to check if a vertical wall is already crossed.
     298!--       ( Required to obtain correct indices. )
     299          x_wall_reached = .FALSE.
     300          y_wall_reached = .FALSE.
     301!
     302!--       Initialize time array
     303          t     = 0.0_wp
     304!
     305!--       Check if particle can reach any wall. This case, calculate the
     306!--       fractional time needed to reach this wall. Store this fractional
     307!--       timestep in array t. Moreover, store indices for these grid
     308!--       boxes where the respective wall belongs to. 
     309!--       Start with x-direction.
     310          t_index    = 1
     311          t(t_index) = ( xwall - pos_x_old )                                   &
     312                     / MERGE( MAX( prt_x - pos_x_old,  1E-30_wp ),             &
     313                              MIN( prt_x - pos_x_old, -1E-30_wp ),             &
     314                              prt_x > pos_x_old )
     315          x_ind(t_index)   = i2
     316          y_ind(t_index)   = j1
     317          reach_x(t_index) = .TRUE.
     318          reach_y(t_index) = .FALSE.
     319          reach_z(t_index) = .FALSE.
     320!
     321!--       Store these values only if particle really reaches any wall. t must
     322!--       be in a interval between [0:1].
     323          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     324             t_index      = t_index + 1
     325             cross_wall_x = .TRUE.
     326          ENDIF
     327!
     328!--       y-direction
     329          t(t_index) = ( ywall - pos_y_old )                                   &
     330                     / MERGE( MAX( prt_y - pos_y_old,  1E-30_wp ),             &
     331                              MIN( prt_y - pos_y_old, -1E-30_wp ),             &
     332                              prt_y > pos_y_old )
     333          x_ind(t_index)   = i1
     334          y_ind(t_index)   = j2
     335          reach_x(t_index) = .FALSE.
     336          reach_y(t_index) = .TRUE.
     337          reach_z(t_index) = .FALSE.
     338          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     339             t_index      = t_index + 1
     340             cross_wall_y = .TRUE.
     341          ENDIF
     342!
     343!--       z-direction
     344!--       At first, check if particle moves downwards. Only in this case a
     345!--       particle can be reflected vertically.
     346          IF ( prt_z < pos_z_old )  THEN
     347
     348             downwards = .TRUE.
     349             dum       =  1.0_wp / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ), &
     350                                          MIN( prt_z - pos_z_old, -1E-30_wp ), &
     351                                          prt_z > pos_z_old )
     352
     353             t(t_index)       = ( zwall1 - pos_z_old ) * dum
     354             x_ind(t_index)   = i2
     355             y_ind(t_index)   = j2
     356             reach_x(t_index) = .FALSE.
     357             reach_y(t_index) = .FALSE.
     358             reach_z(t_index) = .TRUE.
     359             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     360                t_index = t_index + 1
     361
     362             reach_x(t_index) = .FALSE.
     363             reach_y(t_index) = .FALSE.
     364             reach_z(t_index) = .TRUE.
     365             t(t_index)       = ( zwall2 - pos_z_old ) * dum
     366             x_ind(t_index)   = i1
     367             y_ind(t_index)   = j1
     368             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     369                t_index = t_index + 1
     370
     371             reach_x(t_index) = .FALSE.
     372             reach_y(t_index) = .FALSE.
     373             reach_z(t_index) = .TRUE.
     374             t(t_index)       = ( zwall3 - pos_z_old ) * dum
     375             x_ind(t_index)   = i2
     376             y_ind(t_index)   = j1
     377             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     378                t_index = t_index + 1
     379
     380             reach_x(t_index) = .FALSE.
     381             reach_y(t_index) = .FALSE.
     382             reach_z(t_index) = .TRUE.
     383             t(t_index)       = ( zwall4 - pos_z_old ) * dum
     384             x_ind(t_index)   = i1
     385             y_ind(t_index)   = j2
     386             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     387                t_index = t_index + 1
     388
     389          END IF
     390          t_index_number = t_index - 1
     391!
     392!--       Carry out reflection only if particle reaches any wall
     393          IF ( cross_wall_x .OR. cross_wall_y .OR. downwards )  THEN
     394!
     395!--          Sort fractional timesteps in ascending order. Also sort the
     396!--          corresponding indices and flag according to the time interval a 
     397!--          particle reaches the respective wall.
     398             inc = 1
     399             jr  = 1
     400             DO WHILE ( inc <= t_index_number )
     401                inc = 3 * inc + 1
     402             ENDDO
     403
     404             DO WHILE ( inc > 1 )
     405                inc = inc / 3
     406                DO  ir = inc+1, t_index_number
     407                   tmp_t       = t(ir)
     408                   tmp_x       = x_ind(ir)
     409                   tmp_y       = y_ind(ir)
     410                   tmp_reach_x = reach_x(ir)
     411                   tmp_reach_y = reach_y(ir)
     412                   tmp_reach_z = reach_z(ir)
     413                   jr    = ir
     414                   DO WHILE ( t(jr-inc) > tmp_t )
     415                      t(jr)       = t(jr-inc)
     416                      x_ind(jr)   = x_ind(jr-inc)
     417                      y_ind(jr)   = y_ind(jr-inc)
     418                      reach_x(jr) = reach_x(jr-inc)
     419                      reach_y(jr) = reach_y(jr-inc)
     420                      reach_z(jr) = reach_z(jr-inc)
     421                      jr    = jr - inc
     422                      IF ( jr <= inc )  EXIT
     423                   ENDDO
     424                   t(jr)       = tmp_t
     425                   x_ind(jr)   = tmp_x
     426                   y_ind(jr)   = tmp_y
     427                   reach_x(jr) = tmp_reach_x
     428                   reach_y(jr) = tmp_reach_y
     429                   reach_z(jr) = tmp_reach_z
    253430                ENDDO
    254 
    255                 DO  j = j1, j2
    256                    yline      = j * dy + 0.5_wp * dy
    257                    t(t_index) = ( yline - pos_y_old ) / &
    258                                 ( particles(n)%y - pos_y_old )
    259                    t_index    = t_index + 1
    260                 ENDDO
    261 
    262                 IF ( particles(n)%z < pos_z_old )  THEN
    263                    DO  k = k1, k2 , -1
    264                       zline      = k * dz
    265                       t(t_index) = ( pos_z_old - zline ) / &
    266                                    ( pos_z_old - particles(n)%z )
    267                       t_index    = t_index + 1
    268                    ENDDO
     431             ENDDO
     432!
     433!--          Initialize temporary particle positions
     434             pos_x = pos_x_old
     435             pos_y = pos_y_old
     436             pos_z = pos_z_old
     437!
     438!--          Loop over all times a particle possibly moves into a new grid box
     439             t_old = 0.0_wp
     440             DO t_index = 1, t_index_number
     441!           
     442!--             Calculate intermediate particle position according to the
     443!--             timesteps a particle reaches any wall.
     444                pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle           &
     445                                                       * particles(n)%speed_x
     446                pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle           &
     447                                                       * particles(n)%speed_y
     448                pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle           &
     449                                                       * particles(n)%speed_z
     450!
     451!--             Obtain x/y grid indices for intermediate particle position from
     452!--             sorted index array
     453                i3 = x_ind(t_index)
     454                j3 = y_ind(t_index)
     455!
     456!--             Check which wall is already reached
     457                IF ( .NOT. x_wall_reached )  x_wall_reached = reach_x(t_index)
     458                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index)
     459!
     460!--             Check if a particle needs to be reflected at any yz-wall. If
     461!--             necessary, carry out reflection. Please note, a security
     462!--             constant is required, as the particle position do not
     463!--             necessarily exactly match the wall location due to rounding
     464!--             errors.   
     465                IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
     466                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
     467                     reach_x(t_index)                .AND.                     &
     468                     .NOT. reflect_x )  THEN
     469!
     470!--                Reflection in x-direction.
     471!--                Ensure correct reflection by MIN/MAX functions, depending on
     472!--                direction of particle transport.
     473!--                Due to rounding errors pos_x do not exactly matches the wall
     474!--                location, leading to erroneous reflection.             
     475                   pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
     476                                  MAX( 2.0_wp * xwall - pos_x, xwall ),        &
     477                                  particles(n)%x > xwall )
     478!
     479!--                Change sign of particle speed                     
     480                   particles(n)%speed_x = - particles(n)%speed_x
     481!
     482!--                Change also sign of subgrid-scale particle speed
     483                   particles(n)%rvar1 = - particles(n)%rvar1
     484!
     485!--                Set flag that reflection along x is already done
     486                   reflect_x          = .TRUE.
     487!
     488!--                As particle do not crosses any further yz-wall during
     489!--                this timestep, set further x-indices to the current one.
     490                   x_ind(t_index:t_index_number) = i1
     491!
     492!--             If particle already reached the wall but was not reflected,
     493!--             set further x-indices to the new one.
     494                ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
     495                    x_ind(t_index:t_index_number) = i2
    269496                ENDIF
    270 
    271                 t_index_number = t_index - 1
    272 
    273 !
    274 !--             Sorting t
    275                 inc = 1
    276                 jr  = 1
    277                 DO WHILE ( inc <= t_index_number )
    278                    inc = 3 * inc + 1
    279                 ENDDO
    280 
    281                 DO WHILE ( inc > 1 )
    282                    inc = inc / 3
    283                    DO  ir = inc+1, t_index_number
    284                       tmp_t = t(ir)
    285                       jr    = ir
    286                       DO WHILE ( t(jr-inc) > tmp_t )
    287                          t(jr) = t(jr-inc)
    288                          jr    = jr - inc
    289                          IF ( jr <= inc )  EXIT
    290                       ENDDO
    291                       t(jr) = tmp_t
    292                    ENDDO
    293                 ENDDO
    294 
    295          case1: DO  t_index = 1, t_index_number
    296 
    297                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    298                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    299                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    300 
    301                    i3 = ( pos_x + 0.5_wp * dx ) * ddx   
    302                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    303                    k3 = pos_z / dz + offset_ocean_nzt_m1
    304 
    305                    i5 = pos_x * ddx
    306                    j5 = pos_y * ddy
    307                    k5 = pos_z / dz + offset_ocean_nzt_m1
    308 
    309                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    310                         nzb_s_inner(j5,i3) /= 0 )  THEN
    311 
    312                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    313                            k3 == nzb_s_inner(j5,i3) )  THEN
    314                          reflect_z = .TRUE.
    315                          EXIT case1
    316                       ENDIF
     497!
     498!--             Check if a particle needs to be reflected at any xz-wall. If
     499!--             necessary, carry out reflection.
     500                IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
     501                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
     502                     reach_y(t_index)                .AND.                     &
     503                     .NOT. reflect_y ) THEN
     504
     505                   pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
     506                                  MAX( 2.0_wp * ywall - pos_y, ywall ),        &
     507                                  particles(n)%y > ywall )
     508
     509                   particles(n)%speed_y = - particles(n)%speed_y     
     510                   particles(n)%rvar2   = - particles(n)%rvar2       
     511     
     512                   reflect_y            = .TRUE.
     513                   y_ind(t_index:t_index_number) = j1
     514
     515                ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
     516                   y_ind(t_index:t_index_number) = j2
     517                ENDIF
     518!
     519!--             Check if a particle needs to be reflected at any xy-wall. If
     520!--             necessary, carry out reflection.
     521                IF ( downwards .AND. reach_z(t_index) .AND.                    &
     522                     .NOT. reflect_z )  THEN
     523                   IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN
     524 
     525                      pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z,    &
     526                                   zw(nzb_s_inner(j3,i3)) )
     527
     528                      particles(n)%speed_z = - particles(n)%speed_z
     529                      particles(n)%rvar3   = - particles(n)%rvar3
     530
     531                      reflect_z            = .TRUE.
    317532
    318533                   ENDIF
    319534
    320                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    321                         nzb_s_inner(j3,i5) /= 0 )  THEN
    322 
    323                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    324                            k3 == nzb_s_inner(j3,i5) )  THEN
    325                          reflect_z = .TRUE.
    326                          EXIT case1
    327                       ENDIF
    328 
    329                    ENDIF
    330 
    331                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    332                         nzb_s_inner(j3,i3) /= 0 )  THEN
    333 
    334                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    335                            k3 == nzb_s_inner(j3,i3) )  THEN
    336                          reflect_z = .TRUE.
    337                          EXIT case1
    338                       ENDIF
    339 
    340                       IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    341                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    342                          reflect_y = .TRUE.
    343                          EXIT case1
    344                       ENDIF
    345 
    346                       IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    347                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    348                          reflect_x = .TRUE.
    349                          EXIT case1
    350                       ENDIF
    351 
    352                    ENDIF
    353 
    354                 ENDDO case1
    355 
    356 !
    357 !--          Case 2
    358              ELSEIF ( particles(n)%x > pos_x_old  .AND. &
    359                       particles(n)%y < pos_y_old )  THEN
    360 
    361                 t_index = 1
    362 
    363                 DO  i = i1, i2
    364                    xline      = i * dx + 0.5_wp * dx
    365                    t(t_index) = ( xline - pos_x_old ) / &
    366                                 ( particles(n)%x - pos_x_old )
    367                    t_index    = t_index + 1
    368                 ENDDO
    369 
    370                 DO  j = j1, j2, -1
    371                    yline      = j * dy - 0.5_wp * dy
    372                    t(t_index) = ( pos_y_old - yline ) / &
    373                                 ( pos_y_old - particles(n)%y )
    374                    t_index    = t_index + 1
    375                 ENDDO
    376 
    377                 IF ( particles(n)%z < pos_z_old )  THEN
    378                    DO  k = k1, k2 , -1
    379                       zline      = k * dz
    380                       t(t_index) = ( pos_z_old - zline ) / &
    381                                    ( pos_z_old - particles(n)%z )
    382                       t_index    = t_index + 1
    383                    ENDDO
    384535                ENDIF
    385                 t_index_number = t_index-1
    386 
    387 !
    388 !--             Sorting t
    389                 inc = 1
    390                 jr  = 1
    391                 DO WHILE ( inc <= t_index_number )
    392                    inc = 3 * inc + 1
    393                 ENDDO
    394 
    395                 DO WHILE ( inc > 1 )
    396                    inc = inc / 3
    397                    DO  ir = inc+1, t_index_number
    398                       tmp_t = t(ir)
    399                       jr    = ir
    400                       DO WHILE ( t(jr-inc) > tmp_t )
    401                          t(jr) = t(jr-inc)
    402                          jr    = jr - inc
    403                          IF ( jr <= inc )  EXIT
    404                       ENDDO
    405                       t(jr) = tmp_t
    406                    ENDDO
    407                 ENDDO
    408 
    409          case2: DO  t_index = 1, t_index_number
    410 
    411                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    412                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    413                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    414 
    415                    i3 = ( pos_x + 0.5_wp * dx ) * ddx
    416                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    417                    k3 = pos_z / dz + offset_ocean_nzt_m1
    418 
    419                    i5 = pos_x * ddx
    420                    j5 = pos_y * ddy
    421                    k5 = pos_z / dz + offset_ocean_nzt_m1
    422 
    423                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    424                         nzb_s_inner(j3,i5) /= 0 )  THEN
    425 
    426                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    427                            k3 == nzb_s_inner(j3,i5) )  THEN
    428                          reflect_z = .TRUE.
    429                          EXIT case2
    430                       ENDIF
    431 
    432                    ENDIF
    433 
    434                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    435                         nzb_s_inner(j3,i3) /= 0 )  THEN
    436 
    437                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    438                            k3 == nzb_s_inner(j3,i3) )  THEN
    439                          reflect_z = .TRUE.
    440                          EXIT case2
    441                       ENDIF
    442 
    443                       IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    444                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    445                          reflect_x = .TRUE.
    446                          EXIT case2
    447                       ENDIF
    448 
    449                    ENDIF
    450 
    451                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    452                         nzb_s_inner(j5,i3) /= 0 )  THEN
    453 
    454                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    455                            k3 == nzb_s_inner(j5,i3) )  THEN
    456                          reflect_z = .TRUE.
    457                          EXIT case2
    458                       ENDIF
    459 
    460                       IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    461                            pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    462                          reflect_y = .TRUE.
    463                          EXIT case2
    464                       ENDIF
    465 
    466                    ENDIF
    467 
    468                 ENDDO case2
    469 
    470 !
    471 !--          Case 3
    472              ELSEIF ( particles(n)%x < pos_x_old  .AND. &
    473                       particles(n)%y > pos_y_old )  THEN
    474 
    475                 t_index = 1
    476 
    477                 DO  i = i1, i2, -1
    478                    xline      = i * dx - 0.5_wp * dx
    479                    t(t_index) = ( pos_x_old - xline ) / &
    480                                 ( pos_x_old - particles(n)%x )
    481                    t_index    = t_index + 1
    482                 ENDDO
    483 
    484                 DO  j = j1, j2
    485                    yline      = j * dy + 0.5_wp * dy
    486                    t(t_index) = ( yline - pos_y_old ) / &
    487                                 ( particles(n)%y - pos_y_old )
    488                    t_index    = t_index + 1
    489                 ENDDO
    490 
    491                 IF ( particles(n)%z < pos_z_old )  THEN
    492                    DO  k = k1, k2 , -1
    493                       zline      = k * dz
    494                       t(t_index) = ( pos_z_old - zline ) / &
    495                                    ( pos_z_old - particles(n)%z )
    496                       t_index    = t_index + 1
    497                    ENDDO
    498                 ENDIF
    499                 t_index_number = t_index - 1
    500 
    501 !
    502 !--             Sorting t
    503                 inc = 1
    504                 jr  = 1
    505 
    506                 DO WHILE ( inc <= t_index_number )
    507                    inc = 3 * inc + 1
    508                 ENDDO
    509 
    510                 DO WHILE ( inc > 1 )
    511                    inc = inc / 3
    512                    DO  ir = inc+1, t_index_number
    513                       tmp_t = t(ir)
    514                       jr    = ir
    515                       DO WHILE ( t(jr-inc) > tmp_t )
    516                          t(jr) = t(jr-inc)
    517                          jr    = jr - inc
    518                          IF ( jr <= inc )  EXIT
    519                       ENDDO
    520                       t(jr) = tmp_t
    521                    ENDDO
    522                 ENDDO
    523 
    524          case3: DO  t_index = 1, t_index_number
    525 
    526                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    527                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    528                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    529 
    530                    i3 = ( pos_x + 0.5_wp * dx ) * ddx
    531                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    532                    k3 = pos_z / dz + offset_ocean_nzt_m1
    533 
    534                    i5 = pos_x * ddx
    535                    j5 = pos_y * ddy
    536                    k5 = pos_z / dz + offset_ocean_nzt_m1
    537 
    538                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    539                         nzb_s_inner(j5,i3) /= 0 )  THEN
    540 
    541                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    542                            k3 == nzb_s_inner(j5,i3) )  THEN
    543                          reflect_z = .TRUE.
    544                          EXIT case3
    545                       ENDIF
    546 
    547                    ENDIF
    548 
    549                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    550                         nzb_s_inner(j3,i3) /= 0 )  THEN
    551 
    552                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    553                            k3 == nzb_s_inner(j3,i3) )  THEN
    554                          reflect_z = .TRUE.
    555                          EXIT case3
    556                       ENDIF
    557 
    558                       IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    559                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    560                          reflect_y = .TRUE.
    561                          EXIT case3
    562                       ENDIF
    563 
    564                    ENDIF
    565 
    566                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    567                         nzb_s_inner(j3,i5) /= 0 )  THEN
    568 
    569                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    570                            k3 == nzb_s_inner(j3,i5) )  THEN
    571                          reflect_z = .TRUE.
    572                          EXIT case3
    573                       ENDIF
    574 
    575                       IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    576                            pos_z < nzb_s_inner(j3,i5) * dz )  THEN
    577                          reflect_x = .TRUE.
    578                          EXIT case3
    579                       ENDIF
    580 
    581                    ENDIF
    582 
    583                 ENDDO case3
    584 
    585 !
    586 !--          Case 4
    587              ELSEIF ( particles(n)%x < pos_x_old  .AND. &
    588                       particles(n)%y < pos_y_old )  THEN
    589 
    590                 t_index = 1
    591 
    592                 DO  i = i1, i2, -1
    593                    xline      = i * dx - 0.5_wp * dx
    594                    t(t_index) = ( pos_x_old - xline ) / &
    595                                 ( pos_x_old - particles(n)%x )
    596                    t_index    = t_index + 1
    597                 ENDDO
    598 
    599                 DO  j = j1, j2, -1
    600                    yline      = j * dy - 0.5_wp * dy
    601                    t(t_index) = ( pos_y_old - yline ) / &
    602                                 ( pos_y_old - particles(n)%y )
    603                    t_index    = t_index + 1
    604                 ENDDO
    605 
    606                 IF ( particles(n)%z < pos_z_old )  THEN
    607                    DO  k = k1, k2 , -1
    608                       zline      = k * dz
    609                       t(t_index) = ( pos_z_old - zline ) / &
    610                                    ( pos_z_old-particles(n)%z )
    611                       t_index    = t_index + 1
    612                    ENDDO
    613                 ENDIF
    614                 t_index_number = t_index-1
    615 
    616 !
    617 !--             Sorting t
    618                 inc = 1
    619                 jr  = 1
    620 
    621                 DO WHILE ( inc <= t_index_number )
    622                    inc = 3 * inc + 1
    623                 ENDDO
    624 
    625                 DO WHILE ( inc > 1 )
    626                    inc = inc / 3
    627                    DO  ir = inc+1, t_index_number
    628                       tmp_t = t(ir)
    629                       jr = ir
    630                       DO WHILE ( t(jr-inc) > tmp_t )
    631                          t(jr) = t(jr-inc)
    632                          jr    = jr - inc
    633                          IF ( jr <= inc )  EXIT
    634                       ENDDO
    635                       t(jr) = tmp_t
    636                    ENDDO
    637                 ENDDO
    638 
    639          case4: DO  t_index = 1, t_index_number
    640 
    641                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    642                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    643                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    644 
    645                    i3 = ( pos_x + 0.5_wp * dx ) * ddx   
    646                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    647                    k3 = pos_z / dz + offset_ocean_nzt_m1
    648 
    649                    i5 = pos_x * ddx
    650                    j5 = pos_y * ddy
    651                    k5 = pos_z / dz + offset_ocean_nzt_m1
    652 
    653                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    654                         nzb_s_inner(j3,i3) /= 0 )  THEN
    655 
    656                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    657                            k3 == nzb_s_inner(j3,i3) )  THEN
    658                          reflect_z = .TRUE.
    659                          EXIT case4
    660                       ENDIF
    661 
    662                    ENDIF
    663 
    664                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    665                         nzb_s_inner(j3,i5) /= 0 )  THEN
    666 
    667                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    668                            k3 == nzb_s_inner(j3,i5) )  THEN
    669                          reflect_z = .TRUE.
    670                          EXIT case4
    671                       ENDIF
    672 
    673                       IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    674                            nzb_s_inner(j3,i5) /=0  .AND.          &
    675                            pos_z < nzb_s_inner(j3,i5) * dz )  THEN
    676                          reflect_x = .TRUE.
    677                          EXIT case4
    678                       ENDIF
    679 
    680                    ENDIF
    681 
    682                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    683                         nzb_s_inner(j5,i3) /= 0 )  THEN
    684 
    685                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    686                            k5 == nzb_s_inner(j5,i3) )  THEN
    687                          reflect_z = .TRUE.
    688                          EXIT case4
    689                       ENDIF
    690 
    691                       IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    692                            nzb_s_inner(j5,i3) /= 0  .AND.         &
    693                            pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    694                          reflect_y = .TRUE.
    695                          EXIT case4
    696                       ENDIF
    697 
    698                    ENDIF
    699  
    700                 ENDDO case4
     536!
     537!--             Swap time
     538                t_old = t(t_index)
     539
     540             ENDDO
     541!
     542!--          If a particle was reflected, calculate final position from last
     543!--          intermediate position.
     544             IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
     545
     546                particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
     547                                                         * particles(n)%speed_x
     548                particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle      &
     549                                                         * particles(n)%speed_y
     550                particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle      &
     551                                                         * particles(n)%speed_z
    701552
    702553             ENDIF
    703554
    704           ENDIF      ! Check, if particle position is below the surface
    705 
    706 !
    707 !--       Do the respective reflection, in case that one of the above conditions
    708 !--       is found to be true
    709           IF ( reflect_z )  THEN
    710 
    711              particles(n)%z       = 2.0_wp * pos_z - prt_z
    712              particles(n)%speed_z = - particles(n)%speed_z
    713 
    714              IF ( use_sgs_for_particles )  THEN
    715                 particles(n)%rvar3 = - particles(n)%rvar3
    716              ENDIF
    717              reflect_z = .FALSE.
    718 
    719           ELSEIF ( reflect_y )  THEN
    720 
    721              particles(n)%y       = 2.0_wp * pos_y - prt_y
    722              particles(n)%speed_y = - particles(n)%speed_y
    723 
    724              IF ( use_sgs_for_particles )  THEN
    725                 particles(n)%rvar2 = - particles(n)%rvar2
    726              ENDIF
    727              reflect_y = .FALSE.
    728 
    729           ELSEIF ( reflect_x )  THEN
    730 
    731              particles(n)%x       = 2.0_wp * pos_x - prt_x
    732              particles(n)%speed_x = - particles(n)%speed_x
    733 
    734              IF ( use_sgs_for_particles )  THEN
    735                 particles(n)%rvar1 = - particles(n)%rvar1
    736              ENDIF
    737              reflect_x = .FALSE.
    738 
    739           ENDIF       
    740    
     555          ENDIF
     556
    741557       ENDDO
    742558
  • 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
  • palm/trunk/SOURCE/lpm_init.f90

    r1891 r1929  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix in determining initial particle height and grid index in case of
     22! seed_follows_topography.
     23! Bugfix concerning random positions, ensure that particles do not move more
     24! than one grid length.
     25! Bugfix logarithmic interpolation.
     26! Initial setting of sgs_wf_part.
    2227!
    2328! Former revisions:
     
    3136!
    3237! 1873 2016-04-18 14:50:06Z maronga
    33 ! Module renamed (removed _mod)
     38! Module renamed (removed _mod
     39
    3440!
    3541! 1871 2016-04-15 11:46:09Z hoffmann
     
    162168                prt_count, psb, psl, psn, psr, pss, pst,                       &
    163169                radius, random_start_position, read_particles_from_restartfile,&
    164                 seed_follows_topography, sort_count,                           &
     170                seed_follows_topography, sgs_wf_part, sort_count,              &
    165171                total_number_of_particles,                                     &
    166172                use_sgs_for_particles,                                         &
     
    190196    PUBLIC lpm_init, lpm_create_particle
    191197
    192 CONTAINS
     198 CONTAINS
    193199
    194200!------------------------------------------------------------------------------!
     
    290296
    291297!
    292 !-- Allocate arrays required for calculating particle SGS velocities
     298!-- Allocate arrays required for calculating particle SGS velocities.
     299!-- Initialize prefactor required for stoachastic Weil equation.
    293300    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    294301       ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    295302                 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    296303                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     304
     305       sgs_wf_part = 1.0_wp / 3.0_wp   
    297306    ENDIF
    298307
     
    334343!
    335344!--    Precalculate LOG(z/z0)
    336        height_p    = 0.0_wp
     345       height_p    = z0_av_global
    337346       DO  k = 1, number_of_sublayers
    338347
     
    341350
    342351       ENDDO
    343 
    344352
    345353    ENDIF
     
    444452                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
    445453                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0, 0, &
    446                                       0, .FALSE., -1)
     454                                      0, .FALSE., -1 )
    447455
    448456       particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
     
    477485
    478486       CALL lpm_create_particle (PHASE_INIT)
    479 
    480487!
    481488!--    User modification of initial particles
     
    499506    number_of_particles = prt_count(nzb+1,nys,nxl)
    500507    particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
    501 
    502508!
    503509!-- Formats
     
    520526
    521527    USE particle_attributes,                                                   &
    522         ONLY: monodisperse_aerosols
     528        ONLY: deleted_particles, monodisperse_aerosols
    523529
    524530    IMPLICIT  NONE
    525531
    526     INTEGER(iwp)               ::  alloc_size  !<
    527     INTEGER(iwp)               ::  i           !<
    528     INTEGER(iwp)               ::  ip          !<
    529     INTEGER(iwp)               ::  j           !<
    530     INTEGER(iwp)               ::  jp          !<
    531     INTEGER(iwp)               ::  kp          !<
    532     INTEGER(iwp)               ::  loop_stride !<
    533     INTEGER(iwp)               ::  n           !<
    534     INTEGER(iwp)               ::  new_size    !<
    535 
    536     INTEGER(iwp), INTENT(IN)   ::  phase       !<
    537 
    538     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_count !<
    539     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_start !<
    540 
    541     LOGICAL                    ::  first_stride !<
    542 
    543     REAL(wp)                   ::  pos_x !<
    544     REAL(wp)                   ::  pos_y !<
    545     REAL(wp)                   ::  pos_z !<
    546 
    547     TYPE(particle_type),TARGET ::  tmp_particle !<
     532    INTEGER(iwp)               ::  alloc_size  !< relative increase of allocated memory for particles
     533    INTEGER(iwp)               ::  i           !< loop variable ( particle groups )
     534    INTEGER(iwp)               ::  ip          !< index variable along x
     535    INTEGER(iwp)               ::  j           !< loop variable ( particles per point )
     536    INTEGER(iwp)               ::  jp          !< index variable along y
     537    INTEGER(iwp)               ::  kp          !< index variable along z
     538    INTEGER(iwp)               ::  loop_stride !< loop variable for initialization
     539    INTEGER(iwp)               ::  n           !< loop variable ( number of particles )
     540    INTEGER(iwp)               ::  new_size    !< new size of allocated memory for particles
     541
     542    INTEGER(iwp), INTENT(IN)   ::  phase       !< mode of inititialization
     543
     544    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_count !< start address of new particle
     545    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_start !< start address of new particle
     546
     547    LOGICAL                    ::  first_stride !< flag for initialization
     548
     549    REAL(wp)                   ::  pos_x      !< increment for particle position in x     
     550    REAL(wp)                   ::  pos_y      !< increment for particle position in y 
     551    REAL(wp)                   ::  pos_z      !< increment for particle position in z     
     552    REAL(wp)                   ::  rand_contr !< dummy argument for random position
     553
     554    TYPE(particle_type),TARGET ::  tmp_particle !< temporary particle used for initialization
    548555
    549556!
     
    618625                            tmp_particle%tail_id       = 0     !unused
    619626
     627
    620628!
    621629!--                         Determine the grid indices of the particle position
     
    628636!--                            Particle height is given relative to topography
    629637                               kp = kp + nzb_w_inner(jp,ip)
    630                                tmp_particle%z = tmp_particle%z + zw(kp)
     638                               tmp_particle%z = tmp_particle%z +               &
     639                                                         zw(nzb_w_inner(jp,ip))
    631640                               IF ( kp > nzt )  THEN
    632641                                  pos_x = pos_x + pdx(i)
    633642                                  CYCLE xloop
    634643                               ENDIF
     644                            ELSEIF ( .NOT. seed_follows_topography .AND.       &
     645                                     tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
     646                               pos_x = pos_x + pdx(i)
     647                               CYCLE xloop                               
    635648                            ENDIF
    636649
     
    644657                               ENDIF
    645658                               grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
     659
    646660                            ENDIF
    647661                         ENDDO
     
    691705                      ENDIF
    692706                   ENDIF
     707
    693708                ENDDO
    694709             ENDDO
    695710          ENDDO
    696711       ENDIF
     712
    697713    ENDDO
    698714
     
    707723
    708724!
    709 !-- Add random fluctuation to particle positions
     725!-- Add random fluctuation to particle positions.
    710726    IF ( random_start_position )  THEN
    711727       DO  ip = nxl, nxr
     
    715731                IF ( number_of_particles <= 0 )  CYCLE
    716732                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    717 
    718                 DO  n = local_start(kp,jp,ip), number_of_particles              !Move only new particles
     733!
     734!--             Move only new particles. Moreover, limit random fluctuation
     735!--             in order to prevent that particles move more than one grid box,
     736!--             which would lead to problems concerning particle exchange
     737!--             between processors in case pdx/pdy are larger than dx/dy,
     738!--             respectively. 
     739                DO  n = local_start(kp,jp,ip), number_of_particles
    719740                   IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
     741                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     742                                     pdx(particles(n)%group)
    720743                      particles(n)%x = particles(n)%x +                        &
    721                                    ( random_function( iran_part ) - 0.5_wp ) * &
    722                                    pdx(particles(n)%group)
     744                              MERGE( rand_contr, SIGN( dx, rand_contr ), &
     745                                     ABS( rand_contr ) < dx                    &
     746                                   )
    723747                   ENDIF
    724748                   IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
     749                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     750                                     pdy(particles(n)%group)
    725751                      particles(n)%y = particles(n)%y +                        &
    726                                    ( random_function( iran_part ) - 0.5_wp ) * &
    727                                    pdy(particles(n)%group)
     752                              MERGE( rand_contr, SIGN( dy, rand_contr ), &
     753                                     ABS( rand_contr ) < dy                    &
     754                                   )
    728755                   ENDIF
    729756                   IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
     757                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     758                                     pdz(particles(n)%group)
    730759                      particles(n)%z = particles(n)%z +                        &
    731                                    ( random_function( iran_part ) - 0.5_wp ) * &
    732                                    pdz(particles(n)%group)
     760                              MERGE( rand_contr, SIGN( dz, rand_contr ), &
     761                                     ABS( rand_contr ) < dz                    &
     762                                   )
    733763                   ENDIF
    734764                ENDDO
    735765!
    736 !--             Identify particles located outside the model domain
     766!--             Identify particles located outside the model domain and reflect
     767!--             or absorb them if necessary.
    737768                CALL lpm_boundary_conds( 'bottom/top' )
     769!
     770!--             Furthermore, remove particles located in topography. Note, as
     771!--             the particle speed is still zero at this point, wall
     772!--             reflection boundary conditions will not work in this case.
     773                particles =>                                                   &
     774                       grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     775                DO  n = local_start(kp,jp,ip), number_of_particles
     776                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
     777                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
     778                   IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) )  THEN
     779                      particles(n)%particle_mask = .FALSE.
     780                      deleted_particles = deleted_particles + 1
     781                   ENDIF
     782                ENDDO
    738783             ENDDO
    739784          ENDDO
  • palm/trunk/SOURCE/lpm_init_sgs_tke.f90

    r1823 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! sgs_wfu_par, sgs_wfv_par and sgs_wfw_par are replaced by sgs_wf_par
    2222!
    2323! Former revisions:
     
    7272
    7373    USE particle_attributes,                                                   &
    74         ONLY:  sgs_wfu_part, sgs_wfv_part, sgs_wfw_part
     74        ONLY:  sgs_wf_part
    7575
    7676    USE pegrid
     
    9595                  k  > nzb_s_inner(j,i+1) )                                    &
    9696             THEN
    97                 de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
     97                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
    9898                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
    9999             ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)     &
    100100                      .AND.  k <= nzb_s_inner(j,i+1) )                         &
    101101             THEN
    102                 de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
     102                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
    103103                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
    104104             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) )    &
     
    109109                de_dx(k,j,i) = 0.0_wp
    110110             ELSE
    111                 de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
     111                de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
    112112             ENDIF
    113113
     
    115115                  k  > nzb_s_inner(j+1,i) )                                    &
    116116             THEN
    117                 de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
     117                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
    118118                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
    119119             ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  k  > nzb_s_inner(j,i)    &
    120120                      .AND.  k <= nzb_s_inner(j+1,i) )                         &
    121121             THEN
    122                 de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
     122                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
    123123                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
    124124             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) )    &
     
    129129                de_dy(k,j,i) = 0.0_wp
    130130             ELSE
    131                 de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
     131                de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
    132132             ENDIF
    133133
     
    142142
    143143          DO  k = nzb_s_inner(j,i)+2, nzt-1
    144              de_dz(k,j,i)  = 2.0_wp * sgs_wfw_part *                           &
     144             de_dz(k,j,i)  = 2.0_wp * sgs_wf_part *                            &
    145145                             ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
    146146          ENDDO
     
    148148          k = nzb_s_inner(j,i)
    149149          de_dz(nzb:k,j,i) = 0.0_wp
    150           de_dz(k+1,j,i)   = 2.0_wp * sgs_wfw_part *                           &
     150          de_dz(k+1,j,i)   = 2.0_wp * sgs_wf_part *                            &
    151151                             ( e(k+2,j,i) - e(k+1,j,i) ) / ( zu(k+2) - zu(k+1) )
    152152          de_dz(nzt,j,i)   = 0.0_wp
  • palm/trunk/SOURCE/mod_particle_attributes.f90

    r1872 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! -sgs_wfu_par, sgs_wfv_par, sgs_wfw_par
     22! + sgs_wf_par
    2223!
    2324! Former revisions:
     
    109110                 initial_weighting_factor = 1.0_wp,                            &
    110111                 particle_advection_start = 0.0_wp,                            &
    111                  sgs_wfu_part = 0.3333333_wp, sgs_wfv_part = 0.3333333_wp,     &
    112                  sgs_wfw_part = 0.3333333_wp, time_prel = 0.0_wp,              &
    113                  time_sort_particles = 0.0_wp,                                 &
     112                 sgs_wf_part, time_prel = 0.0_wp, time_sort_particles = 0.0_wp,&
    114113                 time_write_particle_data = 0.0_wp, z0_av_global
    115114
  • palm/trunk/SOURCE/pres.f90

    r1919 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Bugfix: weight_substep for initial call, replace by local variable
    2222!
    2323! Former revisions:
     
    162162    REAL(wp)     ::  threadsum      !<
    163163    REAL(wp)     ::  weight_pres_l  !<
     164    REAL(wp)     ::  weight_substep_l !<
    164165
    165166    REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
     
    178179!
    179180!--    If pres is called before initial time step
    180        weight_pres_l = 1.0_wp
    181        d_weight_pres = 1.0_wp
     181       weight_pres_l    = 1.0_wp
     182       d_weight_pres    = 1.0_wp
     183       weight_substep_l = 1.0_wp
    182184    ELSE
    183        weight_pres_l = weight_pres(intermediate_timestep_count)
    184        d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
     185       weight_pres_l    = weight_pres(intermediate_timestep_count)
     186       d_weight_pres    = 1.0_wp / weight_pres(intermediate_timestep_count)
     187       weight_substep_l = weight_substep(intermediate_timestep_count)
    185188    ENDIF
    186189
     
    591594       !$OMP PARALLEL PRIVATE (i,j,k)
    592595       !$OMP DO
    593        !$acc kernels present( p, tend, weight_substep )
     596       !$acc kernels present( p, tend, weight_substep_l )
    594597       !$acc loop independent
    595598       DO  i = nxl-1, nxr+1
     
    599602             DO  k = nzb, nzt+1
    600603                p(k,j,i) = tend(k,j,i) * &
    601                            weight_substep(intermediate_timestep_count)
     604                           weight_substep_l
    602605             ENDDO
    603606          ENDDO
     
    609612       !$OMP PARALLEL PRIVATE (i,j,k)
    610613       !$OMP DO
    611        !$acc kernels present( p, tend, weight_substep )
     614       !$acc kernels present( p, tend, weight_substep_l )
    612615       !$acc loop independent
    613616       DO  i = nxl-1, nxr+1
     
    617620             DO  k = nzb, nzt+1
    618621                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
    619                            weight_substep(intermediate_timestep_count)
     622                           weight_substep_l
    620623             ENDDO
    621624          ENDDO
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r1921 r1929  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Bugfix: avoid segmentation fault in case one grid point is horizontally
     22! completely surrounded by topography
    2223!
    2324! Former revisions:
     
    726727!--             consequence would result in very large shear stresses and very
    727728!--             small momentum fluxes (both are generally unrealistic).
    728                 IF ( ( z_mo / ol(j,i) ) < zeta_min )  ol(j,i) = z_mo / zeta_min
    729                 IF ( ( z_mo / ol(j,i) ) > zeta_max )  ol(j,i) = z_mo / zeta_max
     729                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min )            &
     730                   ol(j,i) = z_mo / zeta_min
     731                IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max )            &
     732                   ol(j,i) = z_mo / zeta_max
    730733             ENDDO
    731734          ENDDO
  • palm/trunk/SOURCE/wind_turbine_model_mod.f90

    r1917 r1929  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix: added preprocessor directives for parallel and serial mode
    2222!
    2323! Former revisions:
     
    20582058
    20592059!
    2060 !--          Transfer of information to the other cpus           
     2060!--          Transfer of information to the other cpus
     2061#if defined( __parallel )         
    20612062             CALL MPI_ALLREDUCE( u_inflow_l, u_inflow, nturbines, MPI_REAL,    &
    20622063                                 MPI_SUM, comm2d, ierr )
     
    20652066             CALL MPI_ALLREDUCE( phi_yaw_l, phi_yaw, nturbines, MPI_REAL,      &
    20662067                                 MPI_SUM, comm2d, ierr )
    2067              
     2068#else
     2069             u_inflow = u_inflow_l
     2070             wdir     = wdir_l
     2071             phi_yaw  = phi_yaw_l
     2072#endif
    20682073             DO inot = 1, nturbines
    20692074!             
     
    20802085!              CALL MPI_ALLREDUCE( omega_gen, omega_gen_old, nturbines,        &
    20812086!                                  MPI_REAL,MPI_SUM, comm2d, ierr )
     2087#if defined( __parallel )   
    20822088             CALL MPI_ALLREDUCE( torque_gen, torque_gen_old, nturbines,        &
    20832089                                 MPI_REAL, MPI_SUM, comm2d, ierr )
     
    20862092             CALL MPI_ALLREDUCE( omega_gen_f, omega_gen_f_old, nturbines,      &
    20872093                                 MPI_REAL, MPI_SUM, comm2d, ierr )
     2094#else
     2095             torque_gen_old  = torque_gen
     2096             omega_rot       = omega_rot_l
     2097             omega_gen_f_old = omega_gen_f
     2098#endif
    20882099           
    20892100          ENDIF
Note: See TracChangeset for help on using the changeset viewer.