Index: /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4144)
+++ /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4145)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Some reformatting
+! 
+! 4144 2019-08-06 09:11:47Z raasch
 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
 ! 
@@ -387,6 +390,6 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
 
-    INTEGER(iwp), PARAMETER         :: PHASE_INIT    = 1  !<
-    INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2  !<
+    INTEGER(iwp), PARAMETER         ::  PHASE_INIT    = 1  !<
+    INTEGER(iwp), PARAMETER, PUBLIC ::  PHASE_RELEASE = 2  !<
 
     SAVE
@@ -871,6 +874,6 @@
 !-- it can be combined as the sgs-velocites for active particles are
 !-- calculated differently, i.e. no subboxes are needed.
-    IF ( .NOT. TRIM(particle_advection_interpolation) == 'trilinear'  .AND.              &
-       use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
+    IF ( .NOT. TRIM( particle_advection_interpolation ) == 'trilinear'  .AND.  &
+       use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
           message_string = 'subrgrid scale velocities in combination with ' // &
                            'simple interpolation method is not '            // &
@@ -1076,13 +1079,13 @@
 !
 !-- Check which particle interpolation method should be used
-    IF ( TRIM(particle_advection_interpolation)  ==  'trilinear' )  THEN
+    IF ( TRIM( particle_advection_interpolation )  ==  'trilinear' )  THEN
        interpolation_simple_corrector = .FALSE.
        interpolation_simple_predictor = .FALSE.
        interpolation_trilinear        = .TRUE.
-    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_corrector' )  THEN
+    ELSEIF ( TRIM( particle_advection_interpolation )  ==  'simple_corrector' )  THEN
        interpolation_simple_corrector = .TRUE.
        interpolation_simple_predictor = .FALSE.
        interpolation_trilinear        = .FALSE.
-    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_predictor' )  THEN
+    ELSEIF ( TRIM( particle_advection_interpolation )  ==  'simple_predictor' )  THEN
        interpolation_simple_corrector = .FALSE.
        interpolation_simple_predictor = .TRUE.
@@ -1334,6 +1337,6 @@
 !
 !--    Calculate initial_weighting_factor diagnostically
-       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp )  THEN
-          initial_weighting_factor =  number_concentration  *                         &
+       IF ( number_concentration /= -1.0_wp  .AND.  number_concentration > 0.0_wp )  THEN
+          initial_weighting_factor =  number_concentration  *                           &
                                       pdx(1) * pdy(1) * pdz(1)
        END IF
@@ -1471,5 +1474,5 @@
                             alloc_factor / 100.0_wp ) ), 1 )
                          IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) )  THEN
-                            CALL realloc_particles_array(ip,jp,kp,alloc_size)
+                            CALL realloc_particles_array( ip, jp, kp, alloc_size )
                          ENDIF
                       ENDIF
@@ -1510,5 +1513,5 @@
 !-- Initialize aerosol background spectrum
     IF ( curvature_solution_effects )  THEN
-       CALL lpm_init_aerosols(local_start)
+       CALL lpm_init_aerosols( local_start )
     ENDIF
 !
@@ -1629,35 +1632,35 @@
  SUBROUTINE lpm_init_aerosols(local_start)
 
-    REAL(wp)  :: afactor            !< curvature effects
-    REAL(wp)  :: bfactor            !< solute effects
-    REAL(wp)  :: dlogr              !< logarithmic width of radius bin
-    REAL(wp)  :: e_a                !< vapor pressure
-    REAL(wp)  :: e_s                !< saturation vapor pressure
-    REAL(wp)  :: rmin = 0.005e-6_wp !< minimum aerosol radius
-    REAL(wp)  :: rmax = 10.0e-6_wp  !< maximum aerosol radius
-    REAL(wp)  :: r_mid              !< mean radius of bin
-    REAL(wp)  :: r_l                !< left radius of bin
-    REAL(wp)  :: r_r                !< right radius of bin
-    REAL(wp)  :: sigma              !< surface tension
-    REAL(wp)  :: t_int              !< temperature
+    REAL(wp) ::  afactor            !< curvature effects
+    REAL(wp) ::  bfactor            !< solute effects
+    REAL(wp) ::  dlogr              !< logarithmic width of radius bin
+    REAL(wp) ::  e_a                !< vapor pressure
+    REAL(wp) ::  e_s                !< saturation vapor pressure
+    REAL(wp) ::  rmin = 0.005e-6_wp !< minimum aerosol radius
+    REAL(wp) ::  rmax = 10.0e-6_wp  !< maximum aerosol radius
+    REAL(wp) ::  r_mid              !< mean radius of bin
+    REAL(wp) ::  r_l                !< left radius of bin
+    REAL(wp) ::  r_r                !< right radius of bin
+    REAL(wp) ::  sigma              !< surface tension
+    REAL(wp) ::  t_int              !< temperature
 
     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
 
-    INTEGER(iwp)  :: n              !<
-    INTEGER(iwp)  :: ip             !<
-    INTEGER(iwp)  :: jp             !<
-    INTEGER(iwp)  :: kp             !<
+    INTEGER(iwp) ::  n              !<
+    INTEGER(iwp) ::  ip             !<
+    INTEGER(iwp) ::  jp             !<
+    INTEGER(iwp) ::  kp             !<
 
 !
 !-- Set constants for different aerosol species
-    IF ( TRIM(aero_species) == 'nacl' )  THEN
+    IF ( TRIM( aero_species ) == 'nacl' )  THEN
        molecular_weight_of_solute = 0.05844_wp 
        rho_s                      = 2165.0_wp
        vanthoff                   = 2.0_wp
-    ELSEIF ( TRIM(aero_species) == 'c3h4o4' )  THEN
+    ELSEIF ( TRIM( aero_species ) == 'c3h4o4' )  THEN
        molecular_weight_of_solute = 0.10406_wp 
        rho_s                      = 1600.0_wp
        vanthoff                   = 1.37_wp
-    ELSEIF ( TRIM(aero_species) == 'nh4o3' )  THEN
+    ELSEIF ( TRIM( aero_species ) == 'nh4o3' )  THEN
        molecular_weight_of_solute = 0.08004_wp 
        rho_s                      = 1720.0_wp
@@ -1671,33 +1674,33 @@
 !-- The following typical aerosol spectra are taken from Jaenicke (1993):
 !-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
-    IF ( TRIM(aero_type) == 'polar' )  THEN
+    IF ( TRIM( aero_type ) == 'polar' )  THEN
        na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp
        rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp
        log_sigma = (/ 0.245, 0.300, 0.291 /)
-    ELSEIF ( TRIM(aero_type) == 'background' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'background' )  THEN
        na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp
        rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp
        log_sigma = (/ 0.645, 0.253, 0.425 /)
-    ELSEIF ( TRIM(aero_type) == 'maritime' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'maritime' )  THEN
        na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp
        rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp
        log_sigma = (/ 0.657, 0.210, 0.396 /)
-    ELSEIF ( TRIM(aero_type) == 'continental' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'continental' )  THEN
        na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp
        rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp
        log_sigma = (/ 0.161, 0.217, 0.380 /)
-    ELSEIF ( TRIM(aero_type) == 'desert' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'desert' )  THEN
        na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp
        rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp
        log_sigma = (/ 0.247, 0.770, 0.438 /)
-    ELSEIF ( TRIM(aero_type) == 'rural' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'rural' )  THEN
        na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp
        rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp
        log_sigma = (/ 0.225, 0.557, 0.266 /)
-    ELSEIF ( TRIM(aero_type) == 'urban' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'urban' )  THEN
        na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp
        rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp
        log_sigma = (/ 0.245, 0.666, 0.337 /)
-    ELSEIF ( TRIM(aero_type) == 'user' )  THEN
+    ELSEIF ( TRIM( aero_type ) == 'user' )  THEN
        CONTINUE
     ELSE
@@ -2117,6 +2120,6 @@
 !
 !--       If necessary, release new set of particles
-          IF ( ( simulated_time - last_particle_release_time ) >= dt_prel  .AND. end_time_prel > simulated_time ) &
-          THEN
+          IF ( ( simulated_time - last_particle_release_time ) >= dt_prel  .AND.     &
+                 end_time_prel > simulated_time )  THEN
              DO WHILE ( ( simulated_time - last_particle_release_time ) >= dt_prel )
                 CALL lpm_create_particle( PHASE_RELEASE )
@@ -2212,10 +2215,10 @@
 !
 !--                   Particle advection
-                      CALL lpm_advec(i,j,k)
+                      CALL lpm_advec( i, j, k )
 !
 !--                   Particle reflection from walls. Only applied if the particles 
 !--                   are in the vertical range of the topography. (Here, some
 !--                   optimization is still possible.)
-                      IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
+                      IF ( topography /= 'flat'  .AND.  k < nzb_max + 2 )  THEN
                          CALL  lpm_boundary_conds( 'walls', i, j, k )
                       ENDIF
@@ -2223,5 +2226,5 @@
 !--                   User-defined actions after the calculation of the new particle 
 !--                   position
-                      CALL user_lpm_advec(i,j,k)
+                      CALL user_lpm_advec( i, j, k )
 !
 !--                   Apply boundary conditions to those particles that have crossed 
@@ -2284,5 +2287,5 @@
 !
 !--          IF .FALSE., lpm_sort_and_delete is done inside pcmp
-             IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run )   THEN   
+             IF ( .NOT. dt_3d_reached  .OR.  .NOT. nested_run )   THEN
 !
 !--             Pack particles (eliminate those marked for deletion),
@@ -2386,8 +2389,8 @@
  SUBROUTINE lpm_write_exchange_statistics
 
-    INTEGER(iwp) :: ip         !<
-    INTEGER(iwp) :: jp         !<
-    INTEGER(iwp) :: kp         !<
-    INTEGER(iwp) :: tot_number_of_particles
+    INTEGER(iwp) ::  ip         !<
+    INTEGER(iwp) ::  jp         !<
+    INTEGER(iwp) ::  kp         !<
+    INTEGER(iwp) ::  tot_number_of_particles !<
 
 !
@@ -2462,5 +2465,5 @@
     WRITE ( 85 )  simulated_time
     WRITE ( 85 )  prt_count
-          
+
     DO  ip = nxl, nxr
        DO  jp = nys, nyn
@@ -2898,10 +2901,10 @@
     CHARACTER (LEN=10) ::  version_on_file            !<
 
-    INTEGER(iwp) :: alloc_size !<
-    INTEGER(iwp) :: ip         !<
-    INTEGER(iwp) :: jp         !<
-    INTEGER(iwp) :: kp         !<
-
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles !<
+    INTEGER(iwp) ::  alloc_size !<
+    INTEGER(iwp) ::  ip         !<
+    INTEGER(iwp) ::  jp         !<
+    INTEGER(iwp) ::  kp         !<
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles !<
 
 !
@@ -3016,5 +3019,5 @@
     LOGICAL, INTENT(OUT)  ::  found
 
-    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::  tmp_3d   !<
 
 
@@ -3359,5 +3362,5 @@
 !--    Predictor step
        kkw = kp - 1
-       DO n = 1, number_of_particles
+       DO  n = 1, number_of_particles
 
           alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
@@ -3374,9 +3377,9 @@
 !
 !--    Corrector step
-       DO n = 1, number_of_particles
+       DO  n = 1, number_of_particles
 
           IF ( .NOT. particles(n)%particle_mask )  CYCLE
 
-          DO nn = 1, iteration_steps
+          DO  nn = 1, iteration_steps
 
 !
@@ -3428,5 +3431,5 @@
 !--    The particle position for the w velociy is based on the value of kp and kp-1
        kkw = kp - 1
-       DO n = 1, number_of_particles
+       DO  n = 1, number_of_particles
           IF ( .NOT. particles(n)%particle_mask )  CYCLE
 
@@ -3830,5 +3833,5 @@
        ENDDO
 
-       DO nb = 0,7
+       DO  nb = 0,7
           i = ip + block_offset(nb)%i_off
           j = jp + block_offset(nb)%j_off
@@ -3963,13 +3966,13 @@
                 ENDIF
 
-                CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),    &
                                         de_dx_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,1), term_1_2(n) )
 
-                CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),    &
                                         de_dy_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,2), term_1_2(n) )
 
-                CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),    &
                                         de_dz_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,3), term_1_2(n) )
@@ -4011,13 +4014,13 @@
                 ENDIF
 
-                CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),    &
                                         de_dx_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,1), term_1_2(n) )
 
-                CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),    &
                                         de_dy_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,2), term_1_2(n) )
 
-                CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),& 
+                CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),    &
                                         de_dz_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,3), term_1_2(n) )
@@ -4420,7 +4423,7 @@
     INTEGER(iwp) ::  tmp_z          !< dummy for sorting algorithm
 
-    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 
-    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
-    INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
+    INTEGER(iwp), DIMENSION(0:10) ::  x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
+    INTEGER(iwp), DIMENSION(0:10) ::  y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
+    INTEGER(iwp), DIMENSION(0:10) ::  z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
 
     LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
@@ -4556,7 +4559,7 @@
           i2 = particles(n)%x * ddx
           j2 = particles(n)%y * ddy
-          IF (zw(k)   < particles(n)%z ) k2 = k + 1
-          IF (zw(k)   > particles(n)%z .AND. zw(k-1) < particles(n)%z ) k2 = k
-          IF (zw(k-1) > particles(n)%z ) k2 = k - 1 
+          IF ( zw(k)   < particles(n)%z ) k2 = k + 1
+          IF ( zw(k)   > particles(n)%z  .AND.  zw(k-1) < particles(n)%z ) k2 = k
+          IF ( zw(k-1) > particles(n)%z ) k2 = k - 1
 !
 !--       Save current particle positions
@@ -4590,5 +4593,5 @@
 !--       Wall to the north
           IF ( prt_y > pos_y_old )  THEN
-             ywall = ( j1 +1 ) * dy
+             ywall = ( j1 + 1 ) * dy
 !--       Wall to the south
           ELSE
@@ -4645,5 +4648,5 @@
 !--       Store these values only if particle really reaches any wall. t must 
 !--       be in a interval between [0:1].
-          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
+          IF ( t(t_index) <= 1.0_wp  .AND.  t(t_index) >= 0.0_wp )  THEN
              t_index      = t_index + 1
              cross_wall_x = .TRUE.
@@ -4661,5 +4664,5 @@
           reach_y(t_index) = .TRUE.
           reach_z(t_index) = .FALSE.
-          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
+          IF ( t(t_index) <= 1.0_wp  .AND.  t(t_index) >= 0.0_wp )  THEN
              t_index      = t_index + 1
              cross_wall_y = .TRUE.
@@ -4678,5 +4681,5 @@
           reach_y(t_index) = .FALSE.
           reach_z(t_index) = .TRUE.
-          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp)  THEN
+          IF( t(t_index) <= 1.0_wp  .AND.  t(t_index) >= 0.0_wp)  THEN
              t_index      = t_index + 1
              cross_wall_z = .TRUE.
@@ -4686,5 +4689,5 @@
 !
 !--       Carry out reflection only if particle reaches any wall
-          IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z )  THEN
+          IF ( cross_wall_x  .OR.  cross_wall_y  .OR.  cross_wall_z )  THEN
 !
 !--          Sort fractional timesteps in ascending order. Also sort the 
@@ -4885,5 +4888,5 @@
 !--          If a particle was reflected, calculate final position from last
 !--          intermediate position.
-             IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
+             IF ( reflect_x  .OR.  reflect_y  .OR.  reflect_z )  THEN
 
                 particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
@@ -4922,8 +4925,8 @@
  SUBROUTINE lpm_droplet_condensation (i,j,k)
 
-    INTEGER(iwp), INTENT(IN) :: i              !<
-    INTEGER(iwp), INTENT(IN) :: j              !<
-    INTEGER(iwp), INTENT(IN) :: k              !<
-    INTEGER(iwp) :: n                          !<
+    INTEGER(iwp), INTENT(IN) ::  i              !<
+    INTEGER(iwp), INTENT(IN) ::  j              !<
+    INTEGER(iwp), INTENT(IN) ::  k              !<
+    INTEGER(iwp) ::  n                          !<
 
     REAL(wp) ::  afactor                       !< curvature effects
@@ -4954,8 +4957,8 @@
 !
 !-- Parameters for Rosenbrock method (see Verwer et al., 1999)
-    REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
-    REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
-    REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
-    REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
+    REAL(wp), PARAMETER ::  prec = 1.0E-3_wp     !< precision of Rosenbrock solution
+    REAL(wp), PARAMETER ::  q_increase = 1.5_wp  !< increase factor in timestep
+    REAL(wp), PARAMETER ::  q_decrease = 0.9_wp  !< decrease factor in timestep
+    REAL(wp), PARAMETER ::  gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
 !
 !-- Parameters for terminal velocity
@@ -5152,6 +5155,6 @@
 !--    case, the model timestep might be too long.
        delta_r = new_r(n) - particles(n)%radius
-       IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
-          WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,             &
+       IF ( delta_r < 0.0_wp  .AND.  new_r(n) < 0.0_wp )  THEN
+          WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,                &
                        ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,               &
                        '&delta_r=',delta_r,                                    &
@@ -6158,10 +6161,10 @@
  SUBROUTINE turb_enhance_eff
 
-    INTEGER(iwp) :: i  !<
-    INTEGER(iwp) :: iq !<
-    INTEGER(iwp) :: ir !<
-    INTEGER(iwp) :: j  !<
-    INTEGER(iwp) :: k  !<
-    INTEGER(iwp) :: kk !<
+    INTEGER(iwp) ::  i  !<
+    INTEGER(iwp) ::  iq !<
+    INTEGER(iwp) ::  ir !<
+    INTEGER(iwp) ::  j  !<
+    INTEGER(iwp) ::  k  !<
+    INTEGER(iwp) ::  kk !<
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
@@ -6348,7 +6351,7 @@
     INTEGER(iwp) ::  np               !< 
     INTEGER(iwp) ::  old_size         !< old particle array size
-    
+
     INTEGER(iwp), PARAMETER ::  n_max = 100 !< number of radii bin for splitting functions    
-    
+
     LOGICAL ::  first_loop_stride_sp = .TRUE. !< flag to calculate constants only once
 
@@ -6384,5 +6387,5 @@
     REAL(wp), DIMENSION(0:n_max-1) ::  r_bin_mid  !< mass weighted mean radius of a bin
     REAL(wp), DIMENSION(0:n_max)   ::  r_bin      !< boundaries of a radius bin
-    
+
     TYPE(particle_type) ::  tmp_particle   !< temporary particle TYPE
 
@@ -6442,5 +6445,5 @@
 !--                   Reallocate particle array if necessary. 
                       IF ( new_size > SIZE(particles) )  THEN 
-                         CALL realloc_particles_array(i,j,k,new_size)
+                         CALL realloc_particles_array( i, j, k, new_size )
                       ENDIF
                       old_size = prt_count(k,j,i)
@@ -6652,5 +6655,5 @@
 !--                      Reallocate particle array if necessary. 
                          IF ( new_size > SIZE(particles) )  THEN 
-                            CALL realloc_particles_array(i,j,k,new_size)
+                            CALL realloc_particles_array( i, j, k, new_size )
                          ENDIF
                          old_size  = prt_count(k,j,i)
@@ -6819,5 +6822,5 @@
 !--                      Reallocate particle array if necessary.
                          IF ( new_size > SIZE(particles) )  THEN 
-                            CALL realloc_particles_array(i,j,k,new_size)
+                            CALL realloc_particles_array( i, j, k, new_size )
                          ENDIF
 !
@@ -6829,7 +6832,7 @@
 !
 !--                      Create splitting_factor-1 new particles.
-                         DO jpp = 1, splitting_factor-1
-                            grid_particles(k,j,i)%particles(jpp+old_size) =    &
-                               tmp_particle                  
+                         DO  jpp = 1, splitting_factor-1
+                            grid_particles(k,j,i)%particles( jpp + old_size ) = &
+                               tmp_particle
                          ENDDO
 !
@@ -6840,5 +6843,5 @@
                                                splitting_factor - 1
                       ENDIF
-                   ENDDO  
+                   ENDDO
                 ENDDO
              ENDDO
@@ -6846,5 +6849,5 @@
        ENDDO
     ENDIF
-        
+
     CALL cpu_log( log_point_s(80), 'lpm_splitting', 'stop' )
 
@@ -7610,5 +7613,5 @@
     pack_done     = .FALSE.
 
-    DO n = 1, SIZE(particle_array)
+    DO  n = 1, SIZE(particle_array)
 
        IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
@@ -7633,5 +7636,5 @@
           IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
              IF ( pack_done )  THEN
-                CALL realloc_particles_array (ip,jp,kp)
+                CALL realloc_particles_array ( ip, jp, kp )
              ELSE
                 CALL lpm_pack
@@ -7639,5 +7642,5 @@
                 pindex = prt_count(kp,jp,ip)+1
                 IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
-                   CALL realloc_particles_array (ip,jp,kp)
+                   CALL realloc_particles_array ( ip, jp, kp )
                 ENDIF
                 pack_done = .TRUE.
@@ -7813,5 +7816,5 @@
                    IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )     &
                    THEN
-                      CALL realloc_particles_array(i,j,k)
+                      CALL realloc_particles_array( i, j, k )
                    ENDIF
 
@@ -7855,5 +7858,5 @@
              IF ( number_of_particles <= 0 )  CYCLE
              particles => grid_particles(k,j,i)%particles(1:number_of_particles)         
-             DO n = 1, number_of_particles
+             DO  n = 1, number_of_particles
 !
 !--             Note, check for CFL does not work at first particle timestep 
@@ -7950,14 +7953,14 @@
 
  
-    INTEGER(iwp) ::  i
-    INTEGER(iwp) ::  j
-    INTEGER(iwp) ::  k
-    INTEGER(iwp) :: old_size        !<
-    INTEGER(iwp) :: new_size        !<
-
-    LOGICAL                                        :: dealloc  
-
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
-    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !<
+    INTEGER(iwp) ::  i               !<
+    INTEGER(iwp) ::  j               !<
+    INTEGER(iwp) ::  k               !<
+    INTEGER(iwp) ::  old_size        !<
+    INTEGER(iwp) ::  new_size        !<
+
+    LOGICAL ::  dealloc
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles_d !<
+    TYPE(particle_type), DIMENSION(500)            ::  tmp_particles_s !<
 
     DO  i = nxl, nxr
@@ -8082,7 +8085,7 @@
 !--                the number_of_particles to the actual value.
                    nn = 0
-                   DO is = 0,7
+                   DO  is = 0,7
                       grid_particles(kp,jp,ip)%start_index(is) = nn + 1
-                      DO n = 1,sort_count(is)
+                      DO  n = 1,sort_count(is)
                          nn = nn + 1
                          particles(nn) = sort_particles(n,is)
@@ -8175,16 +8178,16 @@
     SUBROUTINE lpm_sort_timeloop_done
 
-       INTEGER(iwp) :: end_index     !< particle end index for each sub-box
-       INTEGER(iwp) :: i             !< index of particle grid box in x-direction
-       INTEGER(iwp) :: j             !< index of particle grid box in y-direction
-       INTEGER(iwp) :: k             !< index of particle grid box in z-direction
-       INTEGER(iwp) :: n             !< running index for number of particles
-       INTEGER(iwp) :: nb            !< index of subgrid boux
-       INTEGER(iwp) :: nf            !< indices for particles in each sub-box that already finalized their substeps
-       INTEGER(iwp) :: nnf           !< indices for particles in each sub-box that need further treatment
-       INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps
-       INTEGER(iwp) :: start_index   !< particle start index for each sub-box
-
-       TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles  !< temporary particle array 
+       INTEGER(iwp) ::  end_index     !< particle end index for each sub-box
+       INTEGER(iwp) ::  i             !< index of particle grid box in x-direction
+       INTEGER(iwp) ::  j             !< index of particle grid box in y-direction
+       INTEGER(iwp) ::  k             !< index of particle grid box in z-direction
+       INTEGER(iwp) ::  n             !< running index for number of particles
+       INTEGER(iwp) ::  nb            !< index of subgrid boux
+       INTEGER(iwp) ::  nf            !< indices for particles in each sub-box that already finalized their substeps
+       INTEGER(iwp) ::  nnf           !< indices for particles in each sub-box that need further treatment
+       INTEGER(iwp) ::  num_finalized !< number of particles in each sub-box that already finalized their substeps
+       INTEGER(iwp) ::  start_index   !< particle start index for each sub-box
+
+       TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  sort_particles  !< temporary particle array
 
        DO  i = nxl, nxr
