Changeset 2417 for palm/trunk/SOURCE


Ignore:
Timestamp:
Sep 6, 2017 3:22:27 PM (7 years ago)
Author:
suehring
Message:

Major bugfix in modeling SGS particle speeds.

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r2263 r2417  
    2525! -----------------
    2626! $Id$
     27! Major bugfixes in modeling SGS particle speeds (since revision 1359).
     28! Particle sorting added to distinguish between already completed and
     29! non-completed particles.
     30!
     31! 2263 2017-06-08 14:59:01Z schwenkel
    2732! Implemented splitting and merging algorithm
    2833!
     
    140145
    141146    USE lpm_pack_arrays_mod,                                                   &
    142         ONLY:  lpm_pack_all_arrays
     147        ONLY:  lpm_pack_all_arrays, lpm_sort
    143148
    144149    USE particle_attributes,                                                   &
     
    259264       CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
    260265       
    261        grid_particles(:,:,:)%time_loop_done = .TRUE.
    262266!
    263267!--    If particle advection includes SGS velocity components, calculate the
     
    265269!--    horizontally averaged profiles of the SGS TKE and the resolved-scale
    266270!--    velocity variances)
    267 
    268271       IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    269272          CALL lpm_init_sgs_tke
    270273       ENDIF
     274
     275!
     276!--    In case SGS-particle speed is considered, particles may carry out
     277!--    several particle timesteps. In order to prevent unnecessary
     278!--    treatment of particles that already reached the final time level,
     279!--    particles are sorted into contiguous blocks of finished and
     280!--    not-finished particles, in addition to their already sorting
     281!--    according to their sub-boxes.
     282       IF ( .NOT. first_loop_stride  .AND.  use_sgs_for_particles )            &
     283          CALL lpm_sort
    271284
    272285       DO  i = nxl, nxr
     
    323336!--             optimization is still possible.)
    324337                IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
    325                    CALL lpm_boundary_conds( 'walls' )
     338                   CALL lpm_boundary_conds( 'walls' )!, i, j, k )
    326339                ENDIF
    327340!
     
    333346!--             the top or bottom boundary and delete those particles, which are
    334347!--             older than allowed
    335                 CALL lpm_boundary_conds( 'bottom/top' )
     348                CALL lpm_boundary_conds( 'bottom/top' ) !, i, j, k )
    336349!
    337350!---            If not all particles of the actual grid cell have reached the
    338 !--             LES timestep, this cell has to to another loop iteration. Due to
    339 !--             the fact that particles can move into neighboring grid cell,
    340 !--             these neighbor cells also have to perform another loop iteration
     351!--             LES timestep, this cell has to do another loop iteration. Due to
     352!--             the fact that particles can move into neighboring grid cells,
     353!--             these neighbor cells also have to perform another loop iteration.
     354!--             Please note, this realization does not work properly if
     355!--             particles move into another subdomain.
    341356                IF ( .NOT. dt_3d_reached_l )  THEN
    342                    ks = MAX(nzb+1,k)
    343                    ke = MIN(nzt,k)
    344                    js = MAX(nys,j)
    345                    je = MIN(nyn,j)
    346                    is = MAX(nxl,i)
    347                    ie = MIN(nxr,i)
     357                   ks = MAX(nzb+1,k-1)
     358                   ke = MIN(nzt,k+1)
     359                   js = MAX(nys,j-1)
     360                   je = MIN(nyn,j+1)
     361                   is = MAX(nxl,i-1)
     362                   ie = MIN(nxr,i+1)
    348363                   grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE.
     364                ELSE
     365                   grid_particles(k,j,i)%time_loop_done = .TRUE.
    349366                ENDIF
    350367
     
    392409!--    determine new number of particles
    393410       CALL lpm_pack_all_arrays
    394 
    395411!
    396412!--    Initialize variables for the next (sub-) timestep, i.e., for marking
  • palm/trunk/SOURCE/lpm_advec.f90

    r2318 r2417  
    2525! -----------------
    2626! $Id$
     27! Particle loops adapted for sub-box structure, i.e. for each sub-box the
     28! particle loop runs from start_index up to end_index instead from 1 to
     29! number_of_particles. This way, it is possible to skip unnecessary
     30! computations for particles that already completed the LES timestep.
     31!
     32! 2318 2017-07-20 17:27:44Z suehring
    2733! Get topography top index via Function call
    2834!
     
    267273
    268274    DO  nb = 0, 7
    269 
     275       
    270276       i = ip
    271277       j = jp + block_offset(nb)%j_off
     
    287293!
    288294!--       Determine vertical index of topography top
    289           k_wall = get_topography_top_index( jlog,ilog, 's' )
     295          k_wall = get_topography_top_index( jlog, ilog, 's' )
    290296
    291297          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     
    648654       ELSE ! non-flat topography, e.g., buildings
    649655
    650           DO  n = 1, number_of_particles
    651              i = particles(n)%x * ddx
    652              j = particles(n)%y * ddy
    653              k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
    654                  + offset_ocean_nzt                      ! only exact if eq.dist
    655 !
    656 !--          In case that there are buildings it has to be determined
    657 !--          how many of the gridpoints defining the particle box are
    658 !--          situated within a building
    659 !--          gp_outside_of_building(1): i,j,k
    660 !--          gp_outside_of_building(2): i,j+1,k
    661 !--          gp_outside_of_building(3): i,j,k+1
    662 !--          gp_outside_of_building(4): i,j+1,k+1
    663 !--          gp_outside_of_building(5): i+1,j,k
    664 !--          gp_outside_of_building(6): i+1,j+1,k
    665 !--          gp_outside_of_building(7): i+1,j,k+1
    666 !--          gp_outside_of_building(8): i+1,j+1,k+1
    667 
    668              gp_outside_of_building = 0
    669              location = 0.0_wp
    670              num_gp = 0
    671 
    672 !
    673 !--          Determine vertical index of topography top at (j,i)
    674              k_wall = get_topography_top_index( j, i, 's' )
    675 !
    676 !--          To do: Reconsider order of computations in order to avoid
    677 !--          unnecessary index calculations.
    678              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    679                 num_gp = num_gp + 1
    680                 gp_outside_of_building(1) = 1
    681                 location(num_gp,1) = i * dx
    682                 location(num_gp,2) = j * dy
    683                 location(num_gp,3) = k * dz - 0.5_wp * dz
    684                 ei(num_gp)     = e(k,j,i)
    685                 dissi(num_gp)  = diss(k,j,i)
    686                 de_dxi(num_gp) = de_dx(k,j,i)
    687                 de_dyi(num_gp) = de_dy(k,j,i)
    688                 de_dzi(num_gp) = de_dz(k,j,i)
    689              ENDIF
    690 
    691 !
    692 !--          Determine vertical index of topography top at (j+1,i)
    693              k_wall = get_topography_top_index( j+1, i, 's' )
    694 
    695              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    696                 num_gp = num_gp + 1
    697                 gp_outside_of_building(2) = 1
    698                 location(num_gp,1) = i * dx
    699                 location(num_gp,2) = (j+1) * dy
    700                 location(num_gp,3) = k * dz - 0.5_wp * dz
    701                 ei(num_gp)     = e(k,j+1,i)
    702                 dissi(num_gp)  = diss(k,j+1,i)
    703                 de_dxi(num_gp) = de_dx(k,j+1,i)
    704                 de_dyi(num_gp) = de_dy(k,j+1,i)
    705                 de_dzi(num_gp) = de_dz(k,j+1,i)
    706              ENDIF
    707 
    708 !
    709 !--          Determine vertical index of topography top at (j,i)
    710              k_wall = get_topography_top_index( j, i, 's' )
    711 
    712              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    713                 num_gp = num_gp + 1
    714                 gp_outside_of_building(3) = 1
    715                 location(num_gp,1) = i * dx
    716                 location(num_gp,2) = j * dy
    717                 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    718                 ei(num_gp)     = e(k+1,j,i)
    719                 dissi(num_gp)  = diss(k+1,j,i)
    720                 de_dxi(num_gp) = de_dx(k+1,j,i)
    721                 de_dyi(num_gp) = de_dy(k+1,j,i)
    722                 de_dzi(num_gp) = de_dz(k+1,j,i)
    723              ENDIF
    724 
    725 !
    726 !--          Determine vertical index of topography top at (j+1,i)
    727              k_wall = get_topography_top_index( j+1, i, 's' )
    728              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    729                 num_gp = num_gp + 1
    730                 gp_outside_of_building(4) = 1
    731                 location(num_gp,1) = i * dx
    732                 location(num_gp,2) = (j+1) * dy
    733                 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    734                 ei(num_gp)     = e(k+1,j+1,i)
    735                 dissi(num_gp)  = diss(k+1,j+1,i)
    736                 de_dxi(num_gp) = de_dx(k+1,j+1,i)
    737                 de_dyi(num_gp) = de_dy(k+1,j+1,i)
    738                 de_dzi(num_gp) = de_dz(k+1,j+1,i)
    739              ENDIF
    740 
    741 !
    742 !--          Determine vertical index of topography top at (j,i+1)
    743              k_wall = get_topography_top_index( j, i+1, 's' )
    744              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    745                 num_gp = num_gp + 1
    746                 gp_outside_of_building(5) = 1
    747                 location(num_gp,1) = (i+1) * dx
    748                 location(num_gp,2) = j * dy
    749                 location(num_gp,3) = k * dz - 0.5_wp * dz
    750                 ei(num_gp)     = e(k,j,i+1)
    751                 dissi(num_gp)  = diss(k,j,i+1)
    752                 de_dxi(num_gp) = de_dx(k,j,i+1)
    753                 de_dyi(num_gp) = de_dy(k,j,i+1)
    754                 de_dzi(num_gp) = de_dz(k,j,i+1)
    755              ENDIF
    756 
    757 !
    758 !--          Determine vertical index of topography top at (j+1,i+1)
    759              k_wall = get_topography_top_index( j+1, i+1, 's' )
    760 
    761              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    762                 num_gp = num_gp + 1
    763                 gp_outside_of_building(6) = 1
    764                 location(num_gp,1) = (i+1) * dx
    765                 location(num_gp,2) = (j+1) * dy
    766                 location(num_gp,3) = k * dz - 0.5_wp * dz
    767                 ei(num_gp)     = e(k,j+1,i+1)
    768                 dissi(num_gp)  = diss(k,j+1,i+1)
    769                 de_dxi(num_gp) = de_dx(k,j+1,i+1)
    770                 de_dyi(num_gp) = de_dy(k,j+1,i+1)
    771                 de_dzi(num_gp) = de_dz(k,j+1,i+1)
    772              ENDIF
    773 
    774 !
    775 !--          Determine vertical index of topography top at (j,i+1)
    776              k_wall = get_topography_top_index( j, i+1, 's' )
    777 
    778              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    779                 num_gp = num_gp + 1
    780                 gp_outside_of_building(7) = 1
    781                 location(num_gp,1) = (i+1) * dx
    782                 location(num_gp,2) = j * dy
    783                 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    784                 ei(num_gp)     = e(k+1,j,i+1)
    785                 dissi(num_gp)  = diss(k+1,j,i+1)
    786                 de_dxi(num_gp) = de_dx(k+1,j,i+1)
    787                 de_dyi(num_gp) = de_dy(k+1,j,i+1)
    788                 de_dzi(num_gp) = de_dz(k+1,j,i+1)
    789              ENDIF
    790 
    791 !
    792 !--          Determine vertical index of topography top at (j+1,i+1)
    793              k_wall = get_topography_top_index( j+1, i+1, 's' )
    794 
    795              IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
    796                 num_gp = num_gp + 1
    797                 gp_outside_of_building(8) = 1
    798                 location(num_gp,1) = (i+1) * dx
    799                 location(num_gp,2) = (j+1) * dy
    800                 location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    801                 ei(num_gp)     = e(k+1,j+1,i+1)
    802                 dissi(num_gp)  = diss(k+1,j+1,i+1)
    803                 de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    804                 de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    805                 de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    806              ENDIF
    807 !
    808 !--          If all gridpoints are situated outside of a building, then the
    809 !--          ordinary interpolation scheme can be used.
    810              IF ( num_gp == 8 )  THEN
    811 
    812                 x  = particles(n)%x - i * dx
    813                 y  = particles(n)%y - j * dy
    814                 aa = x**2          + y**2
    815                 bb = ( dx - x )**2 + y**2
    816                 cc = x**2          + ( dy - y )**2
    817                 dd = ( dx - x )**2 + ( dy - y )**2
    818                 gg = aa + bb + cc + dd
    819  
    820                 e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
    821                           + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
    822                           ) / ( 3.0_wp * gg )
    823    
    824                 IF ( k == nzt )  THEN
    825                    e_int(n) = e_int_l
    826                 ELSE
    827                    e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    828                                ( gg - bb ) * e(k+1,j,i+1) + &
    829                                ( gg - cc ) * e(k+1,j+1,i) + &
    830                                ( gg - dd ) * e(k+1,j+1,i+1) &
    831                              ) / ( 3.0_wp * gg )
    832                    e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
    833                                        ( e_int_u - e_int_l )
    834                 ENDIF
    835 !
    836 !--             Needed to avoid NaN particle velocities (this might not be
    837 !--             required any more)
    838                 IF ( e_int(n) <= 0.0_wp )  THEN
    839                    e_int(n) = 1.0E-20_wp
    840                 ENDIF
    841 !
    842 !--             Interpolate the TKE gradient along x (adopt incides i,j,k
    843 !--             and all position variables from above (TKE))
    844                 de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    845                                 ( gg - bb ) * de_dx(k,j,i+1) + &
    846                                 ( gg - cc ) * de_dx(k,j+1,i) + &
    847                                 ( gg - dd ) * de_dx(k,j+1,i+1) &
    848                               ) / ( 3.0_wp * gg )
    849 
    850                 IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
    851                    de_dx_int(n) = de_dx_int_l
    852                 ELSE
    853                    de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    854                                    ( gg - bb ) * de_dx(k+1,j,i+1) + &
    855                                    ( gg - cc ) * de_dx(k+1,j+1,i) + &
    856                                    ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    857                                  ) / ( 3.0_wp * gg )
    858                    de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
    859                                            dz * ( de_dx_int_u - de_dx_int_l )
    860                 ENDIF
    861 
    862 !
    863 !--             Interpolate the TKE gradient along y
    864                 de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    865                                 ( gg - bb ) * de_dy(k,j,i+1) + &
    866                                 ( gg - cc ) * de_dy(k,j+1,i) + &
    867                                 ( gg - dd ) * de_dy(k,j+1,i+1) &
    868                               ) / ( 3.0_wp * gg )
    869                 IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    870                    de_dy_int(n) = de_dy_int_l
    871                 ELSE
    872                    de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    873                                    ( gg - bb ) * de_dy(k+1,j,i+1) + &
    874                                    ( gg - cc ) * de_dy(k+1,j+1,i) + &
    875                                    ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    876                                  ) / ( 3.0_wp * gg )
    877                    de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
    878                                            dz * ( de_dy_int_u - de_dy_int_l )
    879                 ENDIF
    880 
    881 !
    882 !--             Interpolate the TKE gradient along z
    883                 IF ( zv(n) < 0.5_wp * dz )  THEN
    884                    de_dz_int(n) = 0.0_wp
    885                 ELSE
    886                    de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    887                                    ( gg - bb ) * de_dz(k,j,i+1) + &
    888                                    ( gg - cc ) * de_dz(k,j+1,i) + &
    889                                    ( gg - dd ) * de_dz(k,j+1,i+1) &
    890                                  ) / ( 3.0_wp * gg )
    891  
    892                    IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    893                       de_dz_int(n) = de_dz_int_l
    894                    ELSE
    895                       de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    896                                       ( gg - bb ) * de_dz(k+1,j,i+1) + &
    897                                       ( gg - cc ) * de_dz(k+1,j+1,i) + &
    898                                       ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    899                                     ) / ( 3.0_wp * gg )
    900                       de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
    901                                            dz * ( de_dz_int_u - de_dz_int_l )
    902                    ENDIF
    903                 ENDIF
    904 
    905 !
    906 !--             Interpolate the dissipation of TKE
    907                 diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    908                                ( gg - bb ) * diss(k,j,i+1) + &
    909                                ( gg - cc ) * diss(k,j+1,i) + &
    910                                ( gg - dd ) * diss(k,j+1,i+1) &
    911                              ) / ( 3.0_wp * gg )
    912  
    913                 IF ( k == nzt )  THEN
    914                    diss_int(n) = diss_int_l
    915                 ELSE
    916                    diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    917                                   ( gg - bb ) * diss(k+1,j,i+1) + &
    918                                   ( gg - cc ) * diss(k+1,j+1,i) + &
    919                                   ( gg - dd ) * diss(k+1,j+1,i+1) &
    920                                 ) / ( 3.0_wp * gg )
    921                    diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
    922                                            ( diss_int_u - diss_int_l )
    923                 ENDIF
    924 !
    925 !--             Set flag for stochastic equation.
    926                 term_1_2(n) = 1.0_wp
    927  
    928              ELSE
    929  
    930 !
    931 !--             If wall between gridpoint 1 and gridpoint 5, then
    932 !--             Neumann boundary condition has to be applied
    933                 IF ( gp_outside_of_building(1) == 1  .AND. &
    934                      gp_outside_of_building(5) == 0 )  THEN
     656          DO  nb = 0, 7
     657
     658             DO  n = start_index(nb), end_index(nb)
     659
     660                i = particles(n)%x * ddx
     661                j = particles(n)%y * ddy
     662                k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz  &
     663                    + offset_ocean_nzt                      ! only exact if eq.dist
     664!
     665!--             In case that there are buildings it has to be determined
     666!--             how many of the gridpoints defining the particle box are
     667!--             situated within a building
     668!--             gp_outside_of_building(1): i,j,k
     669!--             gp_outside_of_building(2): i,j+1,k
     670!--             gp_outside_of_building(3): i,j,k+1
     671!--             gp_outside_of_building(4): i,j+1,k+1
     672!--             gp_outside_of_building(5): i+1,j,k
     673!--             gp_outside_of_building(6): i+1,j+1,k
     674!--             gp_outside_of_building(7): i+1,j,k+1
     675!--             gp_outside_of_building(8): i+1,j+1,k+1
     676
     677                gp_outside_of_building = 0
     678                location = 0.0_wp
     679                num_gp = 0
     680
     681!
     682!--             Determine vertical index of topography top at (j,i)
     683                k_wall = get_topography_top_index( j, i, 's' )
     684!
     685!--             To do: Reconsider order of computations in order to avoid
     686!--             unnecessary index calculations.
     687                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    935688                   num_gp = num_gp + 1
    936                    location(num_gp,1) = i * dx + 0.5_wp * dx
     689                   gp_outside_of_building(1) = 1
     690                   location(num_gp,1) = i * dx
    937691                   location(num_gp,2) = j * dy
    938692                   location(num_gp,3) = k * dz - 0.5_wp * dz
    939693                   ei(num_gp)     = e(k,j,i)
    940694                   dissi(num_gp)  = diss(k,j,i)
    941                    de_dxi(num_gp) = 0.0_wp
     695                   de_dxi(num_gp) = de_dx(k,j,i)
    942696                   de_dyi(num_gp) = de_dy(k,j,i)
    943697                   de_dzi(num_gp) = de_dz(k,j,i)
    944698                ENDIF
    945699
    946                 IF ( gp_outside_of_building(5) == 1  .AND. &
    947                      gp_outside_of_building(1) == 0 )  THEN
     700!
     701!--             Determine vertical index of topography top at (j+1,i)
     702                k_wall = get_topography_top_index( j+1, i, 's' )
     703
     704                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    948705                   num_gp = num_gp + 1
    949                    location(num_gp,1) = i * dx + 0.5_wp * dx
     706                   gp_outside_of_building(2) = 1
     707                   location(num_gp,1) = i * dx
     708                   location(num_gp,2) = (j+1) * dy
     709                   location(num_gp,3) = k * dz - 0.5_wp * dz
     710                   ei(num_gp)     = e(k,j+1,i)
     711                   dissi(num_gp)  = diss(k,j+1,i)
     712                   de_dxi(num_gp) = de_dx(k,j+1,i)
     713                   de_dyi(num_gp) = de_dy(k,j+1,i)
     714                   de_dzi(num_gp) = de_dz(k,j+1,i)
     715                ENDIF
     716
     717!
     718!--             Determine vertical index of topography top at (j,i)
     719                k_wall = get_topography_top_index( j, i, 's' )
     720
     721                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
     722                   num_gp = num_gp + 1
     723                   gp_outside_of_building(3) = 1
     724                   location(num_gp,1) = i * dx
     725                   location(num_gp,2) = j * dy
     726                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
     727                   ei(num_gp)     = e(k+1,j,i)
     728                   dissi(num_gp)  = diss(k+1,j,i)
     729                   de_dxi(num_gp) = de_dx(k+1,j,i)
     730                   de_dyi(num_gp) = de_dy(k+1,j,i)
     731                   de_dzi(num_gp) = de_dz(k+1,j,i)
     732                ENDIF
     733
     734!
     735!--             Determine vertical index of topography top at (j+1,i)
     736                k_wall = get_topography_top_index( j+1, i, 's' )
     737                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
     738                   num_gp = num_gp + 1
     739                   gp_outside_of_building(4) = 1
     740                   location(num_gp,1) = i * dx
     741                   location(num_gp,2) = (j+1) * dy
     742                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
     743                   ei(num_gp)     = e(k+1,j+1,i)
     744                   dissi(num_gp)  = diss(k+1,j+1,i)
     745                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
     746                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
     747                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
     748                ENDIF
     749
     750!
     751!--             Determine vertical index of topography top at (j,i+1)
     752                k_wall = get_topography_top_index( j, i+1, 's' )
     753                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
     754                   num_gp = num_gp + 1
     755                   gp_outside_of_building(5) = 1
     756                   location(num_gp,1) = (i+1) * dx
    950757                   location(num_gp,2) = j * dy
    951758                   location(num_gp,3) = k * dz - 0.5_wp * dz
    952759                   ei(num_gp)     = e(k,j,i+1)
    953760                   dissi(num_gp)  = diss(k,j,i+1)
    954                    de_dxi(num_gp) = 0.0_wp
     761                   de_dxi(num_gp) = de_dx(k,j,i+1)
    955762                   de_dyi(num_gp) = de_dy(k,j,i+1)
    956763                   de_dzi(num_gp) = de_dz(k,j,i+1)
     
    958765
    959766!
    960 !--             If wall between gridpoint 5 and gridpoint 6, then
    961 !--             then Neumann boundary condition has to be applied
    962                 IF ( gp_outside_of_building(5) == 1  .AND. &
    963                      gp_outside_of_building(6) == 0 )  THEN
     767!--             Determine vertical index of topography top at (j+1,i+1)
     768                k_wall = get_topography_top_index( j+1, i+1, 's' )
     769
     770                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    964771                   num_gp = num_gp + 1
     772                   gp_outside_of_building(6) = 1
    965773                   location(num_gp,1) = (i+1) * dx
    966                    location(num_gp,2) = j * dy + 0.5_wp * dy
    967                    location(num_gp,3) = k * dz - 0.5_wp * dz
    968                    ei(num_gp)     = e(k,j,i+1)
    969                    dissi(num_gp)  = diss(k,j,i+1)
    970                    de_dxi(num_gp) = de_dx(k,j,i+1)
    971                    de_dyi(num_gp) = 0.0_wp
    972                    de_dzi(num_gp) = de_dz(k,j,i+1)
    973                 ENDIF
    974 
    975                 IF ( gp_outside_of_building(6) == 1  .AND. &
    976                      gp_outside_of_building(5) == 0 )  THEN
    977                    num_gp = num_gp + 1
    978                    location(num_gp,1) = (i+1) * dx
    979                    location(num_gp,2) = j * dy + 0.5_wp * dy
     774                   location(num_gp,2) = (j+1) * dy
    980775                   location(num_gp,3) = k * dz - 0.5_wp * dz
    981776                   ei(num_gp)     = e(k,j+1,i+1)
    982777                   dissi(num_gp)  = diss(k,j+1,i+1)
    983778                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
    984                    de_dyi(num_gp) = 0.0_wp
    985                    de_dzi(num_gp) = de_dz(k,j+1,i+1)
    986                 ENDIF
    987 
    988 !
    989 !--             If wall between gridpoint 2 and gridpoint 6, then
    990 !--             Neumann boundary condition has to be applied
    991                 IF ( gp_outside_of_building(2) == 1  .AND. &
    992                      gp_outside_of_building(6) == 0 )  THEN
    993                    num_gp = num_gp + 1
    994                    location(num_gp,1) = i * dx + 0.5_wp * dx
    995                    location(num_gp,2) = (j+1) * dy
    996                    location(num_gp,3) = k * dz - 0.5_wp * dz
    997                    ei(num_gp)     = e(k,j+1,i)
    998                    dissi(num_gp)  = diss(k,j+1,i)
    999                    de_dxi(num_gp) = 0.0_wp
    1000                    de_dyi(num_gp) = de_dy(k,j+1,i)
    1001                    de_dzi(num_gp) = de_dz(k,j+1,i)
    1002                 ENDIF
    1003 
    1004                 IF ( gp_outside_of_building(6) == 1  .AND. &
    1005                      gp_outside_of_building(2) == 0 )  THEN
    1006                    num_gp = num_gp + 1
    1007                    location(num_gp,1) = i * dx + 0.5_wp * dx
    1008                    location(num_gp,2) = (j+1) * dy
    1009                    location(num_gp,3) = k * dz - 0.5_wp * dz
    1010                    ei(num_gp)     = e(k,j+1,i+1)
    1011                    dissi(num_gp)  = diss(k,j+1,i+1)
    1012                    de_dxi(num_gp) = 0.0_wp
    1013779                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
    1014780                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
     
    1016782
    1017783!
    1018 !--             If wall between gridpoint 1 and gridpoint 2, then
    1019 !--             Neumann boundary condition has to be applied
    1020                 IF ( gp_outside_of_building(1) == 1  .AND. &
    1021                      gp_outside_of_building(2) == 0 )  THEN
     784!--             Determine vertical index of topography top at (j,i+1)
     785                k_wall = get_topography_top_index( j, i+1, 's' )
     786
     787                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    1022788                   num_gp = num_gp + 1
    1023                    location(num_gp,1) = i * dx
    1024                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1025                    location(num_gp,3) = k * dz - 0.5_wp * dz
    1026                    ei(num_gp)     = e(k,j,i)
    1027                    dissi(num_gp)  = diss(k,j,i)
    1028                    de_dxi(num_gp) = de_dx(k,j,i)
    1029                    de_dyi(num_gp) = 0.0_wp
    1030                    de_dzi(num_gp) = de_dz(k,j,i) 
    1031                 ENDIF
    1032 
    1033                 IF ( gp_outside_of_building(2) == 1  .AND. &
    1034                      gp_outside_of_building(1) == 0 )  THEN
    1035                    num_gp = num_gp + 1
    1036                    location(num_gp,1) = i * dx
    1037                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1038                    location(num_gp,3) = k * dz - 0.5_wp * dz
    1039                    ei(num_gp)     = e(k,j+1,i)
    1040                    dissi(num_gp)  = diss(k,j+1,i)
    1041                    de_dxi(num_gp) = de_dx(k,j+1,i)
    1042                    de_dyi(num_gp) = 0.0_wp
    1043                    de_dzi(num_gp) = de_dz(k,j+1,i)
    1044                 ENDIF
    1045 
    1046 !
    1047 !--             If wall between gridpoint 3 and gridpoint 7, then
    1048 !--             Neumann boundary condition has to be applied
    1049                 IF ( gp_outside_of_building(3) == 1  .AND. &
    1050                      gp_outside_of_building(7) == 0 )  THEN
    1051                    num_gp = num_gp + 1
    1052                    location(num_gp,1) = i * dx + 0.5_wp * dx
    1053                    location(num_gp,2) = j * dy
    1054                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1055                    ei(num_gp)     = e(k+1,j,i)
    1056                    dissi(num_gp)  = diss(k+1,j,i)
    1057                    de_dxi(num_gp) = 0.0_wp
    1058                    de_dyi(num_gp) = de_dy(k+1,j,i)
    1059                    de_dzi(num_gp) = de_dz(k+1,j,i)
    1060                 ENDIF
    1061 
    1062                 IF ( gp_outside_of_building(7) == 1  .AND. &
    1063                      gp_outside_of_building(3) == 0 )  THEN
    1064                    num_gp = num_gp + 1
    1065                    location(num_gp,1) = i * dx + 0.5_wp * dx
    1066                    location(num_gp,2) = j * dy
    1067                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1068                    ei(num_gp)     = e(k+1,j,i+1)
    1069                    dissi(num_gp)  = diss(k+1,j,i+1)
    1070                    de_dxi(num_gp) = 0.0_wp
    1071                    de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1072                    de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1073                 ENDIF
    1074 
    1075 !
    1076 !--             If wall between gridpoint 7 and gridpoint 8, then
    1077 !--             Neumann boundary condition has to be applied
    1078                 IF ( gp_outside_of_building(7) == 1  .AND. &
    1079                      gp_outside_of_building(8) == 0 )  THEN
    1080                    num_gp = num_gp + 1
    1081                    location(num_gp,1) = (i+1) * dx
    1082                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1083                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1084                    ei(num_gp)     = e(k+1,j,i+1)
    1085                    dissi(num_gp)  = diss(k+1,j,i+1)
    1086                    de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1087                    de_dyi(num_gp) = 0.0_wp
    1088                    de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1089                 ENDIF
    1090 
    1091                 IF ( gp_outside_of_building(8) == 1  .AND. &
    1092                      gp_outside_of_building(7) == 0 )  THEN
    1093                    num_gp = num_gp + 1
    1094                    location(num_gp,1) = (i+1) * dx
    1095                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1096                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1097                    ei(num_gp)     = e(k+1,j+1,i+1)
    1098                    dissi(num_gp)  = diss(k+1,j+1,i+1)
    1099                    de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1100                    de_dyi(num_gp) = 0.0_wp
    1101                    de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1102                 ENDIF
    1103 
    1104 !
    1105 !--             If wall between gridpoint 4 and gridpoint 8, then
    1106 !--             Neumann boundary condition has to be applied
    1107                 IF ( gp_outside_of_building(4) == 1  .AND. &
    1108                      gp_outside_of_building(8) == 0 )  THEN
    1109                    num_gp = num_gp + 1
    1110                    location(num_gp,1) = i * dx + 0.5_wp * dx
    1111                    location(num_gp,2) = (j+1) * dy
    1112                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1113                    ei(num_gp)     = e(k+1,j+1,i)
    1114                    dissi(num_gp)  = diss(k+1,j+1,i)
    1115                    de_dxi(num_gp) = 0.0_wp
    1116                    de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1117                    de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1118                 ENDIF
    1119 
    1120                 IF ( gp_outside_of_building(8) == 1  .AND. &
    1121                      gp_outside_of_building(4) == 0 )  THEN
    1122                    num_gp = num_gp + 1
    1123                    location(num_gp,1) = i * dx + 0.5_wp * dx
    1124                    location(num_gp,2) = (j+1) * dy
    1125                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1126                    ei(num_gp)     = e(k+1,j+1,i+1)
    1127                    dissi(num_gp)  = diss(k+1,j+1,i+1)
    1128                    de_dxi(num_gp) = 0.0_wp
    1129                    de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1130                    de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1131                 ENDIF
    1132 
    1133 !
    1134 !--             If wall between gridpoint 3 and gridpoint 4, then
    1135 !--             Neumann boundary condition has to be applied
    1136                 IF ( gp_outside_of_building(3) == 1  .AND. &
    1137                      gp_outside_of_building(4) == 0 )  THEN
    1138                    num_gp = num_gp + 1
    1139                    location(num_gp,1) = i * dx
    1140                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1141                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1142                    ei(num_gp)     = e(k+1,j,i)
    1143                    dissi(num_gp)  = diss(k+1,j,i)
    1144                    de_dxi(num_gp) = de_dx(k+1,j,i)
    1145                    de_dyi(num_gp) = 0.0_wp
    1146                    de_dzi(num_gp) = de_dz(k+1,j,i)
    1147                 ENDIF
    1148 
    1149                 IF ( gp_outside_of_building(4) == 1  .AND. &
    1150                      gp_outside_of_building(3) == 0 )  THEN
    1151                    num_gp = num_gp + 1
    1152                    location(num_gp,1) = i * dx
    1153                    location(num_gp,2) = j * dy + 0.5_wp * dy
    1154                    location(num_gp,3) = k * dz + 0.5_wp * dz
    1155                    ei(num_gp)     = e(k+1,j+1,i)
    1156                    dissi(num_gp)  = diss(k+1,j+1,i)
    1157                    de_dxi(num_gp) = de_dx(k+1,j+1,i)
    1158                    de_dyi(num_gp) = 0.0_wp
    1159                    de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1160                 ENDIF
    1161 
    1162 !
    1163 !--             If wall between gridpoint 1 and gridpoint 3, then
    1164 !--             Neumann boundary condition has to be applied
    1165 !--             (only one case as only building beneath is possible)
    1166                 IF ( gp_outside_of_building(1) == 0  .AND. &
    1167                      gp_outside_of_building(3) == 1 )  THEN
    1168                    num_gp = num_gp + 1
    1169                    location(num_gp,1) = i * dx
    1170                    location(num_gp,2) = j * dy
    1171                    location(num_gp,3) = k * dz
    1172                    ei(num_gp)     = e(k+1,j,i)
    1173                    dissi(num_gp)  = diss(k+1,j,i)
    1174                    de_dxi(num_gp) = de_dx(k+1,j,i)
    1175                    de_dyi(num_gp) = de_dy(k+1,j,i)
    1176                    de_dzi(num_gp) = 0.0_wp
    1177                 ENDIF
    1178 
    1179 !
    1180 !--             If wall between gridpoint 5 and gridpoint 7, then
    1181 !--             Neumann boundary condition has to be applied
    1182 !--             (only one case as only building beneath is possible)
    1183                 IF ( gp_outside_of_building(5) == 0  .AND. &
    1184                      gp_outside_of_building(7) == 1 )  THEN
    1185                    num_gp = num_gp + 1
     789                   gp_outside_of_building(7) = 1
    1186790                   location(num_gp,1) = (i+1) * dx
    1187791                   location(num_gp,2) = j * dy
    1188                    location(num_gp,3) = k * dz
     792                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    1189793                   ei(num_gp)     = e(k+1,j,i+1)
    1190794                   dissi(num_gp)  = diss(k+1,j,i+1)
    1191795                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1192796                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1193                    de_dzi(num_gp) = 0.0_wp
    1194                 ENDIF
    1195 
    1196 !
    1197 !--             If wall between gridpoint 2 and gridpoint 4, then
    1198 !--             Neumann boundary condition has to be applied
    1199 !--             (only one case as only building beneath is possible)
    1200                 IF ( gp_outside_of_building(2) == 0  .AND. &
    1201                      gp_outside_of_building(4) == 1 )  THEN
     797                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
     798                ENDIF
     799
     800!
     801!--             Determine vertical index of topography top at (j+1,i+1)
     802                k_wall = get_topography_top_index( j+1, i+1, 's' )
     803
     804                IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
    1202805                   num_gp = num_gp + 1
    1203                    location(num_gp,1) = i * dx
    1204                    location(num_gp,2) = (j+1) * dy
    1205                    location(num_gp,3) = k * dz
    1206                    ei(num_gp)     = e(k+1,j+1,i)
    1207                    dissi(num_gp)  = diss(k+1,j+1,i)
    1208                    de_dxi(num_gp) = de_dx(k+1,j+1,i)
    1209                    de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1210                    de_dzi(num_gp) = 0.0_wp
    1211                 ENDIF
    1212 
    1213 !
    1214 !--             If wall between gridpoint 6 and gridpoint 8, then
    1215 !--             Neumann boundary condition has to be applied
    1216 !--             (only one case as only building beneath is possible)
    1217                 IF ( gp_outside_of_building(6) == 0  .AND. &
    1218                      gp_outside_of_building(8) == 1 )  THEN
    1219                    num_gp = num_gp + 1
     806                   gp_outside_of_building(8) = 1
    1220807                   location(num_gp,1) = (i+1) * dx
    1221808                   location(num_gp,2) = (j+1) * dy
    1222                    location(num_gp,3) = k * dz
     809                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    1223810                   ei(num_gp)     = e(k+1,j+1,i+1)
    1224811                   dissi(num_gp)  = diss(k+1,j+1,i+1)
    1225812                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1226813                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1227                    de_dzi(num_gp) = 0.0_wp
    1228                 ENDIF
    1229  
    1230 !
    1231 !--             Carry out the interpolation
    1232                 IF ( num_gp == 1 )  THEN
     814                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     815                ENDIF
     816!
     817!--             If all gridpoints are situated outside of a building, then the
     818!--             ordinary interpolation scheme can be used.
     819                IF ( num_gp == 8 )  THEN
     820
     821                   x  = particles(n)%x - i * dx
     822                   y  = particles(n)%y - j * dy
     823                   aa = x**2          + y**2
     824                   bb = ( dx - x )**2 + y**2
     825                   cc = x**2          + ( dy - y )**2
     826                   dd = ( dx - x )**2 + ( dy - y )**2
     827                   gg = aa + bb + cc + dd
     828     
     829                   e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
     830                             + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
     831                             ) / ( 3.0_wp * gg )
     832     
     833                   IF ( k == nzt )  THEN
     834                      e_int(n) = e_int_l
     835                   ELSE
     836                      e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
     837                                  ( gg - bb ) * e(k+1,j,i+1) + &
     838                                  ( gg - cc ) * e(k+1,j+1,i) + &
     839                                  ( gg - dd ) * e(k+1,j+1,i+1) &
     840                                ) / ( 3.0_wp * gg )
     841                      e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
     842                                          ( e_int_u - e_int_l )
     843                   ENDIF
    1233844!
    1234 !--                If only one of the gridpoints is situated outside of the
    1235 !--                building, it follows that the values at the particle
    1236 !--                location are the same as the gridpoint values
    1237                    e_int(n)     = ei(num_gp)   
    1238                    diss_int(n)  = dissi(num_gp)
    1239                    de_dx_int(n) = de_dxi(num_gp)
    1240                    de_dy_int(n) = de_dyi(num_gp)
    1241                    de_dz_int(n) = de_dzi(num_gp)
     845!--                Needed to avoid NaN particle velocities (this might not be
     846!--                required any more)
     847                   IF ( e_int(n) <= 0.0_wp )  THEN
     848                      e_int(n) = 1.0E-20_wp
     849                   ENDIF
     850!
     851!--                Interpolate the TKE gradient along x (adopt incides i,j,k
     852!--                and all position variables from above (TKE))
     853                   de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
     854                                   ( gg - bb ) * de_dx(k,j,i+1) + &
     855                                   ( gg - cc ) * de_dx(k,j+1,i) + &
     856                                   ( gg - dd ) * de_dx(k,j+1,i+1) &
     857                                 ) / ( 3.0_wp * gg )
     858
     859                   IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
     860                      de_dx_int(n) = de_dx_int_l
     861                   ELSE
     862                      de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
     863                                      ( gg - bb ) * de_dx(k+1,j,i+1) + &
     864                                      ( gg - cc ) * de_dx(k+1,j+1,i) + &
     865                                      ( gg - dd ) * de_dx(k+1,j+1,i+1) &
     866                                    ) / ( 3.0_wp * gg )
     867                      de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
     868                                              dz * ( de_dx_int_u - de_dx_int_l )
     869                   ENDIF
     870
     871!
     872!--                Interpolate the TKE gradient along y
     873                   de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
     874                                   ( gg - bb ) * de_dy(k,j,i+1) + &
     875                                   ( gg - cc ) * de_dy(k,j+1,i) + &
     876                                   ( gg - dd ) * de_dy(k,j+1,i+1) &
     877                                 ) / ( 3.0_wp * gg )
     878                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     879                      de_dy_int(n) = de_dy_int_l
     880                   ELSE
     881                      de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
     882                                      ( gg - bb ) * de_dy(k+1,j,i+1) + &
     883                                      ( gg - cc ) * de_dy(k+1,j+1,i) + &
     884                                      ( gg - dd ) * de_dy(k+1,j+1,i+1) &
     885                                    ) / ( 3.0_wp * gg )
     886                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
     887                                              dz * ( de_dy_int_u - de_dy_int_l )
     888                   ENDIF
     889
     890!
     891!--                Interpolate the TKE gradient along z
     892                   IF ( zv(n) < 0.5_wp * dz )  THEN
     893                      de_dz_int(n) = 0.0_wp
     894                   ELSE
     895                      de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
     896                                      ( gg - bb ) * de_dz(k,j,i+1) + &
     897                                      ( gg - cc ) * de_dz(k,j+1,i) + &
     898                                      ( gg - dd ) * de_dz(k,j+1,i+1) &
     899                                    ) / ( 3.0_wp * gg )
     900   
     901                      IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     902                         de_dz_int(n) = de_dz_int_l
     903                      ELSE
     904                         de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
     905                                         ( gg - bb ) * de_dz(k+1,j,i+1) + &
     906                                         ( gg - cc ) * de_dz(k+1,j+1,i) + &
     907                                         ( gg - dd ) * de_dz(k+1,j+1,i+1) &
     908                                       ) / ( 3.0_wp * gg )
     909                         de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
     910                                              dz * ( de_dz_int_u - de_dz_int_l )
     911                      ENDIF
     912                   ENDIF
     913
     914!
     915!--                Interpolate the dissipation of TKE
     916                   diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
     917                                  ( gg - bb ) * diss(k,j,i+1) + &
     918                                  ( gg - cc ) * diss(k,j+1,i) + &
     919                                  ( gg - dd ) * diss(k,j+1,i+1) &
     920                                ) / ( 3.0_wp * gg )
     921   
     922                   IF ( k == nzt )  THEN
     923                      diss_int(n) = diss_int_l
     924                   ELSE
     925                      diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
     926                                     ( gg - bb ) * diss(k+1,j,i+1) + &
     927                                     ( gg - cc ) * diss(k+1,j+1,i) + &
     928                                     ( gg - dd ) * diss(k+1,j+1,i+1) &
     929                                   ) / ( 3.0_wp * gg )
     930                      diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
     931                                              ( diss_int_u - diss_int_l )
     932                   ENDIF
     933!
     934!--                Set flag for stochastic equation.
     935                   term_1_2(n) = 1.0_wp
     936   
     937                ELSE
     938   
     939!
     940!--                If wall between gridpoint 1 and gridpoint 5, then
     941!--                Neumann boundary condition has to be applied
     942                   IF ( gp_outside_of_building(1) == 1  .AND. &
     943                        gp_outside_of_building(5) == 0 )  THEN
     944                      num_gp = num_gp + 1
     945                      location(num_gp,1) = i * dx + 0.5_wp * dx
     946                      location(num_gp,2) = j * dy
     947                      location(num_gp,3) = k * dz - 0.5_wp * dz
     948                      ei(num_gp)     = e(k,j,i)
     949                      dissi(num_gp)  = diss(k,j,i)
     950                      de_dxi(num_gp) = 0.0_wp
     951                      de_dyi(num_gp) = de_dy(k,j,i)
     952                      de_dzi(num_gp) = de_dz(k,j,i)
     953                   ENDIF
     954
     955                   IF ( gp_outside_of_building(5) == 1  .AND. &
     956                        gp_outside_of_building(1) == 0 )  THEN
     957                      num_gp = num_gp + 1
     958                      location(num_gp,1) = i * dx + 0.5_wp * dx
     959                      location(num_gp,2) = j * dy
     960                      location(num_gp,3) = k * dz - 0.5_wp * dz
     961                      ei(num_gp)     = e(k,j,i+1)
     962                      dissi(num_gp)  = diss(k,j,i+1)
     963                      de_dxi(num_gp) = 0.0_wp
     964                      de_dyi(num_gp) = de_dy(k,j,i+1)
     965                      de_dzi(num_gp) = de_dz(k,j,i+1)
     966                   ENDIF
     967
     968!
     969!--                If wall between gridpoint 5 and gridpoint 6, then
     970!--                then Neumann boundary condition has to be applied
     971                   IF ( gp_outside_of_building(5) == 1  .AND. &
     972                        gp_outside_of_building(6) == 0 )  THEN
     973                      num_gp = num_gp + 1
     974                      location(num_gp,1) = (i+1) * dx
     975                      location(num_gp,2) = j * dy + 0.5_wp * dy
     976                      location(num_gp,3) = k * dz - 0.5_wp * dz
     977                      ei(num_gp)     = e(k,j,i+1)
     978                      dissi(num_gp)  = diss(k,j,i+1)
     979                      de_dxi(num_gp) = de_dx(k,j,i+1)
     980                      de_dyi(num_gp) = 0.0_wp
     981                      de_dzi(num_gp) = de_dz(k,j,i+1)
     982                   ENDIF
     983
     984                   IF ( gp_outside_of_building(6) == 1  .AND. &
     985                        gp_outside_of_building(5) == 0 )  THEN
     986                      num_gp = num_gp + 1
     987                      location(num_gp,1) = (i+1) * dx
     988                      location(num_gp,2) = j * dy + 0.5_wp * dy
     989                      location(num_gp,3) = k * dz - 0.5_wp * dz
     990                      ei(num_gp)     = e(k,j+1,i+1)
     991                      dissi(num_gp)  = diss(k,j+1,i+1)
     992                      de_dxi(num_gp) = de_dx(k,j+1,i+1)
     993                      de_dyi(num_gp) = 0.0_wp
     994                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
     995                   ENDIF
     996
     997!
     998!--                If wall between gridpoint 2 and gridpoint 6, then
     999!--                Neumann boundary condition has to be applied
     1000                   IF ( gp_outside_of_building(2) == 1  .AND. &
     1001                        gp_outside_of_building(6) == 0 )  THEN
     1002                      num_gp = num_gp + 1
     1003                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1004                      location(num_gp,2) = (j+1) * dy
     1005                      location(num_gp,3) = k * dz - 0.5_wp * dz
     1006                      ei(num_gp)     = e(k,j+1,i)
     1007                      dissi(num_gp)  = diss(k,j+1,i)
     1008                      de_dxi(num_gp) = 0.0_wp
     1009                      de_dyi(num_gp) = de_dy(k,j+1,i)
     1010                      de_dzi(num_gp) = de_dz(k,j+1,i)
     1011                   ENDIF
     1012
     1013                   IF ( gp_outside_of_building(6) == 1  .AND. &
     1014                        gp_outside_of_building(2) == 0 )  THEN
     1015                      num_gp = num_gp + 1
     1016                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1017                      location(num_gp,2) = (j+1) * dy
     1018                      location(num_gp,3) = k * dz - 0.5_wp * dz
     1019                      ei(num_gp)     = e(k,j+1,i+1)
     1020                      dissi(num_gp)  = diss(k,j+1,i+1)
     1021                      de_dxi(num_gp) = 0.0_wp
     1022                      de_dyi(num_gp) = de_dy(k,j+1,i+1)
     1023                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
     1024                   ENDIF
     1025
     1026!
     1027!--                If wall between gridpoint 1 and gridpoint 2, then
     1028!--                Neumann boundary condition has to be applied
     1029                   IF ( gp_outside_of_building(1) == 1  .AND. &
     1030                        gp_outside_of_building(2) == 0 )  THEN
     1031                      num_gp = num_gp + 1
     1032                      location(num_gp,1) = i * dx
     1033                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1034                      location(num_gp,3) = k * dz - 0.5_wp * dz
     1035                      ei(num_gp)     = e(k,j,i)
     1036                      dissi(num_gp)  = diss(k,j,i)
     1037                      de_dxi(num_gp) = de_dx(k,j,i)
     1038                      de_dyi(num_gp) = 0.0_wp
     1039                      de_dzi(num_gp) = de_dz(k,j,i) 
     1040                   ENDIF
     1041
     1042                   IF ( gp_outside_of_building(2) == 1  .AND. &
     1043                        gp_outside_of_building(1) == 0 )  THEN
     1044                      num_gp = num_gp + 1
     1045                      location(num_gp,1) = i * dx
     1046                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1047                      location(num_gp,3) = k * dz - 0.5_wp * dz
     1048                      ei(num_gp)     = e(k,j+1,i)
     1049                      dissi(num_gp)  = diss(k,j+1,i)
     1050                      de_dxi(num_gp) = de_dx(k,j+1,i)
     1051                      de_dyi(num_gp) = 0.0_wp
     1052                      de_dzi(num_gp) = de_dz(k,j+1,i)
     1053                   ENDIF
     1054
     1055!
     1056!--                If wall between gridpoint 3 and gridpoint 7, then
     1057!--                Neumann boundary condition has to be applied
     1058                   IF ( gp_outside_of_building(3) == 1  .AND. &
     1059                        gp_outside_of_building(7) == 0 )  THEN
     1060                      num_gp = num_gp + 1
     1061                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1062                      location(num_gp,2) = j * dy
     1063                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1064                      ei(num_gp)     = e(k+1,j,i)
     1065                      dissi(num_gp)  = diss(k+1,j,i)
     1066                      de_dxi(num_gp) = 0.0_wp
     1067                      de_dyi(num_gp) = de_dy(k+1,j,i)
     1068                      de_dzi(num_gp) = de_dz(k+1,j,i)
     1069                   ENDIF
     1070
     1071                   IF ( gp_outside_of_building(7) == 1  .AND. &
     1072                        gp_outside_of_building(3) == 0 )  THEN
     1073                      num_gp = num_gp + 1
     1074                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1075                      location(num_gp,2) = j * dy
     1076                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1077                      ei(num_gp)     = e(k+1,j,i+1)
     1078                      dissi(num_gp)  = diss(k+1,j,i+1)
     1079                      de_dxi(num_gp) = 0.0_wp
     1080                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
     1081                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
     1082                   ENDIF
     1083
     1084!
     1085!--                If wall between gridpoint 7 and gridpoint 8, then
     1086!--                Neumann boundary condition has to be applied
     1087                   IF ( gp_outside_of_building(7) == 1  .AND. &
     1088                        gp_outside_of_building(8) == 0 )  THEN
     1089                      num_gp = num_gp + 1
     1090                      location(num_gp,1) = (i+1) * dx
     1091                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1092                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1093                      ei(num_gp)     = e(k+1,j,i+1)
     1094                      dissi(num_gp)  = diss(k+1,j,i+1)
     1095                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
     1096                      de_dyi(num_gp) = 0.0_wp
     1097                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
     1098                   ENDIF
     1099
     1100                   IF ( gp_outside_of_building(8) == 1  .AND. &
     1101                        gp_outside_of_building(7) == 0 )  THEN
     1102                      num_gp = num_gp + 1
     1103                      location(num_gp,1) = (i+1) * dx
     1104                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1105                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1106                      ei(num_gp)     = e(k+1,j+1,i+1)
     1107                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1108                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
     1109                      de_dyi(num_gp) = 0.0_wp
     1110                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     1111                   ENDIF
     1112
     1113!
     1114!--                If wall between gridpoint 4 and gridpoint 8, then
     1115!--                Neumann boundary condition has to be applied
     1116                   IF ( gp_outside_of_building(4) == 1  .AND. &
     1117                        gp_outside_of_building(8) == 0 )  THEN
     1118                      num_gp = num_gp + 1
     1119                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1120                      location(num_gp,2) = (j+1) * dy
     1121                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1122                      ei(num_gp)     = e(k+1,j+1,i)
     1123                      dissi(num_gp)  = diss(k+1,j+1,i)
     1124                      de_dxi(num_gp) = 0.0_wp
     1125                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
     1126                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
     1127                   ENDIF
     1128
     1129                   IF ( gp_outside_of_building(8) == 1  .AND. &
     1130                        gp_outside_of_building(4) == 0 )  THEN
     1131                      num_gp = num_gp + 1
     1132                      location(num_gp,1) = i * dx + 0.5_wp * dx
     1133                      location(num_gp,2) = (j+1) * dy
     1134                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1135                      ei(num_gp)     = e(k+1,j+1,i+1)
     1136                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1137                      de_dxi(num_gp) = 0.0_wp
     1138                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     1139                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     1140                   ENDIF
     1141
     1142!
     1143!--                If wall between gridpoint 3 and gridpoint 4, then
     1144!--                Neumann boundary condition has to be applied
     1145                   IF ( gp_outside_of_building(3) == 1  .AND. &
     1146                        gp_outside_of_building(4) == 0 )  THEN
     1147                      num_gp = num_gp + 1
     1148                      location(num_gp,1) = i * dx
     1149                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1150                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1151                      ei(num_gp)     = e(k+1,j,i)
     1152                      dissi(num_gp)  = diss(k+1,j,i)
     1153                      de_dxi(num_gp) = de_dx(k+1,j,i)
     1154                      de_dyi(num_gp) = 0.0_wp
     1155                      de_dzi(num_gp) = de_dz(k+1,j,i)
     1156                   ENDIF
     1157
     1158                   IF ( gp_outside_of_building(4) == 1  .AND. &
     1159                        gp_outside_of_building(3) == 0 )  THEN
     1160                      num_gp = num_gp + 1
     1161                      location(num_gp,1) = i * dx
     1162                      location(num_gp,2) = j * dy + 0.5_wp * dy
     1163                      location(num_gp,3) = k * dz + 0.5_wp * dz
     1164                      ei(num_gp)     = e(k+1,j+1,i)
     1165                      dissi(num_gp)  = diss(k+1,j+1,i)
     1166                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
     1167                      de_dyi(num_gp) = 0.0_wp
     1168                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
     1169                   ENDIF
     1170
     1171!
     1172!--                If wall between gridpoint 1 and gridpoint 3, then
     1173!--                Neumann boundary condition has to be applied
     1174!--                (only one case as only building beneath is possible)
     1175                   IF ( gp_outside_of_building(1) == 0  .AND. &
     1176                        gp_outside_of_building(3) == 1 )  THEN
     1177                      num_gp = num_gp + 1
     1178                      location(num_gp,1) = i * dx
     1179                      location(num_gp,2) = j * dy
     1180                      location(num_gp,3) = k * dz
     1181                      ei(num_gp)     = e(k+1,j,i)
     1182                      dissi(num_gp)  = diss(k+1,j,i)
     1183                      de_dxi(num_gp) = de_dx(k+1,j,i)
     1184                      de_dyi(num_gp) = de_dy(k+1,j,i)
     1185                      de_dzi(num_gp) = 0.0_wp
     1186                   ENDIF
     1187
     1188!
     1189!--                If wall between gridpoint 5 and gridpoint 7, then
     1190!--                Neumann boundary condition has to be applied
     1191!--                (only one case as only building beneath is possible)
     1192                   IF ( gp_outside_of_building(5) == 0  .AND. &
     1193                        gp_outside_of_building(7) == 1 )  THEN
     1194                      num_gp = num_gp + 1
     1195                      location(num_gp,1) = (i+1) * dx
     1196                      location(num_gp,2) = j * dy
     1197                      location(num_gp,3) = k * dz
     1198                      ei(num_gp)     = e(k+1,j,i+1)
     1199                      dissi(num_gp)  = diss(k+1,j,i+1)
     1200                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
     1201                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
     1202                      de_dzi(num_gp) = 0.0_wp
     1203                   ENDIF
     1204
     1205!
     1206!--                If wall between gridpoint 2 and gridpoint 4, then
     1207!--                Neumann boundary condition has to be applied
     1208!--                (only one case as only building beneath is possible)
     1209                   IF ( gp_outside_of_building(2) == 0  .AND. &
     1210                        gp_outside_of_building(4) == 1 )  THEN
     1211                      num_gp = num_gp + 1
     1212                      location(num_gp,1) = i * dx
     1213                      location(num_gp,2) = (j+1) * dy
     1214                      location(num_gp,3) = k * dz
     1215                      ei(num_gp)     = e(k+1,j+1,i)
     1216                      dissi(num_gp)  = diss(k+1,j+1,i)
     1217                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
     1218                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
     1219                      de_dzi(num_gp) = 0.0_wp
     1220                   ENDIF
     1221
     1222!
     1223!--                If wall between gridpoint 6 and gridpoint 8, then
     1224!--                Neumann boundary condition has to be applied
     1225!--                (only one case as only building beneath is possible)
     1226                   IF ( gp_outside_of_building(6) == 0  .AND. &
     1227                        gp_outside_of_building(8) == 1 )  THEN
     1228                      num_gp = num_gp + 1
     1229                      location(num_gp,1) = (i+1) * dx
     1230                      location(num_gp,2) = (j+1) * dy
     1231                      location(num_gp,3) = k * dz
     1232                      ei(num_gp)     = e(k+1,j+1,i+1)
     1233                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1234                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
     1235                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     1236                      de_dzi(num_gp) = 0.0_wp
     1237                   ENDIF
     1238     
     1239!
     1240!--                Carry out the interpolation
     1241                   IF ( num_gp == 1 )  THEN
     1242!
     1243!--                   If only one of the gridpoints is situated outside of the
     1244!--                   building, it follows that the values at the particle
     1245!--                   location are the same as the gridpoint values
     1246                      e_int(n)     = ei(num_gp)   
     1247                      diss_int(n)  = dissi(num_gp)
     1248                      de_dx_int(n) = de_dxi(num_gp)
     1249                      de_dy_int(n) = de_dyi(num_gp)
     1250                      de_dz_int(n) = de_dzi(num_gp)
     1251!
     1252!--                   Set flag for stochastic equation which disables calculation
     1253!--                   of drift and memory term near topography.
     1254                      term_1_2(n) = 0.0_wp
     1255                   ELSE IF ( num_gp > 1 )  THEN
     1256   
     1257                      d_sum = 0.0_wp
     1258!
     1259!--                   Evaluation of the distances between the gridpoints
     1260!--                   contributing to the interpolated values, and the particle
     1261!--                   location
     1262                      DO  agp = 1, num_gp
     1263                         d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
     1264                                      + ( particles(n)%y-location(agp,2) )**2  &
     1265                                      + ( zv(n)-location(agp,3) )**2
     1266                         d_sum        = d_sum + d_gp_pl(agp)
     1267                      ENDDO
     1268   
     1269!
     1270!--                  Finally the interpolation can be carried out
     1271                      e_int(n)     = 0.0_wp
     1272                      diss_int(n)  = 0.0_wp
     1273                      de_dx_int(n) = 0.0_wp 
     1274                      de_dy_int(n) = 0.0_wp 
     1275                      de_dz_int(n) = 0.0_wp 
     1276                      DO  agp = 1, num_gp
     1277                         e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
     1278                                                ei(agp) / ( (num_gp-1) * d_sum )
     1279                         diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
     1280                                             dissi(agp) / ( (num_gp-1) * d_sum )
     1281                         de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
     1282                                            de_dxi(agp) / ( (num_gp-1) * d_sum )
     1283                         de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
     1284                                            de_dyi(agp) / ( (num_gp-1) * d_sum )
     1285                         de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
     1286                                            de_dzi(agp) / ( (num_gp-1) * d_sum )
     1287                      ENDDO
     1288   
     1289                   ENDIF
     1290                   e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
     1291                   diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
     1292                   de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
     1293                   de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) )
     1294                   de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) )
    12421295!
    12431296!--                Set flag for stochastic equation which disables calculation
    12441297!--                of drift and memory term near topography.
    12451298                   term_1_2(n) = 0.0_wp
    1246                 ELSE IF ( num_gp > 1 )  THEN
    1247  
    1248                    d_sum = 0.0_wp
    1249 !
    1250 !--                Evaluation of the distances between the gridpoints
    1251 !--                contributing to the interpolated values, and the particle
    1252 !--                location
    1253                    DO  agp = 1, num_gp
    1254                       d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
    1255                                    + ( particles(n)%y-location(agp,2) )**2  &
    1256                                    + ( zv(n)-location(agp,3) )**2
    1257                       d_sum        = d_sum + d_gp_pl(agp)
    1258                    ENDDO
    1259  
    1260 !
    1261 !--                Finally the interpolation can be carried out
    1262                    e_int(n)     = 0.0_wp
    1263                    diss_int(n)  = 0.0_wp
    1264                    de_dx_int(n) = 0.0_wp 
    1265                    de_dy_int(n) = 0.0_wp 
    1266                    de_dz_int(n) = 0.0_wp 
    1267                    DO  agp = 1, num_gp
    1268                       e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
    1269                                              ei(agp) / ( (num_gp-1) * d_sum )
    1270                       diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
    1271                                           dissi(agp) / ( (num_gp-1) * d_sum )
    1272                       de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1273                                          de_dxi(agp) / ( (num_gp-1) * d_sum )
    1274                       de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1275                                          de_dyi(agp) / ( (num_gp-1) * d_sum )
    1276                       de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1277                                          de_dzi(agp) / ( (num_gp-1) * d_sum )
    1278                    ENDDO
    1279  
    1280                 ENDIF
    1281                 e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
    1282                 diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
    1283                 de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
    1284                 de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) )
    1285                 de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) )
    1286 !
    1287 !--             Set flag for stochastic equation which disables calculation
    1288 !--             of drift and memory term near topography.
    1289                 term_1_2(n) = 0.0_wp
    1290              ENDIF
     1299                ENDIF
     1300             ENDDO
    12911301          ENDDO
    12921302       ENDIF
     
    13451355       ENDDO
    13461356
    1347        DO  n = 1, number_of_particles
    1348 
    1349          rg(n,1) = random_gauss( iran_part, 5.0_wp )
    1350          rg(n,2) = random_gauss( iran_part, 5.0_wp )
    1351          rg(n,3) = random_gauss( iran_part, 5.0_wp )
    1352 
     1357       DO  nb = 0, 7
     1358          DO  n = start_index(nb), end_index(nb)
     1359             rg(n,1) = random_gauss( iran_part, 5.0_wp )
     1360             rg(n,2) = random_gauss( iran_part, 5.0_wp )
     1361             rg(n,3) = random_gauss( iran_part, 5.0_wp )
     1362          ENDDO
    13531363       ENDDO
    13541364
    1355        DO  n = 1, number_of_particles
    1356 !
    1357 !--       Calculate the Lagrangian timescale according to Weil et al. (2004).
    1358           lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
    1359                            ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
    1360 
    1361 !
    1362 !--       Calculate the next particle timestep. dt_gap is the time needed to
    1363 !--       complete the current LES timestep.
    1364           dt_gap = dt_3d - particles(n)%dt_sum
    1365           dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
    1366 
    1367 !
    1368 !--       The particle timestep should not be too small in order to prevent
    1369 !--       the number of particle timesteps of getting too large
    1370           IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
    1371              dt_particle(n) = dt_min_part
    1372           ENDIF
    1373 
    1374 !
    1375 !--       Calculate the SGS velocity components
    1376           IF ( particles(n)%age == 0.0_wp )  THEN
    1377 !
    1378 !--          For new particles the SGS components are derived from the SGS
    1379 !--          TKE. Limit the Gaussian random number to the interval
    1380 !--          [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
    1381 !--          from becoming unrealistically large.
    1382              particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    1383                                   ( rg(n,1) - 1.0_wp )
    1384              particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    1385                                   ( rg(n,2) - 1.0_wp )
    1386              particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) *   &
    1387                                   ( rg(n,3) - 1.0_wp )
    1388 
    1389           ELSE
    1390 !
    1391 !--          Restriction of the size of the new timestep: compared to the
    1392 !--          previous timestep the increase must not exceed 200%
    1393 
    1394              dt_particle_m = particles(n)%age - particles(n)%age_m
    1395              IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
    1396                 dt_particle(n) = 2.0_wp * dt_particle_m
     1365       DO  nb = 0, 7
     1366          DO  n = start_index(nb), end_index(nb)
     1367
     1368!
     1369!--          Calculate the Lagrangian timescale according to Weil et al. (2004).
     1370             lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
     1371                              ( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
     1372
     1373!
     1374!--          Calculate the next particle timestep. dt_gap is the time needed to
     1375!--          complete the current LES timestep.
     1376             dt_gap = dt_3d - particles(n)%dt_sum
     1377             dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
     1378
     1379!
     1380!--          The particle timestep should not be too small in order to prevent
     1381!--          the number of particle timesteps of getting too large
     1382             IF ( dt_particle(n) < dt_min_part  .AND.  dt_min_part < dt_gap )  THEN
     1383                dt_particle(n) = dt_min_part
    13971384             ENDIF
    13981385
    13991386!
    1400 !--          For old particles the SGS components are correlated with the
    1401 !--          values from the previous timestep. Random numbers have also to
    1402 !--          be limited (see above).
    1403 !--          As negative values for the subgrid TKE are not allowed, the
    1404 !--          change of the subgrid TKE with time cannot be smaller than
    1405 !--          -e_int(n)/dt_particle. This value is used as a lower boundary
    1406 !--          value for the change of TKE
    1407 
    1408              de_dt_min = - e_int(n) / dt_particle(n)
    1409 
    1410              de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
    1411 
    1412              IF ( de_dt < de_dt_min )  THEN
    1413                 de_dt = de_dt_min
     1387!--          Calculate the SGS velocity components
     1388             IF ( particles(n)%age == 0.0_wp )  THEN
     1389!
     1390!--             For new particles the SGS components are derived from the SGS
     1391!--             TKE. Limit the Gaussian random number to the interval
     1392!--             [-5.0*sigma, 5.0*sigma] in order to prevent the SGS velocities
     1393!--             from becoming unrealistically large.
     1394                particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
     1395                                          + 1E-20_wp ) *                       &
     1396                                     ( rg(n,1) - 1.0_wp )
     1397                particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
     1398                                          + 1E-20_wp ) *                       &
     1399                                     ( rg(n,2) - 1.0_wp )
     1400                particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n)     &
     1401                                          + 1E-20_wp ) *                       &
     1402                                     ( rg(n,3) - 1.0_wp )
     1403
     1404             ELSE
     1405!
     1406!--             Restriction of the size of the new timestep: compared to the
     1407!--             previous timestep the increase must not exceed 200%. First,
     1408!--             check if age > age_m, in order to prevent that particles get zero
     1409!--             timestep.
     1410                dt_particle_m = MERGE( dt_particle(n),                         &
     1411                                       particles(n)%age - particles(n)%age_m,  &
     1412                                       particles(n)%age - particles(n)%age_m < &
     1413                                       1E-8_wp )
     1414                IF ( dt_particle(n) > 2.0_wp * dt_particle_m )  THEN
     1415                   dt_particle(n) = 2.0_wp * dt_particle_m
     1416                ENDIF
     1417
     1418!
     1419!--             For old particles the SGS components are correlated with the
     1420!--             values from the previous timestep. Random numbers have also to
     1421!--             be limited (see above).
     1422!--             As negative values for the subgrid TKE are not allowed, the
     1423!--             change of the subgrid TKE with time cannot be smaller than
     1424!--             -e_int(n)/dt_particle. This value is used as a lower boundary
     1425!--             value for the change of TKE
     1426                de_dt_min = - e_int(n) / dt_particle(n)
     1427
     1428                de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
     1429
     1430                IF ( de_dt < de_dt_min )  THEN
     1431                   de_dt = de_dt_min
     1432                ENDIF
     1433
     1434                CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),&
     1435                                        de_dx_int(n), de_dt, diss_int(n),       &
     1436                                        dt_particle(n), rg(n,1), term_1_2(n) )
     1437
     1438                CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),&
     1439                                        de_dy_int(n), de_dt, diss_int(n),       &
     1440                                        dt_particle(n), rg(n,2), term_1_2(n) )
     1441
     1442                CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),&
     1443                                        de_dz_int(n), de_dt, diss_int(n),       &
     1444                                        dt_particle(n), rg(n,3), term_1_2(n) )
     1445
    14141446             ENDIF
    14151447
    1416              CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n),  &
    1417                                      de_dx_int(n), de_dt, diss_int(n),         &
    1418                                      dt_particle(n), rg(n,1), term_1_2(n) )
    1419 
    1420              CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n),  &
    1421                                      de_dy_int(n), de_dt, diss_int(n),         &
    1422                                      dt_particle(n), rg(n,2), term_1_2(n) )
    1423 
    1424              CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n),  &
    1425                                      de_dz_int(n), de_dt, diss_int(n),         &
    1426                                      dt_particle(n), rg(n,3), term_1_2(n) )
    1427 
    1428           ENDIF
    1429 
    1430           u_int(n) = u_int(n) + particles(n)%rvar1
    1431           v_int(n) = v_int(n) + particles(n)%rvar2
    1432           w_int(n) = w_int(n) + particles(n)%rvar3
    1433 !
    1434 !--       Store the SGS TKE of the current timelevel which is needed for
    1435 !--       for calculating the SGS particle velocities at the next timestep
    1436           particles(n)%e_m = e_int(n)
     1448             u_int(n) = u_int(n) + particles(n)%rvar1
     1449             v_int(n) = v_int(n) + particles(n)%rvar2
     1450             w_int(n) = w_int(n) + particles(n)%rvar3
     1451!
     1452!--          Store the SGS TKE of the current timelevel which is needed for
     1453!--          for calculating the SGS particle velocities at the next timestep
     1454             particles(n)%e_m = e_int(n)
     1455          ENDDO
    14371456       ENDDO
    14381457
     
    14441463
    14451464    ENDIF
    1446 !
    1447 !-- Store the old age of the particle ( needed to prevent that a
    1448 !-- particle crosses several PEs during one timestep, and for the
    1449 !-- evaluation of the subgrid particle velocity fluctuations )
    1450     particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
    14511465
    14521466    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
    14531467
    14541468    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
    1455        DO  n = 1, number_of_particles
    1456 
    1457 !
    1458 !--       Particle advection
    1459           IF ( dens_ratio(n) == 0.0_wp )  THEN
    1460 !
    1461 !--          Pure passive transport (without particle inertia)
    1462              particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
    1463              particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
    1464              particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
    1465 
    1466              particles(n)%speed_x = u_int(n)
    1467              particles(n)%speed_y = v_int(n)
    1468              particles(n)%speed_z = w_int(n)
    1469 
    1470           ELSE
     1469       DO  nb = 0, 7
     1470          DO  n = start_index(nb), end_index(nb)
     1471
     1472!
     1473!--          Particle advection
     1474             IF ( dens_ratio(n) == 0.0_wp )  THEN
     1475!
     1476!--             Pure passive transport (without particle inertia)
     1477                particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
     1478                particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
     1479                particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
     1480
     1481                particles(n)%speed_x = u_int(n)
     1482                particles(n)%speed_y = v_int(n)
     1483                particles(n)%speed_z = w_int(n)
     1484
     1485             ELSE
     1486!
     1487!--             Transport of particles with inertia
     1488                particles(n)%x = particles(n)%x + particles(n)%speed_x * &
     1489                                                  dt_particle(n)
     1490                particles(n)%y = particles(n)%y + particles(n)%speed_y * &
     1491                                                  dt_particle(n)
     1492                particles(n)%z = particles(n)%z + particles(n)%speed_z * &
     1493                                                  dt_particle(n)
     1494
     1495!
     1496!--             Update of the particle velocity
     1497                IF ( cloud_droplets )  THEN
     1498!
     1499!--                Terminal velocity is computed for vertical direction (Rogers et
     1500!--                al., 1993, J. Appl. Meteorol.)
     1501                   diameter = particles(n)%radius * 2000.0_wp !diameter in mm
     1502                   IF ( diameter <= d0_rog )  THEN
     1503                      w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
     1504                   ELSE
     1505                      w_s = a_rog - b_rog * EXP( -c_rog * diameter )
     1506                   ENDIF
     1507
     1508!
     1509!--                If selected, add random velocities following Soelch and Kaercher
     1510!--                (2010, Q. J. R. Meteorol. Soc.)
     1511                   IF ( use_sgs_for_particles )  THEN
     1512                      lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
     1513                      RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
     1514                      sigma          = SQRT( e(kp,jp,ip) )
     1515
     1516                      rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1517                      rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1518                      rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1519
     1520                      particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
     1521                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg1
     1522                      particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
     1523                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg2
     1524                      particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
     1525                                           SQRT( 1.0_wp - RL**2 ) * sigma * rg3
     1526
     1527                      particles(n)%speed_x = u_int(n) + particles(n)%rvar1
     1528                      particles(n)%speed_y = v_int(n) + particles(n)%rvar2
     1529                      particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
     1530                   ELSE
     1531                      particles(n)%speed_x = u_int(n)
     1532                      particles(n)%speed_y = v_int(n)
     1533                      particles(n)%speed_z = w_int(n) - w_s
     1534                   ENDIF
     1535
     1536                ELSE
     1537
     1538                   IF ( use_sgs_for_particles )  THEN
     1539                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1540                      exp_term = EXP( -exp_arg * dt_particle(n) )
     1541                   ELSE
     1542                      exp_arg  = particle_groups(particles(n)%group)%exp_arg
     1543                      exp_term = particle_groups(particles(n)%group)%exp_term
     1544                   ENDIF
     1545                   particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
     1546                                          u_int(n) * ( 1.0_wp - exp_term )
     1547                   particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
     1548                                          v_int(n) * ( 1.0_wp - exp_term )
     1549                   particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
     1550                                          ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
     1551                                          g / exp_arg ) * ( 1.0_wp - exp_term )
     1552                ENDIF
     1553
     1554             ENDIF
     1555          ENDDO
     1556       ENDDO
     1557   
     1558    ELSE
     1559
     1560       DO  nb = 0, 7
     1561          DO  n = start_index(nb), end_index(nb)
    14711562!
    14721563!--          Transport of particles with inertia
    1473              particles(n)%x = particles(n)%x + particles(n)%speed_x * &
    1474                                                dt_particle(n)
    1475              particles(n)%y = particles(n)%y + particles(n)%speed_y * &
    1476                                                dt_particle(n)
    1477              particles(n)%z = particles(n)%z + particles(n)%speed_z * &
    1478                                                dt_particle(n)
    1479 
     1564             particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
     1565             particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
     1566             particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
    14801567!
    14811568!--          Update of the particle velocity
    14821569             IF ( cloud_droplets )  THEN
    14831570!
    1484 !--             Terminal velocity is computed for vertical direction (Rogers et
    1485 !--             al., 1993, J. Appl. Meteorol.)
     1571!--             Terminal velocity is computed for vertical direction (Rogers et al.,
     1572!--             1993, J. Appl. Meteorol.)
    14861573                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
    14871574                IF ( diameter <= d0_rog )  THEN
     
    14951582!--             (2010, Q. J. R. Meteorol. Soc.)
    14961583                IF ( use_sgs_for_particles )  THEN
    1497                    lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
    1498                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
    1499                    sigma          = SQRT( e(kp,jp,ip) )
    1500 
    1501                    rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1502                    rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1503                    rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1504 
    1505                    particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
    1506                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg1
    1507                    particles(n)%rvar2 = RL * particles(n)%rvar2 +              &
    1508                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg2
    1509                    particles(n)%rvar3 = RL * particles(n)%rvar3 +              &
    1510                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg3
    1511 
    1512                    particles(n)%speed_x = u_int(n) + particles(n)%rvar1
    1513                    particles(n)%speed_y = v_int(n) + particles(n)%rvar2
    1514                    particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
     1584                    lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
     1585                    RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
     1586                    sigma          = SQRT( e(kp,jp,ip) )
     1587
     1588                    rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1589                    rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1590                    rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     1591
     1592                    particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
     1593                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg1
     1594                    particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
     1595                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg2
     1596                    particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
     1597                                         SQRT( 1.0_wp - RL**2 ) * sigma * rg3
     1598
     1599                    particles(n)%speed_x = u_int(n) + particles(n)%rvar1
     1600                    particles(n)%speed_y = v_int(n) + particles(n)%rvar2
     1601                    particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
    15151602                ELSE
    1516                    particles(n)%speed_x = u_int(n)
    1517                    particles(n)%speed_y = v_int(n)
    1518                    particles(n)%speed_z = w_int(n) - w_s
     1603                    particles(n)%speed_x = u_int(n)
     1604                    particles(n)%speed_y = v_int(n)
     1605                    particles(n)%speed_z = w_int(n) - w_s
    15191606                ENDIF
    15201607
     
    15281615                   exp_term = particle_groups(particles(n)%group)%exp_term
    15291616                ENDIF
    1530                 particles(n)%speed_x = particles(n)%speed_x * exp_term +         &
     1617                particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
    15311618                                       u_int(n) * ( 1.0_wp - exp_term )
    1532                 particles(n)%speed_y = particles(n)%speed_y * exp_term +         &
     1619                particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
    15331620                                       v_int(n) * ( 1.0_wp - exp_term )
    1534                 particles(n)%speed_z = particles(n)%speed_z * exp_term +         &
    1535                                        ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
    1536                                        g / exp_arg ) * ( 1.0_wp - exp_term )
     1621                particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
     1622                                       ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
     1623                                       exp_arg ) * ( 1.0_wp - exp_term )
    15371624             ENDIF
    1538 
     1625          ENDDO
     1626       ENDDO
     1627
     1628    ENDIF
     1629
     1630!
     1631!-- Store the old age of the particle ( needed to prevent that a
     1632!-- particle crosses several PEs during one timestep, and for the
     1633!-- evaluation of the subgrid particle velocity fluctuations )
     1634    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
     1635
     1636    DO  nb = 0, 7
     1637       DO  n = start_index(nb), end_index(nb)
     1638!
     1639!--       Increment the particle age and the total time that the particle
     1640!--       has advanced within the particle timestep procedure
     1641          particles(n)%age    = particles(n)%age    + dt_particle(n)
     1642          particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
     1643
     1644!
     1645!--       Check whether there is still a particle that has not yet completed
     1646!--       the total LES timestep
     1647          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
     1648             dt_3d_reached_l = .FALSE.
    15391649          ENDIF
    15401650
    15411651       ENDDO
    1542    
    1543     ELSE
    1544 
    1545        DO  n = 1, number_of_particles
    1546 
    1547 !--       Transport of particles with inertia
    1548           particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
    1549           particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
    1550           particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
    1551 !
    1552 !--       Update of the particle velocity
    1553           IF ( cloud_droplets )  THEN
    1554 !
    1555 !--          Terminal velocity is computed for vertical direction (Rogers et al.,
    1556 !--          1993, J. Appl. Meteorol.)
    1557              diameter = particles(n)%radius * 2000.0_wp !diameter in mm
    1558              IF ( diameter <= d0_rog )  THEN
    1559                 w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
    1560              ELSE
    1561                 w_s = a_rog - b_rog * EXP( -c_rog * diameter )
    1562              ENDIF
    1563 
    1564 !
    1565 !--          If selected, add random velocities following Soelch and Kaercher
    1566 !--          (2010, Q. J. R. Meteorol. Soc.)
    1567              IF ( use_sgs_for_particles )  THEN
    1568                  lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
    1569                  RL             = EXP( -1.0_wp * dt_3d / lagr_timescale )
    1570                  sigma          = SQRT( e(kp,jp,ip) )
    1571 
    1572                  rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1573                  rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1574                  rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    1575 
    1576                  particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
    1577                                       SQRT( 1.0_wp - RL**2 ) * sigma * rg1
    1578                  particles(n)%rvar2 = RL * particles(n)%rvar2 +                &
    1579                                       SQRT( 1.0_wp - RL**2 ) * sigma * rg2
    1580                  particles(n)%rvar3 = RL * particles(n)%rvar3 +                &
    1581                                       SQRT( 1.0_wp - RL**2 ) * sigma * rg3
    1582 
    1583                  particles(n)%speed_x = u_int(n) + particles(n)%rvar1
    1584                  particles(n)%speed_y = v_int(n) + particles(n)%rvar2
    1585                  particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
    1586              ELSE
    1587                  particles(n)%speed_x = u_int(n)
    1588                  particles(n)%speed_y = v_int(n)
    1589                  particles(n)%speed_z = w_int(n) - w_s
    1590              ENDIF
    1591 
    1592           ELSE
    1593 
    1594              IF ( use_sgs_for_particles )  THEN
    1595                 exp_arg  = particle_groups(particles(n)%group)%exp_arg
    1596                 exp_term = EXP( -exp_arg * dt_particle(n) )
    1597              ELSE
    1598                 exp_arg  = particle_groups(particles(n)%group)%exp_arg
    1599                 exp_term = particle_groups(particles(n)%group)%exp_term
    1600              ENDIF
    1601              particles(n)%speed_x = particles(n)%speed_x * exp_term +             &
    1602                                     u_int(n) * ( 1.0_wp - exp_term )
    1603              particles(n)%speed_y = particles(n)%speed_y * exp_term +             &
    1604                                     v_int(n) * ( 1.0_wp - exp_term )
    1605              particles(n)%speed_z = particles(n)%speed_z * exp_term +             &
    1606                                     ( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
    1607                                     exp_arg ) * ( 1.0_wp - exp_term )
    1608           ENDIF
    1609 
    1610        ENDDO
    1611 
    1612     ENDIF
    1613 
    1614     DO  n = 1, number_of_particles
    1615 !
    1616 !--    Increment the particle age and the total time that the particle
    1617 !--    has advanced within the particle timestep procedure
    1618        particles(n)%age    = particles(n)%age    + dt_particle(n)
    1619        particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
    1620 
    1621 !
    1622 !--    Check whether there is still a particle that has not yet completed
    1623 !--    the total LES timestep
    1624        IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
    1625           dt_3d_reached_l = .FALSE.
    1626        ENDIF
    1627 
    16281652    ENDDO
    16291653
     
    16931717!-- subgrid-scale velocity component is set to zero, in order to prevent a
    16941718!-- velocity build-up.
    1695 
    16961719!-- This case, set also previous subgrid-scale component to zero.
    16971720    v_sgs = v_sgs * fac + term1 + term2 + term3
  • palm/trunk/SOURCE/lpm_pack_arrays.f90

    r2101 r2417  
    2525! -----------------
    2626! $Id$
     27! New routine which sorts particles into particles that completed and not
     28! completed the LES timestep.
     29!
     30! 2101 2017-01-05 16:42:31Z suehring
    2731!
    2832! 2000 2016-08-20 18:09:15Z knoop
     
    6165!> deletion and move data with higher index values to these free indices.
    6266!> Determine the new number of particles.
     67!> Moreover, particles are also sorted into groups finished and not finished
     68!> its timestep.
    6369!------------------------------------------------------------------------------!
    6470 MODULE lpm_pack_arrays_mod
     
    7076
    7177    PRIVATE
    72     PUBLIC lpm_pack_all_arrays, lpm_pack_arrays
     78    PUBLIC lpm_pack_all_arrays, lpm_pack_arrays, lpm_sort
    7379
    7480    INTERFACE lpm_pack_all_arrays
     
    8086    END INTERFACE lpm_pack_arrays
    8187
    82 CONTAINS
     88    INTERFACE lpm_sort
     89       MODULE PROCEDURE lpm_sort
     90    END INTERFACE lpm_sort
     91
     92
     93 CONTAINS
    8394
    8495!------------------------------------------------------------------------------!
     
    236247    END SUBROUTINE lpm_pack_and_sort
    237248
     249!------------------------------------------------------------------------------!
     250! Description:
     251! ------------
     252!> Sort particles in each sub-grid box into two groups: particles that already
     253!> completed the LES timestep, and particles that need further timestepping to
     254!> complete the LES timestep.
     255!------------------------------------------------------------------------------!
     256    SUBROUTINE lpm_sort
     257
     258       USE control_parameters,                                                 &
     259           ONLY:  dt_3d
     260
     261       USE indices,                                                            &
     262           ONLY: nxl, nxr, nys, nyn, nzb, nzt
     263
     264       USE kinds
     265
     266       IMPLICIT NONE
     267
     268       INTEGER(iwp) :: end_index     !< particle end index for each sub-box
     269       INTEGER(iwp) :: i             !< index of particle grid box in x-direction
     270       INTEGER(iwp) :: j             !< index of particle grid box in y-direction
     271       INTEGER(iwp) :: k             !< index of particle grid box in z-direction
     272       INTEGER(iwp) :: n             !< running index for number of particles
     273       INTEGER(iwp) :: nb            !< index of subgrid boux
     274       INTEGER(iwp) :: nf            !< indices for particles in each sub-box that already finalized their substeps
     275       INTEGER(iwp) :: nnf           !< indices for particles in each sub-box that need further treatment
     276       INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps
     277       INTEGER(iwp) :: start_index   !< particle start index for each sub-box
     278
     279       TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles  !< temporary particle array
     280
     281       DO  i = nxl, nxr
     282          DO  j = nys, nyn
     283             DO  k = nzb+1, nzt
     284
     285                number_of_particles = prt_count(k,j,i)
     286                IF ( number_of_particles <= 0 )  CYCLE
     287
     288                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     289
     290                DO  nb = 0, 7
     291!
     292!--                Obtain start and end index for each subgrid box
     293                   start_index = grid_particles(k,j,i)%start_index(nb)
     294                   end_index   = grid_particles(k,j,i)%end_index(nb)
     295!
     296!--                Allocate temporary array used for sorting.
     297                   ALLOCATE( sort_particles(start_index:end_index) )
     298!
     299!--                Determine number of particles already completed the LES
     300!--                timestep, and write them into a temporary array.
     301                   nf = start_index
     302                   num_finalized = 0
     303                   DO  n = start_index, end_index
     304                      IF ( dt_3d - particles(n)%dt_sum < 1E-8_wp )  THEN
     305                         sort_particles(nf) = particles(n)
     306                         nf                 = nf + 1
     307                         num_finalized      = num_finalized + 1
     308                      ENDIF
     309                   ENDDO
     310!
     311!--                Determine number of particles that not completed the LES
     312!--                timestep, and write them into a temporary array.
     313                   nnf = nf
     314                   DO  n = start_index, end_index
     315                      IF ( dt_3d - particles(n)%dt_sum > 1E-8_wp )  THEN
     316                         sort_particles(nnf) = particles(n)
     317                         nnf                 = nnf + 1
     318                      ENDIF
     319                   ENDDO
     320!
     321!--                Write back sorted particles
     322                   particles(start_index:end_index) =                          &
     323                                           sort_particles(start_index:end_index)
     324!
     325!--                Determine updated start_index, used to masked already
     326!--                completed particles.
     327                   grid_particles(k,j,i)%start_index(nb) =                     &
     328                                      grid_particles(k,j,i)%start_index(nb)    &
     329                                    + num_finalized
     330!
     331!--                Deallocate dummy array
     332                   DEALLOCATE ( sort_particles )
     333!
     334!--                Finally, if number of non-completed particles is non zero
     335!--                in any of the sub-boxes, set control flag appropriately.
     336                   IF ( nnf > nf )                                             &
     337                      grid_particles(k,j,i)%time_loop_done = .FALSE.
     338
     339                ENDDO
     340             ENDDO
     341          ENDDO
     342       ENDDO
     343
     344    END SUBROUTINE lpm_sort
     345
    238346
    239347 END module lpm_pack_arrays_mod
  • palm/trunk/SOURCE/user_lpm_advec.f90

    r2101 r2417  
    2525! -----------------
    2626! $Id$
     27! Particle loops adapted for sub-box structure, i.e. for each sub-box the
     28! particle loop runs from start_index up to end_index instead from 1 to
     29! number_of_particles.
     30!
     31! 2101 2017-01-05 16:42:31Z suehring
    2732!
    2833! 2000 2016-08-20 18:09:15Z knoop
     
    7075    IMPLICIT NONE
    7176
    72     INTEGER(iwp) ::  ip   !<
    73     INTEGER(iwp) ::  jp   !<
    74     INTEGER(iwp) ::  kp   !<
    75     INTEGER(iwp) ::  n    !<
     77    INTEGER(iwp) ::  ip   !< index of particle grid box, x-direction
     78    INTEGER(iwp) ::  jp   !< index of particle grid box, y-direction
     79    INTEGER(iwp) ::  kp   !< index of particle grid box, z-direction
     80    INTEGER(iwp) ::  n    !< particle index
     81    INTEGER(iwp) ::  nb   !< index of sub-box particles are sorted in
     82
     83    INTEGER(iwp), DIMENSION(0:7)  ::  start_index !< start particle index for current sub-box
     84    INTEGER(iwp), DIMENSION(0:7)  ::  end_index   !< start particle index for current sub-box
    7685
    7786!
     
    7988!   number_of_particles = prt_count(kp,jp,ip)
    8089!   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     90!
     91!   start_index = grid_particles(kp,jp,ip)%start_index
     92!   end_index   = grid_particles(kp,jp,ip)%end_index
     93!
    8194!   IF ( number_of_particles <= 0 )  CYCLE
    8295!   DO  n = 1, number_of_particles
    83 !
     96!   DO  nb = 0, 7
     97!      DO  n = start_index(nb), end_index(nb)
     98!         particles(n)%xxx =
     99!      ENDDO
    84100!   ENDDO
    85101
Note: See TracChangeset for help on using the changeset viewer.