Ignore:
Timestamp:
Jul 26, 2019 10:01:22 AM (5 years ago)
Author:
schwenkel
Message:

Implementation of an simple method for interpolating the velocities to particle position

File:
1 edited

Legend:

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

    r4114 r4121  
    2525! -----------------
    2626! $Id$
     27! Implementation of an simple method for interpolating the velocities to
     28! particle position
     29!
     30! 4114 2019-07-23 14:09:27Z schwenkel
    2731! Bugfix: Added working precision for if statement
    2832!
     
    7276! with scalar grid point of same index.
    7377! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
    74 ! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
     78! lpm_pack_all_arrays -> lpm_sort_and_delete, lpm_pack_arrays -> lpm_pack
    7579! lpm_sort -> lpm_sort_timeloop_done
    7680!
     
    169173! Description:
    170174! ------------
    171 !> 
     175!>
    172176!------------------------------------------------------------------------------!
    173177 MODULE lagrangian_particle_model_mod
     
    178182        ONLY:  de_dx, de_dy, de_dz, dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,      &
    179183               pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p,  &
    180                d_exner
     184               d_exner, u_p, v_p, w_p
    181185 
    182186    USE averaging,                                                             &
     
    263267    IMPLICIT NONE
    264268
    265     CHARACTER(LEN=15) ::  aero_species = 'nacl'                    !< aerosol species
    266     CHARACTER(LEN=15) ::  aero_type    = 'maritime'                !< aerosol type
    267     CHARACTER(LEN=15) ::  bc_par_lr    = 'cyclic'                  !< left/right boundary condition
    268     CHARACTER(LEN=15) ::  bc_par_ns    = 'cyclic'                  !< north/south boundary condition
    269     CHARACTER(LEN=15) ::  bc_par_b     = 'reflect'                 !< bottom boundary condition
    270     CHARACTER(LEN=15) ::  bc_par_t     = 'absorb'                  !< top boundary condition
    271     CHARACTER(LEN=15) ::  collision_kernel   = 'none'              !< collision kernel   
    272 
    273     CHARACTER(LEN=5)  ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
    274     CHARACTER(LEN=5)  ::  splitting_mode     = 'const'             !< splitting mode
     269    CHARACTER(LEN=15) ::  aero_species = 'nacl'                   !< aerosol species
     270    CHARACTER(LEN=15) ::  aero_type    = 'maritime'               !< aerosol type
     271    CHARACTER(LEN=15) ::  bc_par_lr    = 'cyclic'                 !< left/right boundary condition
     272    CHARACTER(LEN=15) ::  bc_par_ns    = 'cyclic'                 !< north/south boundary condition
     273    CHARACTER(LEN=15) ::  bc_par_b     = 'reflect'                !< bottom boundary condition
     274    CHARACTER(LEN=15) ::  bc_par_t     = 'absorb'                 !< top boundary condition
     275    CHARACTER(LEN=15) ::  collision_kernel   = 'none'             !< collision kernel
     276
     277    CHARACTER(LEN=5)  ::  splitting_function = 'gamma'            !< function for calculation critical weighting factor
     278    CHARACTER(LEN=5)  ::  splitting_mode     = 'const'            !< splitting mode
     279
     280    CHARACTER(LEN=25) ::  particle_interpolation = 'trilinear'    !< interpolation method for calculatin the particle
    275281
    276282    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step   
     
    496502    END INTERFACE dealloc_particles_array
    497503
    498     INTERFACE lpm_sort_in_subboxes
    499        MODULE PROCEDURE lpm_sort_in_subboxes
    500     END INTERFACE lpm_sort_in_subboxes
     504    INTERFACE lpm_sort_and_delete
     505       MODULE PROCEDURE lpm_sort_and_delete
     506    END INTERFACE lpm_sort_and_delete
    501507
    502508    INTERFACE lpm_sort_timeloop_done
     
    549555       particles_per_point, &
    550556       particle_advection_start, &
     557       particle_interpolation, &
    551558       particle_maximum_age, &
    552559       pdx, &
     
    608615       particles_per_point, &
    609616       particle_advection_start, &
     617       particle_interpolation, &
    610618       particle_maximum_age, &
    611619       pdx, &
     
    849857    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
    850858
     859!
     860!-- Subgrid scale velocites with the simple interpolation method for resolved
     861!-- velocites is not implemented for passive particles. However, for cloud
     862!-- it can be combined as the sgs-velocites for active particles are
     863!-- calculated differently, i.e. no subboxes are needed.
     864    IF ( .NOT. TRIM(particle_interpolation) == 'trilinear'  .AND.              &
     865       use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
     866          message_string = 'subrgrid scale velocities in combination with ' // &
     867                           'simple interpolation method is not '            // &
     868                           'implemented'
     869          CALL message( 'check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
     870    ENDIF
     871
    851872 END SUBROUTINE
    852873 
     
    10671088       CASE ( 'reflect' )
    10681089          ibc_par_t = 2
    1069          
     1090
    10701091       CASE ( 'nested' )
    10711092          ibc_par_t = 3
     
    10871108       CASE ( 'reflect' )
    10881109          ibc_par_lr = 2
    1089          
     1110
    10901111       CASE ( 'nested' )
    10911112          ibc_par_lr = 3
     
    11071128       CASE ( 'reflect' )
    11081129          ibc_par_ns = 2
    1109          
     1130
    11101131       CASE ( 'nested' )
    11111132          ibc_par_ns = 3
     
    15461567!-- which is needed for a fast interpolation of the LES fields on the particle
    15471568!-- position.
    1548     CALL lpm_sort_in_subboxes
     1569    CALL lpm_sort_and_delete
    15491570!
    15501571!-- Determine the current number of particles
     
    16231644!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
    16241645    IF ( TRIM(aero_type) .EQ. 'polar' )  THEN
    1625        na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6
    1626        rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6
     1646       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp
     1647       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp
    16271648       log_sigma = (/ 0.245, 0.300, 0.291 /)
    16281649    ELSEIF ( TRIM(aero_type) .EQ. 'background' )  THEN
    1629        na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6
    1630        rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6
     1650       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp
     1651       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp
    16311652       log_sigma = (/ 0.645, 0.253, 0.425 /)
    16321653    ELSEIF ( TRIM(aero_type) .EQ. 'maritime' )  THEN
    1633        na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6
    1634        rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6
     1654       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp
     1655       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp
    16351656       log_sigma = (/ 0.657, 0.210, 0.396 /)
    16361657    ELSEIF ( TRIM(aero_type) .EQ. 'continental' )  THEN
    1637        na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6
    1638        rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6
     1658       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp
     1659       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp
    16391660       log_sigma = (/ 0.161, 0.217, 0.380 /)
    16401661    ELSEIF ( TRIM(aero_type) .EQ. 'desert' )  THEN
    1641        na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6
    1642        rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6
     1662       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp
     1663       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp
    16431664       log_sigma = (/ 0.247, 0.770, 0.438 /)
    16441665    ELSEIF ( TRIM(aero_type) .EQ. 'rural' )  THEN
    1645        na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6
    1646        rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6
     1666       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp
     1667       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp
    16471668       log_sigma = (/ 0.225, 0.557, 0.266 /)
    16481669    ELSEIF ( TRIM(aero_type) .EQ. 'urban' )  THEN
    1649        na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6
    1650        rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6
     1670       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp
     1671       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp
    16511672       log_sigma = (/ 0.245, 0.666, 0.337 /)
    16521673    ELSEIF ( TRIM(aero_type) .EQ. 'user' )  THEN
     
    16791700                particles(n)%aux1          = r_mid
    16801701                particles(n)%weight_factor =                                           &
    1681                    ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) *                     &
    1682                      EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) +  &
    1683                      na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) *                     &
    1684                      EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) +  &
    1685                      na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) *                     &
    1686                      EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) )    &
     1702                   ( na(1) / ( SQRT( 2.0_wp * pi ) * log_sigma(1) ) *                     &
     1703                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0_wp * log_sigma(1)**2 ) ) +  &
     1704                     na(2) / ( SQRT( 2.0_wp * pi ) * log_sigma(2) ) *                     &
     1705                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0_wp * log_sigma(2)**2 ) ) +  &
     1706                     na(3) / ( SQRT( 2.0_wp * pi ) * log_sigma(3) ) *                     &
     1707                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0_wp * log_sigma(3)**2 ) )    &
    16871708                   ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dzw(kp) )
    16881709
     
    17001721!
    17011722!--             Unnecessary particles will be deleted
    1702                 IF ( particles(n)%weight_factor .LE. 0.0 )  particles(n)%particle_mask = .FALSE.
     1723                IF ( particles(n)%weight_factor .LE. 0.0_wp )  particles(n)%particle_mask = .FALSE.
    17031724
    17041725             ENDDO
     
    17411762
    17421763 END SUBROUTINE lpm_init_aerosols
    1743  
    1744  
     1764
     1765
    17451766!------------------------------------------------------------------------------!
    17461767! Description:
     
    17521773!------------------------------------------------------------------------------!
    17531774 SUBROUTINE lpm_init_sgs_tke
    1754    
     1775
    17551776    USE statistics,                                                            &
    17561777        ONLY:  flow_statistics_called, hom, sums, sums_l
     
    20132034    LOGICAL            ::  first_loop_stride  !<
    20142035
    2015    
     2036
    20162037    SELECT CASE ( location )
    20172038
     
    20982119             CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
    20992120             CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
    2100              
     2121
    21012122!
    21022123!--          If particle advection includes SGS velocity components, calculate the
     
    21162137             IF ( .NOT. first_loop_stride  .AND.  use_sgs_for_particles )            &
    21172138                CALL lpm_sort_timeloop_done
    2118 
    21192139             DO  i = nxl, nxr
    21202140                DO  j = nys, nyn
     
    21642184!
    21652185!--                   Particle advection
    2166                       CALL lpm_advec(i,j,k)
     2186                      CALL lpm_advec(TRIM(particle_interpolation),i,j,k)
    21672187!
    21682188!--                   Particle reflection from walls. Only applied if the particles
     
    22032223                ENDDO
    22042224             ENDDO
    2205 
    22062225             steps = steps + 1
    22072226             dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
     
    22162235             dt_3d_reached = dt_3d_reached_l
    22172236#endif
    2218 
    22192237             CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
    22202238
     
    22372255
    22382256!
    2239 !--          IF .FALSE., lpm_sort_in_subboxes is done inside pcmp
     2257!--          IF .FALSE., lpm_sort_and_delete is done inside pcmp
    22402258             IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run )   THEN   
    22412259!
    22422260!--             Pack particles (eliminate those marked for deletion),
    22432261!--             determine new number of particles
    2244                 CALL lpm_sort_in_subboxes
    2245 !
     2262                CALL lpm_sort_and_delete
     2263
    22462264!--             Initialize variables for the next (sub-) timestep, i.e., for marking
    22472265!--             those particles to be deleted after the timestep
     
    22532271             first_loop_stride = .FALSE.
    22542272          ENDDO   ! timestep loop
    2255 !   
     2273!
    22562274!--       in case of nested runs do the transfer of particles after every full model time step
    22572275          IF ( nested_run )   THEN
     
    22602278             CALL pmcp_p_delete_particles_in_fine_grid_area
    22612279
    2262              CALL lpm_sort_in_subboxes
     2280             CALL lpm_sort_and_delete
    22632281
    22642282             deleted_particles = 0
     
    22842302
    22852303          CALL cpu_log( log_point(25), 'lpm', 'stop' )
    2286          
     2304
    22872305! !
    22882306! !--       Output of particle time series
     
    22952313!              ENDIF
    22962314!           ENDIF
    2297    
     2315
    22982316       CASE DEFAULT
    22992317          CONTINUE
     
    29402958!-- Must be called to sort particles into blocks, which is needed for a fast
    29412959!-- interpolation of the LES fields on the particle position.
    2942     CALL lpm_sort_in_subboxes
    2943        
     2960    CALL lpm_sort_and_delete
     2961
    29442962
    29452963 END SUBROUTINE lpm_rrd_local_particles
     
    29492967                              nxr_on_file, nynf, nync, nyn_on_file, nysf,  &
    29502968                              nysc, nys_on_file, tmp_3d, found )
    2951                              
    2952        
     2969
     2970
    29532971   USE control_parameters,                                                 &
    2954        ONLY: length, restart_string           
     2972       ONLY: length, restart_string
    29552973
    29562974    INTEGER(iwp) ::  k               !<
     
    29722990    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
    29732991
    2974    
     2992
    29752993    found = .TRUE.
    29762994
    29772995    SELECT CASE ( restart_string(1:length) )
    2978    
     2996
    29792997       CASE ( 'iran' ) ! matching random numbers is still unresolved issue
    29802998          IF ( k == 1 )  READ ( 13 )  iran, iran_part
    2981          
     2999
    29823000        CASE ( 'pc_av' )
    29833001           IF ( .NOT. ALLOCATED( pc_av ) )  THEN
     
    31543172!> grid scale) as well as the boundary conditions of particles. As a next step
    31553173!> this submodule should be excluded as an own file.
    3156 !------------------------------------------------------------------------------!   
    3157  SUBROUTINE lpm_advec (ip,jp,kp)
    3158 
     3174!------------------------------------------------------------------------------!
     3175 SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp)
     3176
     3177    CHARACTER (LEN=*), INTENT(IN) ::  interpolation_method !<
    31593178    LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
    31603179
    31613180    INTEGER(iwp) ::  i                           !< index variable along x
    3162     INTEGER(iwp) ::  ip                          !< index variable along x
    3163     INTEGER(iwp) ::  j                           !< index variable along y
    3164     INTEGER(iwp) ::  jp                          !< index variable along y
    3165     INTEGER(iwp) ::  k                           !< index variable along z
     3181    INTEGER(iwp) ::  i_next                      !< index variable along x
     3182    INTEGER(iwp) ::  ip                          !< index variable along x
     3183    INTEGER(iwp) ::  iteration_steps = 1         !< amount of iterations steps for corrector step
     3184    INTEGER(iwp) ::  j                           !< index variable along y
     3185    INTEGER(iwp) ::  j_next                      !< index variable along y
     3186    INTEGER(iwp) ::  jp                          !< index variable along y
     3187    INTEGER(iwp) ::  k                           !< index variable along z
    31663188    INTEGER(iwp) ::  k_wall                      !< vertical index of topography top
    3167     INTEGER(iwp) ::  kp                          !< index variable along z
     3189    INTEGER(iwp) ::  kp                          !< index variable along z
     3190    INTEGER(iwp) ::  k_next                      !< index variable along z
    31683191    INTEGER(iwp) ::  kw                          !< index variable along z
     3192    INTEGER(iwp) ::  kkw                         !< index variable along z
    31693193    INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
    31703194    INTEGER(iwp) ::  nb                          !< block number particles are sorted in
    3171     INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box
     3195    INTEGER(iwp) ::  particle_end                !< end index for partilce loop
     3196    INTEGER(iwp) ::  particle_start              !< start index for particle loop
     3197    INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box
     3198    INTEGER(iwp) ::  subbox_end                  !< end index for loop over subboxes in particle advection
     3199    INTEGER(iwp) ::  subbox_start                !< start index for loop over subboxes in particle advection
     3200    INTEGER(iwp) ::  nn                          !< loop variable over iterations steps
    31723201
    31733202    INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
    31743203    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
    31753204
    3176     REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
    3177     REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
    3178     REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
     3205    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
     3206    REAL(wp) ::  alpha              !< interpolation facor for x-direction
     3207
     3208    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
     3209    REAL(wp) ::  beta               !< interpolation facor for y-direction
     3210    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
    31793211    REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation
    31803212    REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation
     
    31973229    REAL(wp) ::  exp_arg            !< argument in the exponent - particle radius
    31983230    REAL(wp) ::  exp_term           !< exponent term
     3231    REAL(wp) ::  gamma              !< interpolation facor for z-direction
    31993232    REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation
    32003233    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
     
    32083241    REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
    32093242    REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
     3243    REAL(wp) ::  unext              !< calculated particle u-velocity of corrector step
    32103244    REAL(wp) ::  us_int             !< friction velocity at particle grid box
    32113245    REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
     
    32133247    REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
    32143248    REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
     3249    REAL(wp) ::  vnext              !< calculated particle v-velocity of corrector step
    32153250    REAL(wp) ::  vv_int             !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection
    32163251    REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
    32173252    REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
     3253    REAL(wp) ::  wnext              !< calculated particle w-velocity of corrector step
    32183254    REAL(wp) ::  w_s                !< terminal velocity of droplets
    32193255    REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation
    3220     REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
     3256    REAL(wp) ::  xp                 !< calculated particle position in x of predictor step
     3257    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
     3258    REAL(wp) ::  yp                 !< calculated particle position in y of predictor step
    32213259    REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
     3260    REAL(wp) ::  zp                 !< calculated particle position in z of predictor step
    32223261
    32233262    REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
     
    32493288    REAL(wp), DIMENSION(number_of_particles) ::  zv             !< z-position
    32503289
    3251     REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers 
     3290    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
    32523291
    32533292    CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
    3254 
    32553293!
    32563294!-- Determine height of Prandtl layer and distance between Prandtl-layer
    3257 !-- height and horizontal mean roughness height, which are required for 
    3258 !-- vertical logarithmic interpolation of horizontal particle speeds 
     3295!-- height and horizontal mean roughness height, which are required for
     3296!-- vertical logarithmic interpolation of horizontal particle speeds
    32593297!-- (for particles below first vertical grid level).
    32603298    z_p      = zu(nzb+1) - zw(nzb)
    32613299    d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
    32623300
    3263     start_index = grid_particles(kp,jp,ip)%start_index
    3264     end_index   = grid_particles(kp,jp,ip)%end_index
    3265 
    32663301    xv = particles(1:number_of_particles)%x
    32673302    yv = particles(1:number_of_particles)%y
    32683303    zv = particles(1:number_of_particles)%z
    3269 
    3270     DO  nb = 0, 7
    3271 !
    3272 !--    Interpolate u velocity-component       
    3273        i = ip
    3274        j = jp + block_offset(nb)%j_off
    3275        k = kp + block_offset(nb)%k_off
    3276 
    3277        DO  n = start_index(nb), end_index(nb)
    3278 !
    3279 !--       Interpolation of the u velocity component onto particle position. 
    3280 !--       Particles are interpolation bi-linearly in the horizontal and a
    3281 !--       linearly in the vertical. An exception is made for particles below
    3282 !--       the first vertical grid level in case of a prandtl layer. In this
    3283 !--       case the horizontal particle velocity components are determined using
    3284 !--       Monin-Obukhov relations (if branch).
    3285 !--       First, check if particle is located below first vertical grid level
    3286 !--       above topography (Prandtl-layer height)
    3287 !--       Determine vertical index of topography top
    3288           k_wall = get_topography_top_index_ji( jp, ip, 's' )
    3289 
    3290           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
    3291 !
    3292 !--          Resolved-scale horizontal particle velocity is zero below z0.
    3293              IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
    3294                 u_int(n) = 0.0_wp
    3295              ELSE
    3296 !
    3297 !--             Determine the sublayer. Further used as index.
    3298                 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
    3299                                      * REAL( number_of_sublayers, KIND=wp )    &
    3300                                      * d_z_p_z0
    3301 !
    3302 !--             Calculate LOG(z/z0) for exact particle height. Therefore,   
    3303 !--             interpolate linearly between precalculated logarithm.
    3304                 log_z_z0_int = log_z_z0(INT(height_p))                         &
    3305                                  + ( height_p - INT(height_p) )                &
    3306                                  * ( log_z_z0(INT(height_p)+1)                 &
    3307                                       - log_z_z0(INT(height_p))                &
    3308                                    )
    3309 !
    3310 !--             Get friction velocity and momentum flux from new surface data
    3311 !--             types.
    3312                 IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
    3313                      surf_def_h(0)%end_index(jp,ip) )  THEN
    3314                    surf_start = surf_def_h(0)%start_index(jp,ip)
    3315 !--                Limit friction velocity. In narrow canyons or holes the
    3316 !--                friction velocity can become very small, resulting in a too
    3317 !--                large particle speed.
    3318                    us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
    3319                    usws_int  = surf_def_h(0)%usws(surf_start)
    3320                 ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
    3321                          surf_lsm_h%end_index(jp,ip) )  THEN
    3322                    surf_start = surf_lsm_h%start_index(jp,ip)
    3323                    us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
    3324                    usws_int  = surf_lsm_h%usws(surf_start)
    3325                 ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
    3326                          surf_usm_h%end_index(jp,ip) )  THEN
    3327                    surf_start = surf_usm_h%start_index(jp,ip)
    3328                    us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
    3329                    usws_int  = surf_usm_h%usws(surf_start)
     3304    dt_particle = dt_3d
     3305
     3306
     3307    SELECT CASE ( interpolation_method )
     3308
     3309!
     3310!--    This case uses a simple interpolation method for the particle velocites,
     3311!--    and applying a predictor-corrector method. @attention: for the corrector
     3312!--    step the velocities of t(n+1) are required. However, at this moment of
     3313!--    the time integration they are not free of divergence. This interpolation
     3314!--    method is described in more detail in Grabowski et al., 2018 (GMD).
     3315       CASE ( 'simple_corrector' )
     3316!
     3317!--       Predictor step
     3318          kkw = kp - 1
     3319          DO n = 1, number_of_particles
     3320
     3321             alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3322             u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
     3323
     3324             beta  = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3325             v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
     3326
     3327             gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
     3328                               ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
     3329             w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
     3330
     3331          ENDDO
     3332
     3333!
     3334!--       Corrector step
     3335          DO n = 1, number_of_particles
     3336
     3337             IF ( .NOT. particles(n)%particle_mask ) CYCLE
     3338
     3339             DO nn = 1, iteration_steps
     3340
     3341!
     3342!--             Guess new position
     3343                xp = particles(n)%x + u_int(n) * dt_particle(n)
     3344                yp = particles(n)%y + v_int(n) * dt_particle(n)
     3345                zp = particles(n)%z + w_int(n) * dt_particle(n)
     3346!
     3347!--             x direction
     3348                i_next = FLOOR( xp * ddx , KIND=iwp)
     3349                alpha  = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3350!
     3351!--             y direction
     3352                j_next = FLOOR( yp * ddy )
     3353                beta   = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3354!
     3355!--             z_direction
     3356                k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0)
     3357                gamma = MAX( MIN( ( zp - zw(k_next) ) /                      &
     3358                                  ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp )
     3359!
     3360!--             Calculate part of the corrector step
     3361                unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
     3362                        u_p(k_next+1, j_next,   i_next+1) * alpha
     3363
     3364                vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
     3365                        v_p(k_next+1, j_next+1, i_next  ) * beta
     3366
     3367                wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
     3368                        w_p(k_next+1, j_next, i_next  ) * gamma
     3369
     3370!
     3371!--             Calculate interpolated particle velocity with predictor
     3372!--             corrector step. u_int, v_int and w_int describes the part of
     3373!--             the predictor step. unext, vnext and wnext is the part of the
     3374!--             corrector step. The resulting new position is set below. The
     3375!--             implementation is based on Grabowski et al., 2018 (GMD).
     3376                u_int(n) = 0.5_wp * ( u_int(n) + unext )
     3377                v_int(n) = 0.5_wp * ( v_int(n) + vnext )
     3378                w_int(n) = 0.5_wp * ( w_int(n) + wnext )
     3379
     3380             ENDDO
     3381          ENDDO
     3382
     3383
     3384!
     3385!--    This case uses a simple interpolation method for the particle velocites,
     3386!--    and applying a predictor.
     3387       CASE ( 'simple_predictor' )
     3388!
     3389!--       The particle position for the w velociy is based on the value of kp and kp-1
     3390          kkw = kp - 1
     3391          DO n = 1, number_of_particles
     3392             IF ( .NOT. particles(n)%particle_mask ) CYCLE
     3393
     3394             alpha    = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3395             u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
     3396
     3397             beta     = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3398             v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
     3399
     3400             gamma    = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
     3401                                  ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
     3402             w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
     3403          ENDDO
     3404!
     3405!--    The trilinear interpolation.
     3406       CASE ( 'trilinear' )
     3407
     3408          start_index = grid_particles(kp,jp,ip)%start_index
     3409          end_index   = grid_particles(kp,jp,ip)%end_index
     3410
     3411          DO  nb = 0, 7
     3412!
     3413!--          Interpolate u velocity-component
     3414             i = ip
     3415             j = jp + block_offset(nb)%j_off
     3416             k = kp + block_offset(nb)%k_off
     3417
     3418             DO  n = start_index(nb), end_index(nb)
     3419!
     3420!--             Interpolation of the u velocity component onto particle position.
     3421!--             Particles are interpolation bi-linearly in the horizontal and a
     3422!--             linearly in the vertical. An exception is made for particles below
     3423!--             the first vertical grid level in case of a prandtl layer. In this
     3424!--             case the horizontal particle velocity components are determined using
     3425!--             Monin-Obukhov relations (if branch).
     3426!--             First, check if particle is located below first vertical grid level
     3427!--             above topography (Prandtl-layer height)
     3428!--             Determine vertical index of topography top
     3429                k_wall = get_topography_top_index_ji( jp, ip, 's' )
     3430
     3431                IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     3432!
     3433!--                Resolved-scale horizontal particle velocity is zero below z0.
     3434                   IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
     3435                      u_int(n) = 0.0_wp
     3436                   ELSE
     3437!
     3438!--                   Determine the sublayer. Further used as index.
     3439                      height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
     3440                                           * REAL( number_of_sublayers, KIND=wp )    &
     3441                                           * d_z_p_z0
     3442!
     3443!--                   Calculate LOG(z/z0) for exact particle height. Therefore,
     3444!--                   interpolate linearly between precalculated logarithm.
     3445                      log_z_z0_int = log_z_z0(INT(height_p))                         &
     3446                                       + ( height_p - INT(height_p) )                &
     3447                                       * ( log_z_z0(INT(height_p)+1)                 &
     3448                                            - log_z_z0(INT(height_p))                &
     3449                                         )
     3450!
     3451!--                   Get friction velocity and momentum flux from new surface data
     3452!--                   types.
     3453                      IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3454                           surf_def_h(0)%end_index(jp,ip) )  THEN
     3455                         surf_start = surf_def_h(0)%start_index(jp,ip)
     3456!--                      Limit friction velocity. In narrow canyons or holes the
     3457!--                      friction velocity can become very small, resulting in a too
     3458!--                      large particle speed.
     3459                         us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
     3460                         usws_int  = surf_def_h(0)%usws(surf_start)
     3461                      ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
     3462                               surf_lsm_h%end_index(jp,ip) )  THEN
     3463                         surf_start = surf_lsm_h%start_index(jp,ip)
     3464                         us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
     3465                         usws_int  = surf_lsm_h%usws(surf_start)
     3466                      ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
     3467                               surf_usm_h%end_index(jp,ip) )  THEN
     3468                         surf_start = surf_usm_h%start_index(jp,ip)
     3469                         us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
     3470                         usws_int  = surf_usm_h%usws(surf_start)
     3471                      ENDIF
     3472!
     3473!--                   Neutral solution is applied for all situations, e.g. also for
     3474!--                   unstable and stable situations. Even though this is not exact
     3475!--                   this saves a lot of CPU time since several calls of intrinsic
     3476!--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     3477!--                   as sensitivity studies revealed no significant effect of
     3478!--                   using the neutral solution also for un/stable situations.
     3479                      u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
     3480                                  * log_z_z0_int - u_gtrans
     3481                   ENDIF
     3482!
     3483!--             Particle above the first grid level. Bi-linear interpolation in the
     3484!--             horizontal and linear interpolation in the vertical direction.
     3485                ELSE
     3486                   x  = xv(n) - i * dx
     3487                   y  = yv(n) + ( 0.5_wp - j ) * dy
     3488                   aa = x**2          + y**2
     3489                   bb = ( dx - x )**2 + y**2
     3490                   cc = x**2          + ( dy - y )**2
     3491                   dd = ( dx - x )**2 + ( dy - y )**2
     3492                   gg = aa + bb + cc + dd
     3493
     3494                   u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
     3495                               + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
     3496                               u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     3497
     3498                   IF ( k == nzt )  THEN
     3499                      u_int(n) = u_int_l
     3500                   ELSE
     3501                      u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
     3502                                  + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
     3503                                  u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     3504                      u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *            &
     3505                                 ( u_int_u - u_int_l )
     3506                   ENDIF
    33303507                ENDIF
    3331 
    3332 !
    3333 !--             Neutral solution is applied for all situations, e.g. also for
    3334 !--             unstable and stable situations. Even though this is not exact
    3335 !--             this saves a lot of CPU time since several calls of intrinsic
    3336 !--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    3337 !--             as sensitivity studies revealed no significant effect of
    3338 !--             using the neutral solution also for un/stable situations.
    3339                 u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
    3340                             * log_z_z0_int - u_gtrans
    3341 
    3342              ENDIF
    3343 !
    3344 !--       Particle above the first grid level. Bi-linear interpolation in the
    3345 !--       horizontal and linear interpolation in the vertical direction.
    3346           ELSE
    3347 
    3348              x  = xv(n) - i * dx
    3349              y  = yv(n) + ( 0.5_wp - j ) * dy
    3350              aa = x**2          + y**2
    3351              bb = ( dx - x )**2 + y**2
    3352              cc = x**2          + ( dy - y )**2
    3353              dd = ( dx - x )**2 + ( dy - y )**2
    3354              gg = aa + bb + cc + dd
    3355 
    3356              u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
    3357                          + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
    3358                          u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
    3359 
    3360              IF ( k == nzt )  THEN
    3361                 u_int(n) = u_int_l
    3362              ELSE
    3363                 u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
    3364                             + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
    3365                             u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
    3366                 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
    3367                            ( u_int_u - u_int_l )
    3368              ENDIF
    3369 
    3370           ENDIF
    3371 
    3372        ENDDO
    3373 !
    3374 !--    Same procedure for interpolation of the v velocity-component
    3375        i = ip + block_offset(nb)%i_off
    3376        j = jp
    3377        k = kp + block_offset(nb)%k_off
    3378 
    3379        DO  n = start_index(nb), end_index(nb)
    3380 
    3381 !
    3382 !--       Determine vertical index of topography top
    3383           k_wall = get_topography_top_index_ji( jp,ip, 's' )
    3384 
    3385           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
    3386              IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
    3387 !
    3388 !--             Resolved-scale horizontal particle velocity is zero below z0.
    3389                 v_int(n) = 0.0_wp
    3390              ELSE       
    3391 !
    3392 !--             Determine the sublayer. Further used as index. Please note,
    3393 !--             logarithmus can not be reused from above, as in in case of
    3394 !--             topography particle on u-grid can be above surface-layer height,
    3395 !--             whereas it can be below on v-grid.
    3396                 height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
    3397                                   * REAL( number_of_sublayers, KIND=wp )       &
    3398                                   * d_z_p_z0
    3399 !
    3400 !--             Calculate LOG(z/z0) for exact particle height. Therefore,   
    3401 !--             interpolate linearly between precalculated logarithm.
    3402                 log_z_z0_int = log_z_z0(INT(height_p))                         &
    3403                                  + ( height_p - INT(height_p) )                &
    3404                                  * ( log_z_z0(INT(height_p)+1)                 &
    3405                                       - log_z_z0(INT(height_p))                &
    3406                                    )
    3407 !
    3408 !--             Get friction velocity and momentum flux from new surface data
    3409 !--             types.
    3410                 IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
    3411                      surf_def_h(0)%end_index(jp,ip) )  THEN
    3412                    surf_start = surf_def_h(0)%start_index(jp,ip)
    3413 !--                Limit friction velocity. In narrow canyons or holes the
    3414 !--                friction velocity can become very small, resulting in a too
    3415 !--                large particle speed.
    3416                    us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp ) 
    3417                    vsws_int  = surf_def_h(0)%vsws(surf_start)
    3418                 ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
    3419                          surf_lsm_h%end_index(jp,ip) )  THEN
    3420                    surf_start = surf_lsm_h%start_index(jp,ip)
    3421                    us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp ) 
    3422                    vsws_int  = surf_lsm_h%vsws(surf_start)
    3423                 ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
    3424                          surf_usm_h%end_index(jp,ip) )  THEN
    3425                    surf_start = surf_usm_h%start_index(jp,ip)
    3426                    us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp ) 
    3427                    vsws_int  = surf_usm_h%vsws(surf_start)
    3428                 ENDIF
    3429 !
    3430 !--             Neutral solution is applied for all situations, e.g. also for
    3431 !--             unstable and stable situations. Even though this is not exact
    3432 !--             this saves a lot of CPU time since several calls of intrinsic
    3433 !--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    3434 !--             as sensitivity studies revealed no significant effect of
    3435 !--             using the neutral solution also for un/stable situations.
    3436                 v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
    3437                          * log_z_z0_int - v_gtrans
    3438 
    3439              ENDIF
    3440 
    3441           ELSE
    3442              x  = xv(n) + ( 0.5_wp - i ) * dx
    3443              y  = yv(n) - j * dy
    3444              aa = x**2          + y**2
    3445              bb = ( dx - x )**2 + y**2
    3446              cc = x**2          + ( dy - y )**2
    3447              dd = ( dx - x )**2 + ( dy - y )**2
    3448              gg = aa + bb + cc + dd
    3449 
    3450              v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
    3451                        + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
    3452                        ) / ( 3.0_wp * gg ) - v_gtrans
    3453 
    3454              IF ( k == nzt )  THEN
    3455                 v_int(n) = v_int_l
    3456              ELSE
    3457                 v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
    3458                           + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
    3459                           ) / ( 3.0_wp * gg ) - v_gtrans
    3460                 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
    3461                                   ( v_int_u - v_int_l )
    3462              ENDIF
    3463 
    3464           ENDIF
    3465 
    3466        ENDDO
    3467 !
    3468 !--    Same procedure for interpolation of the w velocity-component
    3469        i = ip + block_offset(nb)%i_off
    3470        j = jp + block_offset(nb)%j_off
    3471        k = kp - 1
    3472 
    3473        DO  n = start_index(nb), end_index(nb)
    3474 
    3475           IF ( vertical_particle_advection(particles(n)%group) )  THEN
    3476 
    3477              x  = xv(n) + ( 0.5_wp - i ) * dx
    3478              y  = yv(n) + ( 0.5_wp - j ) * dy
    3479              aa = x**2          + y**2
    3480              bb = ( dx - x )**2 + y**2
    3481              cc = x**2          + ( dy - y )**2
    3482              dd = ( dx - x )**2 + ( dy - y )**2
    3483              gg = aa + bb + cc + dd
    3484 
    3485              w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
    3486                        + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
    3487                        ) / ( 3.0_wp * gg )
    3488 
    3489              IF ( k == nzt )  THEN
    3490                 w_int(n) = w_int_l
    3491              ELSE
    3492                 w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
    3493                             ( gg-bb ) * w(k+1,j,i+1) + &
    3494                             ( gg-cc ) * w(k+1,j+1,i) + &
    3495                             ( gg-dd ) * w(k+1,j+1,i+1) &
    3496                           ) / ( 3.0_wp * gg )
    3497                 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
    3498                            ( w_int_u - w_int_l )
    3499              ENDIF
    3500 
    3501           ELSE
    3502 
    3503              w_int(n) = 0.0_wp
    3504 
    3505           ENDIF
    3506 
    3507        ENDDO
    3508 
    3509     ENDDO
     3508             ENDDO
     3509!
     3510!--          Same procedure for interpolation of the v velocity-component
     3511             i = ip + block_offset(nb)%i_off
     3512             j = jp
     3513             k = kp + block_offset(nb)%k_off
     3514
     3515             DO  n = start_index(nb), end_index(nb)
     3516!
     3517!--             Determine vertical index of topography top
     3518                k_wall = get_topography_top_index_ji( jp,ip, 's' )
     3519
     3520                IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     3521                   IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
     3522!
     3523!--                   Resolved-scale horizontal particle velocity is zero below z0.
     3524                      v_int(n) = 0.0_wp
     3525                   ELSE
     3526!
     3527!--                   Determine the sublayer. Further used as index. Please note,
     3528!--                   logarithmus can not be reused from above, as in in case of
     3529!--                   topography particle on u-grid can be above surface-layer height,
     3530!--                   whereas it can be below on v-grid.
     3531                      height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
     3532                                        * REAL( number_of_sublayers, KIND=wp )       &
     3533                                        * d_z_p_z0
     3534!
     3535!--                   Calculate LOG(z/z0) for exact particle height. Therefore,
     3536!--                   interpolate linearly between precalculated logarithm.
     3537                      log_z_z0_int = log_z_z0(INT(height_p))                         &
     3538                                       + ( height_p - INT(height_p) )                &
     3539                                       * ( log_z_z0(INT(height_p)+1)                 &
     3540                                            - log_z_z0(INT(height_p))                &
     3541                                         )
     3542!
     3543!--                   Get friction velocity and momentum flux from new surface data
     3544!--                   types.
     3545                      IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3546                           surf_def_h(0)%end_index(jp,ip) )  THEN
     3547                         surf_start = surf_def_h(0)%start_index(jp,ip)
     3548!--                      Limit friction velocity. In narrow canyons or holes the
     3549!--                      friction velocity can become very small, resulting in a too
     3550!--                      large particle speed.
     3551                         us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
     3552                         vsws_int  = surf_def_h(0)%vsws(surf_start)
     3553                      ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
     3554                               surf_lsm_h%end_index(jp,ip) )  THEN
     3555                         surf_start = surf_lsm_h%start_index(jp,ip)
     3556                         us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
     3557                         vsws_int  = surf_lsm_h%vsws(surf_start)
     3558                      ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
     3559                               surf_usm_h%end_index(jp,ip) )  THEN
     3560                         surf_start = surf_usm_h%start_index(jp,ip)
     3561                         us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
     3562                         vsws_int  = surf_usm_h%vsws(surf_start)
     3563                      ENDIF
     3564!
     3565!--                   Neutral solution is applied for all situations, e.g. also for
     3566!--                   unstable and stable situations. Even though this is not exact
     3567!--                   this saves a lot of CPU time since several calls of intrinsic
     3568!--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     3569!--                   as sensitivity studies revealed no significant effect of
     3570!--                   using the neutral solution also for un/stable situations.
     3571                      v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
     3572                               * log_z_z0_int - v_gtrans
     3573
     3574                   ENDIF
     3575                ELSE
     3576                   x  = xv(n) + ( 0.5_wp - i ) * dx
     3577                   y  = yv(n) - j * dy
     3578                   aa = x**2          + y**2
     3579                   bb = ( dx - x )**2 + y**2
     3580                   cc = x**2          + ( dy - y )**2
     3581                   dd = ( dx - x )**2 + ( dy - y )**2
     3582                   gg = aa + bb + cc + dd
     3583
     3584                   v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
     3585                             + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
     3586                             ) / ( 3.0_wp * gg ) - v_gtrans
     3587
     3588                   IF ( k == nzt )  THEN
     3589                      v_int(n) = v_int_l
     3590                   ELSE
     3591                      v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
     3592                                + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
     3593                                ) / ( 3.0_wp * gg ) - v_gtrans
     3594                      v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
     3595                                        ( v_int_u - v_int_l )
     3596                   ENDIF
     3597                ENDIF
     3598             ENDDO
     3599!
     3600!--          Same procedure for interpolation of the w velocity-component
     3601             i = ip + block_offset(nb)%i_off
     3602             j = jp + block_offset(nb)%j_off
     3603             k = kp - 1
     3604
     3605             DO  n = start_index(nb), end_index(nb)
     3606                IF ( vertical_particle_advection(particles(n)%group) )  THEN
     3607                   x  = xv(n) + ( 0.5_wp - i ) * dx
     3608                   y  = yv(n) + ( 0.5_wp - j ) * dy
     3609                   aa = x**2          + y**2
     3610                   bb = ( dx - x )**2 + y**2
     3611                   cc = x**2          + ( dy - y )**2
     3612                   dd = ( dx - x )**2 + ( dy - y )**2
     3613                   gg = aa + bb + cc + dd
     3614
     3615                   w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
     3616                             + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
     3617                             ) / ( 3.0_wp * gg )
     3618
     3619                   IF ( k == nzt )  THEN
     3620                      w_int(n) = w_int_l
     3621                   ELSE
     3622                      w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
     3623                                  ( gg-bb ) * w(k+1,j,i+1) + &
     3624                                  ( gg-cc ) * w(k+1,j+1,i) + &
     3625                                  ( gg-dd ) * w(k+1,j+1,i+1) &
     3626                                ) / ( 3.0_wp * gg )
     3627                      w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
     3628                                 ( w_int_u - w_int_l )
     3629                   ENDIF
     3630                ELSE
     3631                   w_int(n) = 0.0_wp
     3632                ENDIF
     3633             ENDDO
     3634          ENDDO
     3635
     3636       CASE DEFAULT
     3637          WRITE( message_string, * )  'unknown particle velocity interpolation method = "',  &
     3638                                       TRIM( interpolation_method ), '"'
     3639          CALL message( 'lpm_advec', 'PA0660', 1, 2, 0, 6, 0 )
     3640
     3641    END SELECT
    35103642
    35113643!-- Interpolate and calculate quantities needed for calculating the SGS
     
    38193951!--    and calculate SGS velocities again
    38203952       dz_temp = zw(kp)-zw(kp-1)
    3821        
     3953
    38223954       DO  nb = 0, 7
    38233955          DO  n = start_index(nb), end_index(nb)
     
    38263958                 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n))  .OR.   &
    38273959                 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN
    3828                
     3960
    38293961                dt_particle(n) = 0.9_wp * MIN(                                 &
    38303962                                 ( dx / ABS( u_int(n) + rvar1_temp(n) ) ),     &
     
    38373969                rvar2_temp(n) = particles(n)%rvar2
    38383970                rvar3_temp(n) = particles(n)%rvar3
    3839                
     3971
    38403972                de_dt_min = - e_int(n) / dt_particle(n)
    38413973
     
    38573989                                        de_dz_int(n), de_dt, diss_int(n),       &
    38583990                                        dt_particle(n), rg(n,3), term_1_2(n) )
    3859              ENDIF                           
    3860              
     3991             ENDIF
     3992
    38613993!
    38623994!--          Update particle velocites
     
    38734005          ENDDO
    38744006       ENDDO
    3875        
     4007
    38764008    ELSE
    38774009!
     
    38834015
    38844016    dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
    3885 
    38864017    IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
    3887        DO  nb = 0, 7
    3888           DO  n = start_index(nb), end_index(nb)
     4018!
     4019!--    Decide whether the particle loop runs over the subboxes or only over 1,
     4020!--    number_of_particles. This depends on the selected interpolation method.
     4021!--    If particle interpolation method is not trilinear, then the sorting within
     4022!--    subboxes is not required. However, therefore the index start_index(nb) and
     4023!--    end_index(nb) are not defined and the loops are still over
     4024!--    number_of_particles. @todo find a more generic way to write this loop or
     4025!--    delete trilinear interpolation
     4026       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4027          subbox_start = 0
     4028          subbox_end   = 7
     4029       ELSE
     4030          subbox_start = 1
     4031          subbox_end   = 1
     4032       ENDIF
     4033!
     4034!--    loop over subboxes. In case of simple interpolation scheme no subboxes
     4035!--    are introduced, as they are not required. Accordingly, this loops goes
     4036!--    from 1 to 1.
     4037       DO  nb = subbox_start, subbox_end
     4038          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4039             particle_start = start_index(nb)
     4040             particle_end   = end_index(nb)
     4041          ELSE
     4042             particle_start = 1
     4043             particle_end   = number_of_particles
     4044          ENDIF
     4045!
     4046!--         Loop from particle start to particle end
     4047            DO  n = particle_start, particle_end
    38894048
    38904049!
     
    39154074                IF ( cloud_droplets )  THEN
    39164075!
    3917 !--                Terminal velocity is computed for vertical direction (Rogers et 
     4076!--                Terminal velocity is computed for vertical direction (Rogers et
    39184077!--                al., 1993, J. Appl. Meteorol.)
    39194078                   diameter = particles(n)%radius * 2000.0_wp !diameter in mm
     
    39744133          ENDDO
    39754134       ENDDO
    3976    
     4135
    39774136    ELSE
    3978 
    3979        DO  nb = 0, 7
    3980           DO  n = start_index(nb), end_index(nb)
     4137!
     4138!--    Decide whether the particle loop runs over the subboxes or only over 1,
     4139!--    number_of_particles. This depends on the selected interpolation method.
     4140       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4141          subbox_start = 0
     4142          subbox_end   = 7
     4143       ELSE
     4144          subbox_start = 1
     4145          subbox_end   = 1
     4146       ENDIF
     4147!--    loop over subboxes. In case of simple interpolation scheme no subboxes
     4148!--    are introduced, as they are not required. Accordingly, this loops goes
     4149!--    from 1 to 1.
     4150       DO  nb = subbox_start, subbox_end
     4151          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4152             particle_start = start_index(nb)
     4153             particle_end   = end_index(nb)
     4154          ELSE
     4155             particle_start = 1
     4156             particle_end   = number_of_particles
     4157          ENDIF
     4158!
     4159!--         Loop from particle start to particle end
     4160            DO  n = particle_start, particle_end
     4161
    39814162!
    39824163!--          Transport of particles with inertia
     
    39884169             IF ( cloud_droplets )  THEN
    39894170!
    3990 !--             Terminal velocity is computed for vertical direction (Rogers et al., 
     4171!--             Terminal velocity is computed for vertical direction (Rogers et al.,
    39914172!--             1993, J. Appl. Meteorol.)
    39924173                diameter = particles(n)%radius * 2000.0_wp !diameter in mm
     
    40544235    particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
    40554236
    4056     DO  nb = 0, 7
    4057        DO  n = start_index(nb), end_index(nb)
     4237!
     4238!--    loop over subboxes. In case of simple interpolation scheme no subboxes
     4239!--    are introduced, as they are not required. Accordingly, this loops goes
     4240!--    from 1 to 1.
     4241!
     4242!-- Decide whether the particle loop runs over the subboxes or only over 1,
     4243!-- number_of_particles. This depends on the selected interpolation method.
     4244    IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4245       subbox_start = 0
     4246       subbox_end   = 7
     4247    ELSE
     4248       subbox_start = 1
     4249       subbox_end   = 1
     4250    ENDIF
     4251    DO  nb = subbox_start, subbox_end
     4252       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4253          particle_start = start_index(nb)
     4254          particle_end   = end_index(nb)
     4255       ELSE
     4256          particle_start = 1
     4257          particle_end   = number_of_particles
     4258       ENDIF
     4259!
     4260!--    Loop from particle start to particle end
     4261       DO  n = particle_start, particle_end
    40584262!
    40594263!--       Increment the particle age and the total time that the particle
     
    40634267
    40644268!
    4065 !--       Check whether there is still a particle that has not yet completed 
     4269!--       Check whether there is still a particle that has not yet completed
    40664270!--       the total LES timestep
    40674271          IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
     
    40864290 SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
    40874291                                dt_n, rg_n, fac )
    4088                                
     4292
    40894293    REAL(wp) ::  a1      !< dummy argument
    40904294    REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position
     
    41614365    INTEGER(iwp), INTENT(IN) ::  j !< grid index of particle box along y
    41624366    INTEGER(iwp), INTENT(IN) ::  k !< grid index of particle box along z
    4163    
     4367
    41644368    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
    41654369    INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
     
    41844388    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
    41854389    INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
    4186    
     4390
    41874391    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
    41884392    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
     
    42254429       CASE ( 'bottom/top' )
    42264430
    4227 !     IF ( location_bc == 'bottom/top' )  THEN
    4228 
    42294431!
    42304432!--    Apply boundary conditions to those particles that have crossed the top or
     
    42704472             ENDIF
    42714473          ENDIF
    4272          
     4474
    42734475          IF ( particles(n)%z < zw(0)  .AND.  particles(n)%particle_mask )  THEN
    42744476             IF ( ibc_par_b == 1 )  THEN
     
    42904492       ENDDO
    42914493
    4292 !     ELSEIF ( location_bc == 'walls' )  THEN
    42934494      CASE ( 'walls' )
    42944495
     
    43474548          ELSE
    43484549             zwall = zw(k-1)
    4349           ENDIF     
     4550          ENDIF
    43504551!
    43514552!--       Initialize flags to check if particle reflection is necessary
     
    44184619                              MIN( prt_z - pos_z_old, -1E-30_wp ),             &
    44194620                              prt_z > pos_z_old )
    4420                      
     4621
    44214622          x_ind(t_index)   = i1
    44224623          y_ind(t_index)   = j1
     
    44294630             cross_wall_z = .TRUE.
    44304631          ENDIF
    4431          
     4632
    44324633          t_index_number = t_index - 1
    44334634!
     
    44844685             t_old = 0.0_wp
    44854686             DO t_index = 1, t_index_number
    4486 !           
     4687!
    44874688!--             Calculate intermediate particle position according to the
    44884689!--             timesteps a particle reaches any wall.
     
    45834784                    y_ind(t_index:t_index_number) = j2
    45844785                ENDIF !particle reflection in y direction done
    4585                
     4786
    45864787!
    45874788!--             Check if a particle needs to be reflected at any xy-wall. If
     
    46234824                    z_ind(t_index:t_index_number) = k2
    46244825                ENDIF !particle reflection in z direction done               
    4625                
     4826
    46264827!
    46274828!--             Swap time
     
    46554856
    46564857 END SUBROUTINE lpm_boundary_conds
    4657  
    4658  
     4858
     4859
     4860!------------------------------------------------------------------------------!
     4861! Description:
     4862! ------------
     4863!> Calculates change in droplet radius by condensation/evaporation, using
     4864!> either an analytic formula or by numerically integrating the radius growth
     4865!> equation including curvature and solution effects using Rosenbrocks method
     4866!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
     4867!> The analytical formula and growth equation follow those given in
     4868!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
     4869!------------------------------------------------------------------------------!
    46594870 SUBROUTINE lpm_droplet_condensation (i,j,k)
    46604871
     
    48125023             DO
    48135024
    4814                 drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
     5025                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
    48155026                                                          afactor / r_ros +    &
    48165027                                                          bfactor / r_ros**3   &
     
    48185029
    48195030                d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
    4820                                                 (e_a / e_s - 1.0) * r_ros**4 - &
    4821                                                 afactor * r0 * r_ros**2 -      &
    4822                                                 2.0 * afactor * r_ros**3 +     &
    4823                                                 3.0 * bfactor * r0 +           &
    4824                                                 4.0 * bfactor * r_ros          &
     5031                                            (e_a / e_s - 1.0_wp ) * r_ros**4 - &
     5032                                            afactor * r0 * r_ros**2 -          &
     5033                                            2.0_wp * afactor * r_ros**3 +      &
     5034                                            3.0_wp * bfactor * r0 +            &
     5035                                            4.0_wp * bfactor * r_ros           &
    48255036                                                            )                  &
    48265037                          / ( r_ros**4 * ( r_ros + r0 )**2 )
    48275038
    4828                 k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
     5039                k1    = drdt / ( 1.0_wp - gamma * dt_ros * d2rdtdr )
    48295040
    48305041                r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
    48315042                r_err = r_ros
    48325043
    4833                 drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -  &
    4834                                                            afactor / r_ros +   &
    4835                                                            bfactor / r_ros**3  &
     5044                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
     5045                                                           afactor / r_ros +    &
     5046                                                           bfactor / r_ros**3   &
    48365047                                                         ) / ( r_ros + r0 )
    48375048
    48385049                k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
    4839                      ( 1.0 - dt_ros * gamma * d2rdtdr )
    4840 
    4841                 r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
     5050                     ( 1.0_wp - dt_ros * gamma * d2rdtdr )
     5051
     5052                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5_wp * k1 + 0.5_wp * k2), particles(n)%aux1)
    48425053   !
    48435054   !--          Check error of the solution, and reduce dt_ros if necessary.
     
    49095120          particles(n)%class = MAX( particles(n)%class, 1 )
    49105121       ENDIF
    4911  !
    4912  !--   Store new radius to particle features
     5122!
     5123!--    Store new radius to particle features
    49135124       particles(n)%radius = new_r(n)
    49145125
     
    51345345          aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
    51355346                                             particles(1:number_of_particles)%aux1**3       * &
    5136                                              4.0 / 3.0 * pi * rho_s
     5347                                             4.0_wp / 3.0_wp * pi * rho_s
    51375348       ENDIF
    51385349!
     
    51445355!--          For collisions, the weighting factor of at least one super-droplet
    51455356!--          needs to be larger or equal to one.
    5146              IF ( MIN( weight(n), weight(m) ) .LT. 1.0 )  CYCLE
     5357             IF ( MIN( weight(n), weight(m) ) .LT. 1.0_wp )  CYCLE
    51475358!
    51485359!--          Get mass of individual droplets (aerosols)
     
    51755386             ENDIF
    51765387
    5177              IF ( collection_probability .GT. 0.0 )  THEN
     5388             IF ( collection_probability .GT. 0.0_wp )  THEN
    51785389!
    51795390!--             Super-droplet n collects droplets of super-droplet m
     
    53635574                                       ( hkernel(i,j), i = 1,radius_classes )
    53645575          ENDDO
    5365          
     5576
    53665577          DO  k = 1, dissipation_classes
    53675578             DO  i = 1, radius_classes
     
    61546365                     ql(k,j,i) < ql_crit )  CYCLE
    61556366                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    6156 !               
     6367!
    61576368!--             Start splitting operations. Each particle is checked if it
    61586369!--             fulfilled the splitting criterion's. In splitting mode 'const'   
     
    61686379                        particles(n)%radius >= radius_split  .AND.             &
    61696380                        particles(n)%weight_factor >= weight_factor_split )    &
    6170                    THEN         
     6381                   THEN
    61716382!
    61726383!--                   Calculate the new number of particles.
    6173                       new_size = prt_count(k,j,i) + splitting_factor - 1                                   
     6384                      new_size = prt_count(k,j,i) + splitting_factor - 1
    61746385!
    61756386!--                   Cycle if maximum number of particles per grid box
    61766387!--                   is greater than the allowed maximum number.
    6177                       IF ( new_size >= max_number_particles_per_gridbox )  CYCLE                     
     6388                      IF ( new_size >= max_number_particles_per_gridbox )  CYCLE
    61786389!
    61796390!--                   Reallocate particle array if necessary.
     
    61916402                      DO  jpp = 1, splitting_factor-1
    61926403                         grid_particles(k,j,i)%particles(jpp+old_size) =       &
    6193                             tmp_particle         
     6404                            tmp_particle
    61946405                      ENDDO 
    61956406                      new_particles_gb = new_particles_gb + splitting_factor - 1
     
    61976408!--                   Save the new number of super droplets for every grid box.
    61986409                      prt_count(k,j,i) = prt_count(k,j,i) +                    &
    6199                                          splitting_factor - 1         
     6410                                         splitting_factor - 1
    62006411                   ENDIF
    62016412                ENDDO
     
    62166427       m3           = 0.0_wp
    62176428       m3_total     = 0.0_wp
    6218        nr           = 0.0_wp   
     6429       nr           = 0.0_wp
    62196430       nrclgb       = 0.0_wp
    6220        nrclgb_total = 0.0_wp 
     6431       nrclgb_total = 0.0_wp
    62216432       nr_total     = 0.0_wp
    62226433       rm           = 0.0_wp
    62236434       rm_total     = 0.0_wp
    6224        
     6435
    62256436       DO  i = nxl, nxr
    62266437          DO  j = nys, nyn
     
    62316442                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    62326443                nrclgb = nrclgb + 1.0_wp
    6233 !               
    6234 !--             Calculate moments of DSD.               
     6444!
     6445!--             Calculate moments of DSD.
    62356446                DO  n = 1, number_of_particles
    62366447                   IF ( particles(n)%particle_mask  .AND.                      &
     
    62426453                                 particles(n)%weight_factor
    62436454                      IF ( isf == 1 )  THEN           
    6244                          diameter   = particles(n)%radius * 2.0_wp           
     6455                         diameter   = particles(n)%radius * 2.0_wp
    62456456                         lwc = lwc + factor_volume_to_mass *                   &
    62466457                                     particles(n)%radius**3 *                  &
    62476458                                     particles(n)%weight_factor
    6248                          m1  = m1  + particles(n)%weight_factor * diameter                                               
    6249                          m2  = m2  + particles(n)%weight_factor * diameter**2           
     6459                         m1  = m1  + particles(n)%weight_factor * diameter
     6460                         m2  = m2  + particles(n)%weight_factor * diameter**2
    62506461                         m3  = m3  + particles(n)%weight_factor * diameter**3
    6251                       ENDIF   
    6252                    ENDIF                       
     6462                      ENDIF
     6463                   ENDIF
    62536464                ENDDO
    62546465             ENDDO
     
    62876498                            ( nr_total * factor_volume_to_mass )               &
    62886499                          )**0.3333333_wp, 0.0_wp, nrclgb_total > 0.0_wp       &
    6289                        )                         
     6500                       )
    62906501!
    62916502!--    Check which function should be used to approximate the DSD.
     
    63036514                             0.0_wp, nrclgb_total > 0.0_wp                     &
    63046515                           )
    6305           zeta = m1_total * m3_total / m2_total**2                             
     6516          zeta = m1_total * m3_total / m2_total**2
    63066517          mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /                &
    63076518                        ( zeta - 1.0_wp ), 0.0_wp                              &
     
    63096520
    63106521          lambda = ( pirho_l * nr_total / lwc_total *                          &
    6311                      ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &                                         
     6522                     ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &
    63126523                   )**0.3333333_wp
    63136524          nr0 = nr_total / gamma( mu + 1.0_wp ) * lambda**( mu + 1.0_wp )
    6314          
     6525
    63156526          DO  n = 0, n_max-1
    6316              diameter  = r_bin_mid(n) * 2.0_wp           
     6527             diameter  = r_bin_mid(n) * 2.0_wp
    63176528             an_spl(n) = nr0 * diameter**mu * EXP( -lambda * diameter ) *      &
    63186529                         ( r_bin(n+1) - r_bin(n) ) * 2.0_wp
     
    63386549!--    Criterion to avoid super droplets with a weighting factor < 1.0.
    63396550       an_spl = MAX(an_spl, 1.0_wp)
    6340  
     6551
    63416552       DO  i = nxl, nxr
    63426553          DO  j = nys, nyn
    63436554             DO  k = nzb+1, nzt
    63446555                number_of_particles = prt_count(k,j,i)
    6345                 IF ( number_of_particles <= 0  .OR.                            & 
     6556                IF ( number_of_particles <= 0  .OR.                            &
    63466557                     ql(k,j,i) < ql_crit )  CYCLE
    63476558                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    6348                 new_particles_gb = 0         
    6349 !               
     6559                new_particles_gb = 0
     6560!
    63506561!--             Start splitting operations. Each particle is checked if it
    6351 !--             fulfilled the splitting criterion's. In splitting mode 'cl_av'   
    6352 !--             a critical radius (radius_split) and a splitting function must 
     6562!--             fulfilled the splitting criterion's. In splitting mode 'cl_av'
     6563!--             a critical radius (radius_split) and a splitting function must
    63536564!--             be prescribed (see particles_par). The critical weighting factor
    63546565!--             is calculated while approximating a 'gamma', 'log' or 'exp'-
    6355 !--             drop size distribution. In this mode the DSD is calculated as 
    6356 !--             an average over all cloudy grid boxes. Super droplets which 
     6566!--             drop size distribution. In this mode the DSD is calculated as
     6567!--             an average over all cloudy grid boxes. Super droplets which
    63576568!--             have a larger radius and larger weighting factor are split into
    6358 !--             'splitting_factor' super droplets. In this case the splitting 
    6359 !--             factor is calculated of weighting factor of the super droplet 
    6360 !--             and the approximated number concentration for droplet of such 
    6361 !--             a size. Due to the splitting, the weighting factor of the 
    6362 !--             super droplet and all created clones is reduced by the factor 
     6569!--             'splitting_factor' super droplets. In this case the splitting
     6570!--             factor is calculated of weighting factor of the super droplet
     6571!--             and the approximated number concentration for droplet of such
     6572!--             a size. Due to the splitting, the weighting factor of the
     6573!--             super droplet and all created clones is reduced by the factor
    63636574!--             of 'splitting_facor'.
    63646575                DO  n = 1, number_of_particles
     
    63796590                         IF ( splitting_factor < 2 )  CYCLE
    63806591!
    6381 !--                      Calculate the new number of particles.                                                           
     6592!--                      Calculate the new number of particles.
    63826593                         new_size = prt_count(k,j,i) + splitting_factor - 1
    63836594!
    63846595!--                      Cycle if maximum number of particles per grid box
    6385 !--                      is greater than the allowed maximum number.                         
     6596!--                      is greater than the allowed maximum number.
    63866597                         IF ( new_size >= max_number_particles_per_gridbox )   &
    63876598                         CYCLE
     
    63916602                            CALL realloc_particles_array(i,j,k,new_size)
    63926603                         ENDIF
    6393                          old_size  = prt_count(k,j,i)                             
     6604                         old_size  = prt_count(k,j,i)
    63946605                         new_particles_gb = new_particles_gb +                 &
    63956606                                            splitting_factor - 1
     
    64006611                         tmp_particle = particles(n)
    64016612!
    6402 !--                      Create splitting_factor-1 new particles.                                                 
     6613!--                      Create splitting_factor-1 new particles.
    64036614                         DO  jpp = 1, splitting_factor-1
    64046615                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
    6405                                                                     tmp_particle         
     6616                                                                    tmp_particle
    64066617                         ENDDO
    6407 !   
     6618!
    64086619!--                      Save the new number of super droplets.
    64096620                         prt_count(k,j,i) = prt_count(k,j,i) +                 &
     
    64126623                   ENDDO
    64136624                ENDDO
    6414                
     6625
    64156626             ENDDO
    64166627          ENDDO
     
    64226633          DO  j = nys, nyn
    64236634             DO  k = nzb+1, nzt
    6424              
    6425 !
    6426 !--             Initialize summing variables.             
     6635
     6636!
     6637!--             Initialize summing variables.
    64276638                lwc = 0.0_wp
    64286639                m1  = 0.0_wp
     
    64316642                nr  = 0.0_wp
    64326643                rm  = 0.0_wp 
    6433                
     6644
    64346645                new_particles_gb = 0
    64356646                number_of_particles = prt_count(k,j,i)
     
    64376648                     ql(k,j,i) < ql_crit )  CYCLE
    64386649                particles => grid_particles(k,j,i)%particles
    6439 !               
    6440 !--             Calculate moments of DSD.               
     6650!
     6651!--             Calculate moments of DSD.
    64416652                DO  n = 1, number_of_particles
    64426653                   IF ( particles(n)%particle_mask  .AND.                      &
     
    64476658                                 particles(n)%radius**3  *                     &
    64486659                                 particles(n)%weight_factor
    6449                       IF ( isf == 1 )  THEN           
    6450                          diameter   = particles(n)%radius * 2.0_wp           
     6660                      IF ( isf == 1 )  THEN
     6661                         diameter   = particles(n)%radius * 2.0_wp
    64516662                         lwc = lwc + factor_volume_to_mass *                   &
    6452                                      particles(n)%radius**3 *                  & 
     6663                                     particles(n)%radius**3 *                  &
    64536664                                     particles(n)%weight_factor
    64546665                         m1  = m1 + particles(n)%weight_factor * diameter
    64556666                         m2  = m2 + particles(n)%weight_factor * diameter**2
    64566667                         m3  = m3 + particles(n)%weight_factor * diameter**3
    6457                       ENDIF     
    6458                    ENDIF                                           
    6459                 ENDDO 
    6460 
    6461                 IF ( nr <= 0.0  .OR.  rm <= 0.0_wp )  CYCLE
    6462 !
    6463 !--             Calculate mean volume averaged radius.               
     6668                      ENDIF
     6669                   ENDIF
     6670                ENDDO
     6671
     6672                IF ( nr <= 0.0_wp  .OR.  rm <= 0.0_wp )  CYCLE
     6673!
     6674!--             Calculate mean volume averaged radius.
    64646675                rm = ( rm / ( nr * factor_volume_to_mass ) )**0.3333333_wp
    64656676!
    6466 !--             Check which function should be used to approximate the DSD.             
     6677!--             Check which function should be used to approximate the DSD.
    64676678                IF ( isf == 1 )  THEN
    64686679!
     
    64756686                   lambda = ( pirho_l * nr / lwc *                             &
    64766687                              ( mu + 3.0_wp ) * ( mu + 2.0_wp ) *              &
    6477                               ( mu + 1.0_wp )                                  &                                 
     6688                              ( mu + 1.0_wp )                                  &
    64786689                            )**0.3333333_wp
    64796690                   nr0 =  ( nr / (gamma( mu + 1.0_wp ) ) ) *                   &
     
    64816692
    64826693                   DO  n = 0, n_max-1
    6483                       diameter         = r_bin_mid(n) * 2.0_wp           
     6694                      diameter         = r_bin_mid(n) * 2.0_wp
    64846695                      an_spl(n) = nr0 * diameter**mu *                         &
    64856696                                  EXP( -lambda * diameter ) *                  &
     
    64896700!
    64906701!--                Lognormal size distribution to calculate critical
    6491 !--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).   
     6702!--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).
    64926703                   DO  n = 0, n_max-1
    64936704                      an_spl(n) = nr / ( SQRT( 2.0_wp * pi ) *                 &
     
    64956706                                        ) *                                    &
    64966707                                  EXP( -( LOG( r_bin_mid(n) / rm )**2 ) /      &
    6497                                         ( 2.0_wp * LOG(sigma_log)**2 )         & 
    6498                                       ) *                                      & 
     6708                                        ( 2.0_wp * LOG(sigma_log)**2 )         &
     6709                                      ) *                                      &
    64996710                                  ( r_bin(n+1) - r_bin(n) )
    65006711                   ENDDO
     
    65096720                   ENDDO
    65106721                ENDIF
    6511                
    6512 !
    6513 !--             Criterion to avoid super droplets with a weighting factor < 1.0.                                   
     6722
     6723!
     6724!--             Criterion to avoid super droplets with a weighting factor < 1.0.
    65146725                an_spl = MAX(an_spl, 1.0_wp)
    6515 !               
     6726!
    65166727!--             Start splitting operations. Each particle is checked if it
    6517 !--             fulfilled the splitting criterion's. In splitting mode 'gb_av'   
     6728!--             fulfilled the splitting criterion's. In splitting mode 'gb_av'
    65186729!--             a critical radius (radius_split) and a splitting function must
    65196730!--             be prescribed (see particles_par). The critical weighting factor
     
    65466757
    65476758!
    6548 !--                      Calculate the new number of particles.                                                                                         
     6759!--                      Calculate the new number of particles.
    65496760                         new_size = prt_count(k,j,i) + splitting_factor - 1
    65506761!
    65516762!--                      Cycle if maximum number of particles per grid box
    6552 !--                      is greater than the allowed maximum number.                                                 
     6763!--                      is greater than the allowed maximum number.
    65536764                         IF ( new_size >= max_number_particles_per_gridbox )   &
    65546765                         CYCLE
    65556766!
    6556 !--                      Reallocate particle array if necessary.                         
     6767!--                      Reallocate particle array if necessary.
    65576768                         IF ( new_size > SIZE(particles) )  THEN
    65586769                            CALL realloc_particles_array(i,j,k,new_size)
     
    65756786                                               splitting_factor - 1
    65766787                         new_particles_gb    = new_particles_gb +              &
    6577                                                splitting_factor - 1                                       
     6788                                               splitting_factor - 1
    65786789                      ENDIF
    65796790                   ENDDO 
     
    66026813 SUBROUTINE lpm_merging
    66036814
    6604     INTEGER(iwp) ::  i         !<     
    6605     INTEGER(iwp) ::  j         !<       
    6606     INTEGER(iwp) ::  k         !<       
     6815    INTEGER(iwp) ::  i         !<
     6816    INTEGER(iwp) ::  j         !<
     6817    INTEGER(iwp) ::  k         !<
    66076818    INTEGER(iwp) ::  n         !<
    66086819    INTEGER(iwp) ::  merge_drp = 0     !< number of merged droplets
    6609    
    6610    
     6820
     6821
    66116822    REAL(wp) ::  ql_crit = 1.0E-5_wp  !< threshold lwc for cloudy grid cells
    66126823                                      !< (e.g. Siebesma et al 2003, JAS, 60)
    6613    
     6824
    66146825    CALL cpu_log( log_point_s(81), 'lpm_merging', 'start' )
    66156826
    66166827    merge_drp  = 0
    6617    
     6828
    66186829    IF ( weight_factor_merge == -1.0_wp )  THEN
    66196830       weight_factor_merge = 0.5_wp * initial_weighting_factor
     
    66236834       DO  j = nys, nyn
    66246835          DO  k = nzb+1, nzt
    6625    
     6836
    66266837             number_of_particles = prt_count(k,j,i)
    66276838             IF ( number_of_particles <= 0  .OR.                               &
     
    66356846!--          particle array. Tests showed that this simplified method can be
    66366847!--          used because it will only take place outside of cloudy grid
    6637 !--          boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned   
     6848!--          boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned
    66386849!--          and subsequent evaporated super droplets will be merged.
    66396850             DO  n = 1, number_of_particles-1
     
    66416852                     particles(n+1)%particle_mask                  .AND.       &
    66426853                     particles(n)%radius        <= radius_merge    .AND.       &
    6643                      particles(n)%weight_factor <= weight_factor_merge )       &   
     6854                     particles(n)%weight_factor <= weight_factor_merge )       &
    66446855                THEN
    66456856                   particles(n+1)%weight_factor  =                             &
     
    66526863                   deleted_particles          = deleted_particles + 1
    66536864                   merge_drp                  = merge_drp + 1
    6654                
     6865
    66556866                ENDIF
    66566867             ENDDO
     
    66586869       ENDDO
    66596870    ENDDO
    6660        
    6661        
     6871
     6872
    66626873    CALL cpu_log( log_point_s(81), 'lpm_merging', 'stop' )
    66636874
     
    71457356
    71467357       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
    7147 !     
     7358!
    71487359!--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit
    71497360!--    variables in structure particle_type (due to the calculation of par_size)
     
    71647375
    71657376       ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
    7166 !     
     7377!
    71677378!--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit
    71687379!--    variables in structure particle_type (due to the calculation of par_size)
     
    75167727             IF ( np_before_move <= 0 )  CYCLE
    75177728             particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move)
    7518              
     7729
    75197730             DO  n = 1, np_before_move
    75207731                i = particles_before_move(n)%x * ddx
     
    75787789!------------------------------------------------------------------------------!
    75797790 SUBROUTINE lpm_check_cfl 
    7580        
     7791
    75817792    IMPLICIT NONE
    7582    
     7793
    75837794    INTEGER(iwp)  ::  i !< running index, x-direction
    75847795    INTEGER(iwp)  ::  j !< running index, y-direction
     
    77597970! -----------
    77607971!> Routine for the whole processor
    7761 !> Sort all particles into the 8 respective subgrid boxes
     7972!> Sort all particles into the 8 respective subgrid boxes and free space of
     7973!> particles which has been marked for deletion
    77627974!------------------------------------------------------------------------------!
    7763    SUBROUTINE lpm_sort_in_subboxes
     7975   SUBROUTINE lpm_sort_and_delete
    77647976
    77657977       INTEGER(iwp) ::  i  !<
     
    77767988       INTEGER(iwp), DIMENSION(0:7) ::  sort_count  !<
    77777989
    7778        TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE ::  sort_particles  !<
    7779 
    7780        CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'start' )
    7781        DO  ip = nxl, nxr
    7782           DO  jp = nys, nyn
    7783              DO  kp = nzb+1, nzt
    7784                 number_of_particles = prt_count(kp,jp,ip)
    7785                 IF ( number_of_particles <= 0 )  CYCLE
    7786                 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    7787                      
    7788                 nn = 0
    7789                 sort_count = 0
    7790                 ALLOCATE( sort_particles(number_of_particles, 0:7) )
    7791                
    7792                 DO  n = 1, number_of_particles
    7793                    sort_index = 0
    7794 
    7795                    IF ( particles(n)%particle_mask )  THEN
    7796                       nn = nn + 1
    7797 !
    7798 !--                   Sorting particles with a binary scheme
    7799 !--                   sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox
    7800 !--                   sort_index=000_2=0_10 -> particle at the right,north,top subgridbox
    7801 !--                   For this the center of the gridbox is calculated
    7802                       i = (particles(n)%x + 0.5_wp * dx) * ddx
    7803                       j = (particles(n)%y + 0.5_wp * dy) * ddy
    7804                      
    7805                       IF ( i == ip )  sort_index = sort_index + 4
    7806                       IF ( j == jp )  sort_index = sort_index + 2                     
    7807                       IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1
    7808                      
    7809                       sort_count(sort_index) = sort_count(sort_index) + 1
    7810                       m = sort_count(sort_index)
    7811                       sort_particles(m,sort_index) = particles(n)
    7812                       sort_particles(m,sort_index)%block_nr = sort_index
    7813                    ENDIF
     7990       TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE ::  sort_particles    !<
     7991
     7992       CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' )
     7993       IF ( TRIM(particle_interpolation)  == 'trilinear' ) THEN
     7994          DO  ip = nxl, nxr
     7995             DO  jp = nys, nyn
     7996                DO  kp = nzb+1, nzt
     7997                   number_of_particles = prt_count(kp,jp,ip)
     7998                   IF ( number_of_particles <= 0 )  CYCLE
     7999                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     8000                   nn = 0
     8001                   sort_count = 0
     8002                   ALLOCATE( sort_particles(number_of_particles, 0:7) )
     8003
     8004                   DO  n = 1, number_of_particles
     8005                      sort_index = 0
     8006
     8007                      IF ( particles(n)%particle_mask )  THEN
     8008                         nn = nn + 1
     8009!
     8010!--                      Sorting particles with a binary scheme
     8011!--                      sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox
     8012!--                      sort_index=000_2=0_10 -> particle at the right,north,top subgridbox
     8013!--                      For this the center of the gridbox is calculated
     8014                         i = (particles(n)%x + 0.5_wp * dx) * ddx
     8015                         j = (particles(n)%y + 0.5_wp * dy) * ddy
     8016
     8017                         IF ( i == ip )  sort_index = sort_index + 4
     8018                         IF ( j == jp )  sort_index = sort_index + 2
     8019                         IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1
     8020
     8021                         sort_count(sort_index) = sort_count(sort_index) + 1
     8022                         m = sort_count(sort_index)
     8023                         sort_particles(m,sort_index) = particles(n)
     8024                         sort_particles(m,sort_index)%block_nr = sort_index
     8025                      ENDIF
     8026                   ENDDO
     8027
     8028                   nn = 0
     8029                   DO is = 0,7
     8030                      grid_particles(kp,jp,ip)%start_index(is) = nn + 1
     8031                      DO n = 1,sort_count(is)
     8032                         nn = nn + 1
     8033                         particles(nn) = sort_particles(n,is)
     8034                      ENDDO
     8035                      grid_particles(kp,jp,ip)%end_index(is) = nn
     8036                   ENDDO
     8037
     8038                   number_of_particles = nn
     8039                   prt_count(kp,jp,ip) = number_of_particles
     8040                   DEALLOCATE(sort_particles)
    78148041                ENDDO
    7815 
    7816                 nn = 0
    7817                 DO is = 0,7
    7818                    grid_particles(kp,jp,ip)%start_index(is) = nn + 1
    7819                    DO n = 1,sort_count(is)
    7820                       nn = nn + 1
    7821                       particles(nn) = sort_particles(n,is)
    7822                    ENDDO
    7823                    grid_particles(kp,jp,ip)%end_index(is) = nn
    7824                 ENDDO
    7825 
    7826                 number_of_particles = nn   
    7827                 prt_count(kp,jp,ip) = number_of_particles
    7828                 DEALLOCATE(sort_particles)
    78298042             ENDDO
    78308043          ENDDO
    7831        ENDDO
    7832        CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'stop' )
    7833 
    7834     END SUBROUTINE lpm_sort_in_subboxes
     8044
     8045!--    In case of the simple interpolation method the particles must not
     8046!--    be sorted in subboxes but particles marked for deletion must be
     8047!--    deleted and number of particles must be recalculated
     8048       ELSE
     8049
     8050          DO  ip = nxl, nxr
     8051             DO  jp = nys, nyn
     8052                DO  kp = nzb+1, nzt
     8053
     8054                   number_of_particles = prt_count(kp,jp,ip)
     8055                   IF ( number_of_particles <= 0 )  CYCLE
     8056                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
     8057!
     8058!--                repack particles array, i.e. delete particles and recalculate
     8059!--                number of particles
     8060                   CALL lpm_pack
     8061                   prt_count(kp,jp,ip) = number_of_particles
     8062                ENDDO
     8063             ENDDO
     8064          ENDDO
     8065       ENDIF
     8066       CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'stop' )
     8067
     8068    END SUBROUTINE lpm_sort_and_delete
    78358069
    78368070 
     
    78418075!------------------------------------------------------------------------------!
    78428076    SUBROUTINE lpm_pack
    7843    
     8077
    78448078       INTEGER(iwp) ::  n       !<
    78458079       INTEGER(iwp) ::  nn      !<
     
    78748108
    78758109    END SUBROUTINE lpm_pack
    7876                
    7877                
     8110
     8111
    78788112!------------------------------------------------------------------------------!
    78798113! Description:
     
    79048138                number_of_particles = prt_count(k,j,i)
    79058139                IF ( number_of_particles <= 0 )  CYCLE
    7906 
    79078140                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
    79088141
Note: See TracChangeset for help on using the changeset viewer.