Changeset 4143 for palm/trunk


Ignore:
Timestamp:
Aug 5, 2019 3:14:53 PM (5 years ago)
Author:
schwenkel
Message:

Rename variable and change select case to if statement

File:
1 edited

Legend:

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

    r4122 r4143  
    2525! -----------------
    2626! $Id$
     27! Rename variable and change select case to if statement
     28!
     29! 4122 2019-07-26 13:11:56Z schwenkel
    2730! Implement reset method as bottom boundary condition
    2831!
     
    278281    CHARACTER(LEN=5)  ::  splitting_mode     = 'const'            !< splitting mode
    279282
    280     CHARACTER(LEN=25) ::  particle_interpolation = 'trilinear'    !< interpolation method for calculatin the particle
     283    CHARACTER(LEN=25) ::  particle_advection_interpolation = 'trilinear' !< interpolation method for calculatin the particle
    281284
    282285    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step   
     
    315318    LOGICAL ::  use_kernel_tables = .FALSE.               !< parameter, which turns on the use of precalculated collision kernels
    316319    LOGICAL ::  write_particle_statistics = .FALSE.       !< namelist parameter (see documentation)
     320    LOGICAL ::  interpolation_simple_predictor = .FALSE.  !< flag for simple particle advection interpolation with predictor step
     321    LOGICAL ::  interpolation_simple_corrector = .FALSE.  !< flag for simple particle advection interpolation with corrector step
     322    LOGICAL ::  interpolation_trilinear = .FALSE.         !< flag for trilinear particle advection interpolation
    317323
    318324    LOGICAL, DIMENSION(max_number_of_particle_groups) ::   vertical_particle_advection = .TRUE. !< Switch for vertical particle transport
     
    555561       particles_per_point, &
    556562       particle_advection_start, &
    557        particle_interpolation, &
     563       particle_advection_interpolation, &
    558564       particle_maximum_age, &
    559565       pdx, &
     
    615621       particles_per_point, &
    616622       particle_advection_start, &
    617        particle_interpolation, &
     623       particle_advection_interpolation, &
    618624       particle_maximum_age, &
    619625       pdx, &
     
    862868!-- it can be combined as the sgs-velocites for active particles are
    863869!-- calculated differently, i.e. no subboxes are needed.
    864     IF ( .NOT. TRIM(particle_interpolation) == 'trilinear'  .AND.              &
     870    IF ( .NOT. TRIM(particle_advection_interpolation) == 'trilinear'  .AND.              &
    865871       use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
    866872          message_string = 'subrgrid scale velocities in combination with ' // &
     
    10121018       de_dx = 0.0_wp
    10131019       de_dy = 0.0_wp
    1014        de_dz = 0.0_wp             
    1015                  
     1020       de_dz = 0.0_wp
     1021
    10161022       sgs_wf_part = 1.0_wp / 3.0_wp
    10171023    ENDIF
     
    10631069       ENDDO
    10641070
     1071    ENDIF
     1072
     1073!
     1074!-- Check which particle interpolation method should be used
     1075    IF ( TRIM(particle_advection_interpolation)  ==  'trilinear' )  THEN
     1076       interpolation_simple_corrector = .FALSE.
     1077       interpolation_simple_predictor = .FALSE.
     1078       interpolation_trilinear        = .TRUE.
     1079    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_corrector' )  THEN
     1080       interpolation_simple_corrector = .TRUE.
     1081       interpolation_simple_predictor = .FALSE.
     1082       interpolation_trilinear        = .FALSE.
     1083    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_predictor' )  THEN
     1084       interpolation_simple_corrector = .FALSE.
     1085       interpolation_simple_predictor = .TRUE.
     1086       interpolation_trilinear        = .FALSE.
    10651087    ENDIF
    10661088
     
    13091331!
    13101332!--    Calculate initial_weighting_factor diagnostically
    1311        IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
    1312           initial_weighting_factor =  number_concentration  *                        &
     1333       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp )  THEN
     1334          initial_weighting_factor =  number_concentration  *                         &
    13131335                                      pdx(1) * pdy(1) * pdz(1)
    13141336       END IF
     
    13221344                DO WHILE ( pos_y <= psn(i) )
    13231345                   IF ( pos_y >= nys * dy  .AND.                  &
    1324                         pos_y <  ( nyn + 1 ) * dy  ) THEN
     1346                        pos_y <  ( nyn + 1 ) * dy  )  THEN
    13251347                      pos_x = psl(i)
    13261348               xloop: DO WHILE ( pos_x <= psr(i) )
    13271349                         IF ( pos_x >= nxl * dx  .AND.            &
    1328                               pos_x <  ( nxr + 1) * dx ) THEN
     1350                              pos_x <  ( nxr + 1) * dx )  THEN
    13291351                            DO  j = 1, particles_per_point
    13301352                               n = n + 1
     
    13661388!--                            In case of stretching the actual k index is found iteratively
    13671389                               IF ( dz_stretch_level .NE. -9999999.9_wp  .OR.           &
    1368                                     dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN
     1390                                    dz_stretch_level_start(1) .NE. -9999999.9_wp )  THEN
    13691391                                  kp = MINLOC( ABS( tmp_particle%z - zu ), DIM = 1 ) - 1
    13701392                               ELSE
     
    16261648!
    16271649!-- Set constants for different aerosol species
    1628     IF ( TRIM(aero_species) .EQ. 'nacl' ) THEN
     1650    IF ( TRIM(aero_species) .EQ. 'nacl' )  THEN
    16291651       molecular_weight_of_solute = 0.05844_wp
    16301652       rho_s                      = 2165.0_wp
    16311653       vanthoff                   = 2.0_wp
    1632     ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' ) THEN
     1654    ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' )  THEN
    16331655       molecular_weight_of_solute = 0.10406_wp
    16341656       rho_s                      = 1600.0_wp
    16351657       vanthoff                   = 1.37_wp
    1636     ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' ) THEN
     1658    ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' )  THEN
    16371659       molecular_weight_of_solute = 0.08004_wp
    16381660       rho_s                      = 1720.0_wp
     
    21532175
    21542176                      IF ( .NOT. first_loop_stride  .AND.  &
    2155                            grid_particles(k,j,i)%time_loop_done ) CYCLE
     2177                           grid_particles(k,j,i)%time_loop_done )  CYCLE
    21562178
    21572179                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     
    21872209!
    21882210!--                   Particle advection
    2189                       CALL lpm_advec(TRIM(particle_interpolation),i,j,k)
     2211                      CALL lpm_advec(i,j,k)
    21902212!
    21912213!--                   Particle reflection from walls. Only applied if the particles
     
    22432265!--          Apply splitting and merging algorithm
    22442266             IF ( cloud_droplets )  THEN
    2245                 IF ( splitting ) THEN
     2267                IF ( splitting )  THEN
    22462268                   CALL lpm_splitting
    22472269                ENDIF
    2248                 IF ( merging ) THEN
     2270                IF ( merging )  THEN
    22492271                   CALL lpm_merging
    22502272                ENDIF
     
    23922414    CALL close_file( 80 )
    23932415
    2394     IF ( number_of_particles > 0 ) THEN
     2416    IF ( number_of_particles > 0 )  THEN
    23952417        WRITE(9,*) 'number_of_particles ', number_of_particles,                &
    23962418                    current_timestep_number + 1, simulated_time + dt_3d
     
    31313153    WRITE ( 14 )  curvature_solution_effects
    31323154
     3155    CALL wrd_write_string( 'interpolation_simple_corrector' )
     3156    WRITE ( 14 )  interpolation_simple_corrector
     3157
     3158    CALL wrd_write_string( 'interpolation_simple_predictor' )
     3159    WRITE ( 14 )  interpolation_simple_predictor
     3160
     3161    CALL wrd_write_string( 'interpolation_trilinear' )
     3162    WRITE ( 14 )  interpolation_trilinear
     3163
    31333164 END SUBROUTINE lpm_wrd_global
    31343165 
     
    31523183       CASE ( 'curvature_solution_effects' )
    31533184          READ ( 13 )  curvature_solution_effects
    3154          
     3185
     3186       CASE ( 'interpolation_simple_corrector' )
     3187          READ ( 13 )  interpolation_simple_corrector
     3188
     3189       CASE ( 'interpolation_simple_predictor' )
     3190          READ ( 13 )  interpolation_simple_predictor
     3191
     3192       CASE ( 'interpolation_trilinear' )
     3193          READ ( 13 )  interpolation_trilinear
     3194
    31553195!          CASE ( 'global_paramter' )
    31563196!             READ ( 13 )  global_parameter
     
    31763216!> this submodule should be excluded as an own file.
    31773217!------------------------------------------------------------------------------!
    3178  SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp)
    3179 
    3180     CHARACTER (LEN=*), INTENT(IN) ::  interpolation_method !<
     3218 SUBROUTINE lpm_advec (ip,jp,kp)
     3219
    31813220    LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
    31823221
     
    33073346    dt_particle = dt_3d
    33083347
    3309 
    3310     SELECT CASE ( interpolation_method )
    3311 
    3312 !
    3313 !--    This case uses a simple interpolation method for the particle velocites,
    3314 !--    and applying a predictor-corrector method. @attention: for the corrector
    3315 !--    step the velocities of t(n+1) are required. However, at this moment of
    3316 !--    the time integration they are not free of divergence. This interpolation
    3317 !--    method is described in more detail in Grabowski et al., 2018 (GMD).
    3318        CASE ( 'simple_corrector' )
    3319 !
    3320 !--       Predictor step
    3321           kkw = kp - 1
    3322           DO n = 1, number_of_particles
    3323 
    3324              alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
    3325              u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
    3326 
    3327              beta  = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
    3328              v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
    3329 
    3330              gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
     3348!
     3349!-- This case uses a simple interpolation method for the particle velocites,
     3350!-- and applying a predictor-corrector method. @attention: for the corrector
     3351!-- step the velocities of t(n+1) are required. However, at this moment of
     3352!-- the time integration they are not free of divergence. This interpolation
     3353!-- method is described in more detail in Grabowski et al., 2018 (GMD).
     3354    IF ( interpolation_simple_corrector )  THEN
     3355!
     3356!--    Predictor step
     3357       kkw = kp - 1
     3358       DO n = 1, number_of_particles
     3359
     3360          alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3361          u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
     3362
     3363          beta  = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3364          v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
     3365
     3366          gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
     3367                            ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
     3368          w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
     3369
     3370       ENDDO
     3371!
     3372!--    Corrector step
     3373       DO n = 1, number_of_particles
     3374
     3375          IF ( .NOT. particles(n)%particle_mask )  CYCLE
     3376
     3377          DO nn = 1, iteration_steps
     3378
     3379!
     3380!--          Guess new position
     3381             xp = particles(n)%x + u_int(n) * dt_particle(n)
     3382             yp = particles(n)%y + v_int(n) * dt_particle(n)
     3383             zp = particles(n)%z + w_int(n) * dt_particle(n)
     3384!
     3385!--          x direction
     3386             i_next = FLOOR( xp * ddx , KIND=iwp)
     3387             alpha  = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3388!
     3389!--          y direction
     3390             j_next = FLOOR( yp * ddy )
     3391             beta   = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3392!
     3393!--          z_direction
     3394             k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0)
     3395             gamma = MAX( MIN( ( zp - zw(k_next) ) /                      &
     3396                               ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp )
     3397!
     3398!--          Calculate part of the corrector step
     3399             unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
     3400                     u_p(k_next+1, j_next,   i_next+1) * alpha
     3401
     3402             vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
     3403                     v_p(k_next+1, j_next+1, i_next  ) * beta
     3404
     3405             wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
     3406                     w_p(k_next+1, j_next, i_next  ) * gamma
     3407
     3408!
     3409!--          Calculate interpolated particle velocity with predictor
     3410!--          corrector step. u_int, v_int and w_int describes the part of
     3411!--          the predictor step. unext, vnext and wnext is the part of the
     3412!--          corrector step. The resulting new position is set below. The
     3413!--          implementation is based on Grabowski et al., 2018 (GMD).
     3414             u_int(n) = 0.5_wp * ( u_int(n) + unext )
     3415             v_int(n) = 0.5_wp * ( v_int(n) + vnext )
     3416             w_int(n) = 0.5_wp * ( w_int(n) + wnext )
     3417
     3418          ENDDO
     3419       ENDDO
     3420!
     3421!-- This case uses a simple interpolation method for the particle velocites,
     3422!-- and applying a predictor.
     3423    ELSEIF ( interpolation_simple_predictor )  THEN
     3424!
     3425!--    The particle position for the w velociy is based on the value of kp and kp-1
     3426       kkw = kp - 1
     3427       DO n = 1, number_of_particles
     3428          IF ( .NOT. particles(n)%particle_mask )  CYCLE
     3429
     3430          alpha    = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
     3431          u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
     3432
     3433          beta     = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
     3434          v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
     3435
     3436          gamma    = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
    33313437                               ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
    3332              w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
    3333 
     3438          w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
     3439       ENDDO
     3440!
     3441!-- The trilinear interpolation.
     3442    ELSEIF ( interpolation_trilinear )  THEN
     3443
     3444       start_index = grid_particles(kp,jp,ip)%start_index
     3445       end_index   = grid_particles(kp,jp,ip)%end_index
     3446
     3447       DO  nb = 0, 7
     3448!
     3449!--       Interpolate u velocity-component
     3450          i = ip
     3451          j = jp + block_offset(nb)%j_off
     3452          k = kp + block_offset(nb)%k_off
     3453
     3454          DO  n = start_index(nb), end_index(nb)
     3455!
     3456!--          Interpolation of the u velocity component onto particle position.
     3457!--          Particles are interpolation bi-linearly in the horizontal and a
     3458!--          linearly in the vertical. An exception is made for particles below
     3459!--          the first vertical grid level in case of a prandtl layer. In this
     3460!--          case the horizontal particle velocity components are determined using
     3461!--          Monin-Obukhov relations (if branch).
     3462!--          First, check if particle is located below first vertical grid level
     3463!--          above topography (Prandtl-layer height)
     3464!--          Determine vertical index of topography top
     3465             k_wall = get_topography_top_index_ji( jp, ip, 's' )
     3466
     3467             IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     3468!
     3469!--             Resolved-scale horizontal particle velocity is zero below z0.
     3470                IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
     3471                   u_int(n) = 0.0_wp
     3472                ELSE
     3473!
     3474!--                Determine the sublayer. Further used as index.
     3475                   height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
     3476                                        * REAL( number_of_sublayers, KIND=wp )    &
     3477                                        * d_z_p_z0
     3478!
     3479!--                Calculate LOG(z/z0) for exact particle height. Therefore,
     3480!--                interpolate linearly between precalculated logarithm.
     3481                   log_z_z0_int = log_z_z0(INT(height_p))                         &
     3482                                    + ( height_p - INT(height_p) )                &
     3483                                    * ( log_z_z0(INT(height_p)+1)                 &
     3484                                         - log_z_z0(INT(height_p))                &
     3485                                      )
     3486!
     3487!--                Get friction velocity and momentum flux from new surface data
     3488!--                types.
     3489                   IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3490                        surf_def_h(0)%end_index(jp,ip) )  THEN
     3491                      surf_start = surf_def_h(0)%start_index(jp,ip)
     3492!--                   Limit friction velocity. In narrow canyons or holes the
     3493!--                   friction velocity can become very small, resulting in a too
     3494!--                   large particle speed.
     3495                      us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
     3496                      usws_int  = surf_def_h(0)%usws(surf_start)
     3497                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
     3498                            surf_lsm_h%end_index(jp,ip) )  THEN
     3499                      surf_start = surf_lsm_h%start_index(jp,ip)
     3500                      us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
     3501                      usws_int  = surf_lsm_h%usws(surf_start)
     3502                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
     3503                            surf_usm_h%end_index(jp,ip) )  THEN
     3504                      surf_start = surf_usm_h%start_index(jp,ip)
     3505                      us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
     3506                      usws_int  = surf_usm_h%usws(surf_start)
     3507                   ENDIF
     3508!
     3509!--                Neutral solution is applied for all situations, e.g. also for
     3510!--                unstable and stable situations. Even though this is not exact
     3511!--                this saves a lot of CPU time since several calls of intrinsic
     3512!--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     3513!--                as sensitivity studies revealed no significant effect of
     3514!--                using the neutral solution also for un/stable situations.
     3515                   u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
     3516                               * log_z_z0_int - u_gtrans
     3517                ENDIF
     3518!
     3519!--          Particle above the first grid level. Bi-linear interpolation in the
     3520!--          horizontal and linear interpolation in the vertical direction.
     3521             ELSE
     3522                x  = xv(n) - i * dx
     3523                y  = yv(n) + ( 0.5_wp - j ) * dy
     3524                aa = x**2          + y**2
     3525                bb = ( dx - x )**2 + y**2
     3526                cc = x**2          + ( dy - y )**2
     3527                dd = ( dx - x )**2 + ( dy - y )**2
     3528                gg = aa + bb + cc + dd
     3529
     3530                u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
     3531                            + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
     3532                            u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     3533
     3534                IF ( k == nzt )  THEN
     3535                   u_int(n) = u_int_l
     3536                ELSE
     3537                   u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
     3538                               + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
     3539                               u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
     3540                   u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *            &
     3541                              ( u_int_u - u_int_l )
     3542                ENDIF
     3543             ENDIF
    33343544          ENDDO
    3335 
    3336 !
    3337 !--       Corrector step
    3338           DO n = 1, number_of_particles
    3339 
    3340              IF ( .NOT. particles(n)%particle_mask ) CYCLE
    3341 
    3342              DO nn = 1, iteration_steps
    3343 
    3344 !
    3345 !--             Guess new position
    3346                 xp = particles(n)%x + u_int(n) * dt_particle(n)
    3347                 yp = particles(n)%y + v_int(n) * dt_particle(n)
    3348                 zp = particles(n)%z + w_int(n) * dt_particle(n)
    3349 !
    3350 !--             x direction
    3351                 i_next = FLOOR( xp * ddx , KIND=iwp)
    3352                 alpha  = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp )
    3353 !
    3354 !--             y direction
    3355                 j_next = FLOOR( yp * ddy )
    3356                 beta   = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp )
    3357 !
    3358 !--             z_direction
    3359                 k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0)
    3360                 gamma = MAX( MIN( ( zp - zw(k_next) ) /                      &
    3361                                   ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp )
    3362 !
    3363 !--             Calculate part of the corrector step
    3364                 unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
    3365                         u_p(k_next+1, j_next,   i_next+1) * alpha
    3366 
    3367                 vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
    3368                         v_p(k_next+1, j_next+1, i_next  ) * beta
    3369 
    3370                 wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
    3371                         w_p(k_next+1, j_next, i_next  ) * gamma
    3372 
    3373 !
    3374 !--             Calculate interpolated particle velocity with predictor
    3375 !--             corrector step. u_int, v_int and w_int describes the part of
    3376 !--             the predictor step. unext, vnext and wnext is the part of the
    3377 !--             corrector step. The resulting new position is set below. The
    3378 !--             implementation is based on Grabowski et al., 2018 (GMD).
    3379                 u_int(n) = 0.5_wp * ( u_int(n) + unext )
    3380                 v_int(n) = 0.5_wp * ( v_int(n) + vnext )
    3381                 w_int(n) = 0.5_wp * ( w_int(n) + wnext )
    3382 
    3383              ENDDO
     3545!
     3546!--       Same procedure for interpolation of the v velocity-component
     3547          i = ip + block_offset(nb)%i_off
     3548          j = jp
     3549          k = kp + block_offset(nb)%k_off
     3550
     3551          DO  n = start_index(nb), end_index(nb)
     3552!
     3553!--          Determine vertical index of topography top
     3554             k_wall = get_topography_top_index_ji( jp,ip, 's' )
     3555
     3556             IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     3557                IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
     3558!
     3559!--                Resolved-scale horizontal particle velocity is zero below z0.
     3560                   v_int(n) = 0.0_wp
     3561                ELSE
     3562!
     3563!--                Determine the sublayer. Further used as index. Please note,
     3564!--                logarithmus can not be reused from above, as in in case of
     3565!--                topography particle on u-grid can be above surface-layer height,
     3566!--                whereas it can be below on v-grid.
     3567                   height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
     3568                                     * REAL( number_of_sublayers, KIND=wp )       &
     3569                                     * d_z_p_z0
     3570!
     3571!--                Calculate LOG(z/z0) for exact particle height. Therefore,
     3572!--                interpolate linearly between precalculated logarithm.
     3573                   log_z_z0_int = log_z_z0(INT(height_p))                         &
     3574                                    + ( height_p - INT(height_p) )                &
     3575                                    * ( log_z_z0(INT(height_p)+1)                 &
     3576                                         - log_z_z0(INT(height_p))                &
     3577                                      )
     3578!
     3579!--                Get friction velocity and momentum flux from new surface data
     3580!--                types.
     3581                   IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
     3582                        surf_def_h(0)%end_index(jp,ip) )  THEN
     3583                      surf_start = surf_def_h(0)%start_index(jp,ip)
     3584!--                   Limit friction velocity. In narrow canyons or holes the
     3585!--                   friction velocity can become very small, resulting in a too
     3586!--                   large particle speed.
     3587                      us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
     3588                      vsws_int  = surf_def_h(0)%vsws(surf_start)
     3589                   ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
     3590                            surf_lsm_h%end_index(jp,ip) )  THEN
     3591                      surf_start = surf_lsm_h%start_index(jp,ip)
     3592                      us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
     3593                      vsws_int  = surf_lsm_h%vsws(surf_start)
     3594                   ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
     3595                            surf_usm_h%end_index(jp,ip) )  THEN
     3596                      surf_start = surf_usm_h%start_index(jp,ip)
     3597                      us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
     3598                      vsws_int  = surf_usm_h%vsws(surf_start)
     3599                   ENDIF
     3600!
     3601!--                Neutral solution is applied for all situations, e.g. also for
     3602!--                unstable and stable situations. Even though this is not exact
     3603!--                this saves a lot of CPU time since several calls of intrinsic
     3604!--                FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     3605!--                as sensitivity studies revealed no significant effect of
     3606!--                using the neutral solution also for un/stable situations.
     3607                   v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
     3608                            * log_z_z0_int - v_gtrans
     3609
     3610                ENDIF
     3611             ELSE
     3612                x  = xv(n) + ( 0.5_wp - i ) * dx
     3613                y  = yv(n) - j * dy
     3614                aa = x**2          + y**2
     3615                bb = ( dx - x )**2 + y**2
     3616                cc = x**2          + ( dy - y )**2
     3617                dd = ( dx - x )**2 + ( dy - y )**2
     3618                gg = aa + bb + cc + dd
     3619
     3620                v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
     3621                          + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
     3622                          ) / ( 3.0_wp * gg ) - v_gtrans
     3623
     3624                IF ( k == nzt )  THEN
     3625                   v_int(n) = v_int_l
     3626                ELSE
     3627                   v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
     3628                             + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
     3629                             ) / ( 3.0_wp * gg ) - v_gtrans
     3630                   v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
     3631                                     ( v_int_u - v_int_l )
     3632                ENDIF
     3633             ENDIF
    33843634          ENDDO
    3385 
    3386 
    3387 !
    3388 !--    This case uses a simple interpolation method for the particle velocites,
    3389 !--    and applying a predictor.
    3390        CASE ( 'simple_predictor' )
    3391 !
    3392 !--       The particle position for the w velociy is based on the value of kp and kp-1
    3393           kkw = kp - 1
    3394           DO n = 1, number_of_particles
    3395              IF ( .NOT. particles(n)%particle_mask ) CYCLE
    3396 
    3397              alpha    = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
    3398              u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
    3399 
    3400              beta     = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
    3401              v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
    3402 
    3403              gamma    = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
    3404                                   ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
    3405              w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
     3635!
     3636!--       Same procedure for interpolation of the w velocity-component
     3637          i = ip + block_offset(nb)%i_off
     3638          j = jp + block_offset(nb)%j_off
     3639          k = kp - 1
     3640
     3641          DO  n = start_index(nb), end_index(nb)
     3642             IF ( vertical_particle_advection(particles(n)%group) )  THEN
     3643                x  = xv(n) + ( 0.5_wp - i ) * dx
     3644                y  = yv(n) + ( 0.5_wp - j ) * dy
     3645                aa = x**2          + y**2
     3646                bb = ( dx - x )**2 + y**2
     3647                cc = x**2          + ( dy - y )**2
     3648                dd = ( dx - x )**2 + ( dy - y )**2
     3649                gg = aa + bb + cc + dd
     3650
     3651                w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
     3652                          + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
     3653                          ) / ( 3.0_wp * gg )
     3654
     3655                IF ( k == nzt )  THEN
     3656                   w_int(n) = w_int_l
     3657                ELSE
     3658                   w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
     3659                               ( gg-bb ) * w(k+1,j,i+1) + &
     3660                               ( gg-cc ) * w(k+1,j+1,i) + &
     3661                               ( gg-dd ) * w(k+1,j+1,i+1) &
     3662                             ) / ( 3.0_wp * gg )
     3663                   w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
     3664                              ( w_int_u - w_int_l )
     3665                ENDIF
     3666             ELSE
     3667                w_int(n) = 0.0_wp
     3668             ENDIF
    34063669          ENDDO
    3407 !
    3408 !--    The trilinear interpolation.
    3409        CASE ( 'trilinear' )
    3410 
    3411           start_index = grid_particles(kp,jp,ip)%start_index
    3412           end_index   = grid_particles(kp,jp,ip)%end_index
    3413 
    3414           DO  nb = 0, 7
    3415 !
    3416 !--          Interpolate u velocity-component
    3417              i = ip
    3418              j = jp + block_offset(nb)%j_off
    3419              k = kp + block_offset(nb)%k_off
    3420 
    3421              DO  n = start_index(nb), end_index(nb)
    3422 !
    3423 !--             Interpolation of the u velocity component onto particle position.
    3424 !--             Particles are interpolation bi-linearly in the horizontal and a
    3425 !--             linearly in the vertical. An exception is made for particles below
    3426 !--             the first vertical grid level in case of a prandtl layer. In this
    3427 !--             case the horizontal particle velocity components are determined using
    3428 !--             Monin-Obukhov relations (if branch).
    3429 !--             First, check if particle is located below first vertical grid level
    3430 !--             above topography (Prandtl-layer height)
    3431 !--             Determine vertical index of topography top
    3432                 k_wall = get_topography_top_index_ji( jp, ip, 's' )
    3433 
    3434                 IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
    3435 !
    3436 !--                Resolved-scale horizontal particle velocity is zero below z0.
    3437                    IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
    3438                       u_int(n) = 0.0_wp
    3439                    ELSE
    3440 !
    3441 !--                   Determine the sublayer. Further used as index.
    3442                       height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
    3443                                            * REAL( number_of_sublayers, KIND=wp )    &
    3444                                            * d_z_p_z0
    3445 !
    3446 !--                   Calculate LOG(z/z0) for exact particle height. Therefore,
    3447 !--                   interpolate linearly between precalculated logarithm.
    3448                       log_z_z0_int = log_z_z0(INT(height_p))                         &
    3449                                        + ( height_p - INT(height_p) )                &
    3450                                        * ( log_z_z0(INT(height_p)+1)                 &
    3451                                             - log_z_z0(INT(height_p))                &
    3452                                          )
    3453 !
    3454 !--                   Get friction velocity and momentum flux from new surface data
    3455 !--                   types.
    3456                       IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
    3457                            surf_def_h(0)%end_index(jp,ip) )  THEN
    3458                          surf_start = surf_def_h(0)%start_index(jp,ip)
    3459 !--                      Limit friction velocity. In narrow canyons or holes the
    3460 !--                      friction velocity can become very small, resulting in a too
    3461 !--                      large particle speed.
    3462                          us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
    3463                          usws_int  = surf_def_h(0)%usws(surf_start)
    3464                       ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
    3465                                surf_lsm_h%end_index(jp,ip) )  THEN
    3466                          surf_start = surf_lsm_h%start_index(jp,ip)
    3467                          us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
    3468                          usws_int  = surf_lsm_h%usws(surf_start)
    3469                       ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
    3470                                surf_usm_h%end_index(jp,ip) )  THEN
    3471                          surf_start = surf_usm_h%start_index(jp,ip)
    3472                          us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
    3473                          usws_int  = surf_usm_h%usws(surf_start)
    3474                       ENDIF
    3475 !
    3476 !--                   Neutral solution is applied for all situations, e.g. also for
    3477 !--                   unstable and stable situations. Even though this is not exact
    3478 !--                   this saves a lot of CPU time since several calls of intrinsic
    3479 !--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    3480 !--                   as sensitivity studies revealed no significant effect of
    3481 !--                   using the neutral solution also for un/stable situations.
    3482                       u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
    3483                                   * log_z_z0_int - u_gtrans
    3484                    ENDIF
    3485 !
    3486 !--             Particle above the first grid level. Bi-linear interpolation in the
    3487 !--             horizontal and linear interpolation in the vertical direction.
    3488                 ELSE
    3489                    x  = xv(n) - i * dx
    3490                    y  = yv(n) + ( 0.5_wp - j ) * dy
    3491                    aa = x**2          + y**2
    3492                    bb = ( dx - x )**2 + y**2
    3493                    cc = x**2          + ( dy - y )**2
    3494                    dd = ( dx - x )**2 + ( dy - y )**2
    3495                    gg = aa + bb + cc + dd
    3496 
    3497                    u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
    3498                                + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
    3499                                u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
    3500 
    3501                    IF ( k == nzt )  THEN
    3502                       u_int(n) = u_int_l
    3503                    ELSE
    3504                       u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
    3505                                   + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
    3506                                   u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
    3507                       u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *            &
    3508                                  ( u_int_u - u_int_l )
    3509                    ENDIF
    3510                 ENDIF
    3511              ENDDO
    3512 !
    3513 !--          Same procedure for interpolation of the v velocity-component
    3514              i = ip + block_offset(nb)%i_off
    3515              j = jp
    3516              k = kp + block_offset(nb)%k_off
    3517 
    3518              DO  n = start_index(nb), end_index(nb)
    3519 !
    3520 !--             Determine vertical index of topography top
    3521                 k_wall = get_topography_top_index_ji( jp,ip, 's' )
    3522 
    3523                 IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
    3524                    IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
    3525 !
    3526 !--                   Resolved-scale horizontal particle velocity is zero below z0.
    3527                       v_int(n) = 0.0_wp
    3528                    ELSE
    3529 !
    3530 !--                   Determine the sublayer. Further used as index. Please note,
    3531 !--                   logarithmus can not be reused from above, as in in case of
    3532 !--                   topography particle on u-grid can be above surface-layer height,
    3533 !--                   whereas it can be below on v-grid.
    3534                       height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
    3535                                         * REAL( number_of_sublayers, KIND=wp )       &
    3536                                         * d_z_p_z0
    3537 !
    3538 !--                   Calculate LOG(z/z0) for exact particle height. Therefore,
    3539 !--                   interpolate linearly between precalculated logarithm.
    3540                       log_z_z0_int = log_z_z0(INT(height_p))                         &
    3541                                        + ( height_p - INT(height_p) )                &
    3542                                        * ( log_z_z0(INT(height_p)+1)                 &
    3543                                             - log_z_z0(INT(height_p))                &
    3544                                          )
    3545 !
    3546 !--                   Get friction velocity and momentum flux from new surface data
    3547 !--                   types.
    3548                       IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
    3549                            surf_def_h(0)%end_index(jp,ip) )  THEN
    3550                          surf_start = surf_def_h(0)%start_index(jp,ip)
    3551 !--                      Limit friction velocity. In narrow canyons or holes the
    3552 !--                      friction velocity can become very small, resulting in a too
    3553 !--                      large particle speed.
    3554                          us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
    3555                          vsws_int  = surf_def_h(0)%vsws(surf_start)
    3556                       ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
    3557                                surf_lsm_h%end_index(jp,ip) )  THEN
    3558                          surf_start = surf_lsm_h%start_index(jp,ip)
    3559                          us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
    3560                          vsws_int  = surf_lsm_h%vsws(surf_start)
    3561                       ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
    3562                                surf_usm_h%end_index(jp,ip) )  THEN
    3563                          surf_start = surf_usm_h%start_index(jp,ip)
    3564                          us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
    3565                          vsws_int  = surf_usm_h%vsws(surf_start)
    3566                       ENDIF
    3567 !
    3568 !--                   Neutral solution is applied for all situations, e.g. also for
    3569 !--                   unstable and stable situations. Even though this is not exact
    3570 !--                   this saves a lot of CPU time since several calls of intrinsic
    3571 !--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
    3572 !--                   as sensitivity studies revealed no significant effect of
    3573 !--                   using the neutral solution also for un/stable situations.
    3574                       v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
    3575                                * log_z_z0_int - v_gtrans
    3576 
    3577                    ENDIF
    3578                 ELSE
    3579                    x  = xv(n) + ( 0.5_wp - i ) * dx
    3580                    y  = yv(n) - j * dy
    3581                    aa = x**2          + y**2
    3582                    bb = ( dx - x )**2 + y**2
    3583                    cc = x**2          + ( dy - y )**2
    3584                    dd = ( dx - x )**2 + ( dy - y )**2
    3585                    gg = aa + bb + cc + dd
    3586 
    3587                    v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
    3588                              + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
    3589                              ) / ( 3.0_wp * gg ) - v_gtrans
    3590 
    3591                    IF ( k == nzt )  THEN
    3592                       v_int(n) = v_int_l
    3593                    ELSE
    3594                       v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
    3595                                 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
    3596                                 ) / ( 3.0_wp * gg ) - v_gtrans
    3597                       v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
    3598                                         ( v_int_u - v_int_l )
    3599                    ENDIF
    3600                 ENDIF
    3601              ENDDO
    3602 !
    3603 !--          Same procedure for interpolation of the w velocity-component
    3604              i = ip + block_offset(nb)%i_off
    3605              j = jp + block_offset(nb)%j_off
    3606              k = kp - 1
    3607 
    3608              DO  n = start_index(nb), end_index(nb)
    3609                 IF ( vertical_particle_advection(particles(n)%group) )  THEN
    3610                    x  = xv(n) + ( 0.5_wp - i ) * dx
    3611                    y  = yv(n) + ( 0.5_wp - j ) * dy
    3612                    aa = x**2          + y**2
    3613                    bb = ( dx - x )**2 + y**2
    3614                    cc = x**2          + ( dy - y )**2
    3615                    dd = ( dx - x )**2 + ( dy - y )**2
    3616                    gg = aa + bb + cc + dd
    3617 
    3618                    w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
    3619                              + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
    3620                              ) / ( 3.0_wp * gg )
    3621 
    3622                    IF ( k == nzt )  THEN
    3623                       w_int(n) = w_int_l
    3624                    ELSE
    3625                       w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
    3626                                   ( gg-bb ) * w(k+1,j,i+1) + &
    3627                                   ( gg-cc ) * w(k+1,j+1,i) + &
    3628                                   ( gg-dd ) * w(k+1,j+1,i+1) &
    3629                                 ) / ( 3.0_wp * gg )
    3630                       w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
    3631                                  ( w_int_u - w_int_l )
    3632                    ENDIF
    3633                 ELSE
    3634                    w_int(n) = 0.0_wp
    3635                 ENDIF
    3636              ENDDO
    3637           ENDDO
    3638 
    3639        CASE DEFAULT
    3640           WRITE( message_string, * )  'unknown particle velocity interpolation method = "',  &
    3641                                        TRIM( interpolation_method ), '"'
    3642           CALL message( 'lpm_advec', 'PA0660', 1, 2, 0, 6, 0 )
    3643 
    3644     END SELECT
     3670       ENDDO
     3671    ENDIF
    36453672
    36463673!-- Interpolate and calculate quantities needed for calculating the SGS
     
    36533680!
    36543681!--       In case of topography check if subbox is adjacent to a wall
    3655           IF ( .NOT. topography == 'flat' ) THEN
     3682          IF ( .NOT. topography == 'flat' )  THEN
    36563683             i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
    36573684             j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
     
    36643691             ENDIF
    36653692          ENDIF
    3666           IF ( subbox_at_wall ) THEN
     3693          IF ( subbox_at_wall )  THEN
    36673694             e_int(start_index(nb):end_index(nb))     = e(kp,jp,ip)
    36683695             diss_int(start_index(nb):end_index(nb))  = diss(kp,jp,ip)
     
    39603987                (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n))  .OR.   &
    39613988                 ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n))  .OR.   &
    3962                  ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN
     3989                 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n))))  THEN
    39633990
    39643991                dt_particle(n) = 0.9_wp * MIN(                                 &
     
    40274054!--    number_of_particles. @todo find a more generic way to write this loop or
    40284055!--    delete trilinear interpolation
    4029        IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4056       IF ( interpolation_trilinear )  THEN
    40304057          subbox_start = 0
    40314058          subbox_end   = 7
     
    40394066!--    from 1 to 1.
    40404067       DO  nb = subbox_start, subbox_end
    4041           IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4068          IF ( interpolation_trilinear )  THEN
    40424069             particle_start = start_index(nb)
    40434070             particle_end   = end_index(nb)
     
    41414168!--    Decide whether the particle loop runs over the subboxes or only over 1,
    41424169!--    number_of_particles. This depends on the selected interpolation method.
    4143        IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4170       IF ( interpolation_trilinear )  THEN
    41444171          subbox_start = 0
    41454172          subbox_end   = 7
     
    41524179!--    from 1 to 1.
    41534180       DO  nb = subbox_start, subbox_end
    4154           IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4181          IF ( interpolation_trilinear )  THEN
    41554182             particle_start = start_index(nb)
    41564183             particle_end   = end_index(nb)
     
    42454272!-- Decide whether the particle loop runs over the subboxes or only over 1,
    42464273!-- number_of_particles. This depends on the selected interpolation method.
    4247     IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4274    IF ( interpolation_trilinear )  THEN
    42484275       subbox_start = 0
    42494276       subbox_end   = 7
     
    42534280    ENDIF
    42544281    DO  nb = subbox_start, subbox_end
    4255        IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
     4282       IF ( interpolation_trilinear )  THEN
    42564283          particle_start = start_index(nb)
    42574284          particle_end   = end_index(nb)
     
    45054532                particles(n)%z       = reset_top *  ( 1.0  + ( ran_val / 10.0_wp) )
    45064533                particles(n)%speed_z = 0.0_wp
    4507                 IF ( curvature_solution_effects ) THEN
     4534                IF ( curvature_solution_effects )  THEN
    45084535                   particles(n)%radius = particles(n)%aux1
    45094536                ELSE
     
    45664593          ENDIF
    45674594
    4568           IF ( prt_z > pos_z_old ) THEN
     4595          IF ( prt_z > pos_z_old )  THEN
    45694596             zwall = zw(k)
    45704597          ELSE
     
    46484675          reach_y(t_index) = .FALSE.
    46494676          reach_z(t_index) = .TRUE.
    4650           IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN
     4677          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp)  THEN
    46514678             t_index      = t_index + 1
    46524679             cross_wall_z = .TRUE.
     
    50015028    ENDDO
    50025029
    5003     IF( .NOT. curvature_solution_effects ) then
     5030    IF( .NOT. curvature_solution_effects )  THEN
    50045031!
    50055032!--    Use analytic model for diffusional growth including gas-kinetic
     
    61366163
    61376164    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
    6138    
     6165
    61396166    LOGICAL, SAVE ::  first = .TRUE. !<
    61406167
     
    61496176    REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
    61506177    REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
    6151    
     6178
    61526179    REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_100 !<
    61536180    REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_400 !<
     
    70107037          DO  kp = nzb+1, nzt
    70117038             number_of_particles = prt_count(kp,jp,ip)
    7012              IF ( number_of_particles <= 0 ) CYCLE
     7039             IF ( number_of_particles <= 0 )  CYCLE
    70137040             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    70147041             DO  n = 1, number_of_particles
     
    75897616!--    In case of stretching the actual k index must be found
    75907617       IF ( dz_stretch_level .NE. -9999999.9_wp  .OR.         &
    7591             dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN
     7618            dz_stretch_level_start(1) .NE. -9999999.9_wp )  THEN
    75927619          kp = MINLOC( ABS( particle_array(n)%z - zu ), DIM = 1 ) - 1
    75937620       ELSE
     
    79928019! -----------
    79938020!> Routine for the whole processor
    7994 !> Sort all particles into the 8 respective subgrid boxes and free space of
    7995 !> particles which has been marked for deletion
     8021!> Sort all particles into the 8 respective subgrid boxes (in case of trilinear
     8022!> interpolation method) and free space of particles which has been marked for
     8023!> deletion.
    79968024!------------------------------------------------------------------------------!
    79978025   SUBROUTINE lpm_sort_and_delete
     
    80138041
    80148042       CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' )
    8015        IF ( TRIM(particle_interpolation)  == 'trilinear' ) THEN
     8043       IF ( interpolation_trilinear ) THEN
    80168044          DO  ip = nxl, nxr
    80178045             DO  jp = nys, nyn
     
    80478075                      ENDIF
    80488076                   ENDDO
    8049 
     8077!
     8078!--                Delete and resort particles by overwritting and set
     8079!--                the number_of_particles to the actual value.
    80508080                   nn = 0
    80518081                   DO is = 0,7
     
    80668096
    80678097!--    In case of the simple interpolation method the particles must not
    8068 !--    be sorted in subboxes but particles marked for deletion must be
    8069 !--    deleted and number of particles must be recalculated
     8098!--    be sorted in subboxes. Particles marked for deletion however, must be
     8099!--    deleted and number of particles must be recalculated as it is also
     8100!--    done for the trilinear particle advection interpolation method.
    80708101       ELSE
    80718102
     
    80788109                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
    80798110!
    8080 !--                repack particles array, i.e. delete particles and recalculate
     8111!--                Repack particles array, i.e. delete particles and recalculate
    80818112!--                number of particles
    80828113                   CALL lpm_pack
Note: See TracChangeset for help on using the changeset viewer.