Changeset 4145 for palm


Ignore:
Timestamp:
Aug 6, 2019 9:55:22 AM (5 years ago)
Author:
schwenkel
Message:

Some reformatting

File:
1 edited

Legend:

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

    r4144 r4145  
    2525! -----------------
    2626! $Id$
     27! Some reformatting
     28!
     29! 4144 2019-08-06 09:11:47Z raasch
    2730! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    2831!
     
    387390    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
    388391
    389     INTEGER(iwp), PARAMETER         :: PHASE_INIT    = 1  !<
    390     INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2  !<
     392    INTEGER(iwp), PARAMETER         ::  PHASE_INIT    = 1  !<
     393    INTEGER(iwp), PARAMETER, PUBLIC ::  PHASE_RELEASE = 2  !<
    391394
    392395    SAVE
     
    871874!-- it can be combined as the sgs-velocites for active particles are
    872875!-- calculated differently, i.e. no subboxes are needed.
    873     IF ( .NOT. TRIM(particle_advection_interpolation) == 'trilinear'  .AND.              &
    874        use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
     876    IF ( .NOT. TRIM( particle_advection_interpolation ) == 'trilinear'  .AND.  &
     877       use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    875878          message_string = 'subrgrid scale velocities in combination with ' // &
    876879                           'simple interpolation method is not '            // &
     
    10761079!
    10771080!-- Check which particle interpolation method should be used
    1078     IF ( TRIM(particle_advection_interpolation)  ==  'trilinear' )  THEN
     1081    IF ( TRIM( particle_advection_interpolation )  ==  'trilinear' )  THEN
    10791082       interpolation_simple_corrector = .FALSE.
    10801083       interpolation_simple_predictor = .FALSE.
    10811084       interpolation_trilinear        = .TRUE.
    1082     ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_corrector' )  THEN
     1085    ELSEIF ( TRIM( particle_advection_interpolation )  ==  'simple_corrector' )  THEN
    10831086       interpolation_simple_corrector = .TRUE.
    10841087       interpolation_simple_predictor = .FALSE.
    10851088       interpolation_trilinear        = .FALSE.
    1086     ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_predictor' )  THEN
     1089    ELSEIF ( TRIM( particle_advection_interpolation )  ==  'simple_predictor' )  THEN
    10871090       interpolation_simple_corrector = .FALSE.
    10881091       interpolation_simple_predictor = .TRUE.
     
    13341337!
    13351338!--    Calculate initial_weighting_factor diagnostically
    1336        IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp )  THEN
    1337           initial_weighting_factor =  number_concentration  *                         &
     1339       IF ( number_concentration /= -1.0_wp  .AND. number_concentration > 0.0_wp )  THEN
     1340          initial_weighting_factor =  number_concentration  *                           &
    13381341                                      pdx(1) * pdy(1) * pdz(1)
    13391342       END IF
     
    14711474                            alloc_factor / 100.0_wp ) ), 1 )
    14721475                         IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) )  THEN
    1473                             CALL realloc_particles_array(ip,jp,kp,alloc_size)
     1476                            CALL realloc_particles_array( ip, jp, kp, alloc_size )
    14741477                         ENDIF
    14751478                      ENDIF
     
    15101513!-- Initialize aerosol background spectrum
    15111514    IF ( curvature_solution_effects )  THEN
    1512        CALL lpm_init_aerosols(local_start)
     1515       CALL lpm_init_aerosols( local_start )
    15131516    ENDIF
    15141517!
     
    16291632 SUBROUTINE lpm_init_aerosols(local_start)
    16301633
    1631     REAL(wp)  :: afactor            !< curvature effects
    1632     REAL(wp)  :: bfactor            !< solute effects
    1633     REAL(wp)  :: dlogr              !< logarithmic width of radius bin
    1634     REAL(wp)  :: e_a                !< vapor pressure
    1635     REAL(wp)  :: e_s                !< saturation vapor pressure
    1636     REAL(wp)  :: rmin = 0.005e-6_wp !< minimum aerosol radius
    1637     REAL(wp)  :: rmax = 10.0e-6_wp  !< maximum aerosol radius
    1638     REAL(wp)  :: r_mid              !< mean radius of bin
    1639     REAL(wp)  :: r_l                !< left radius of bin
    1640     REAL(wp)  :: r_r                !< right radius of bin
    1641     REAL(wp)  :: sigma              !< surface tension
    1642     REAL(wp)  :: t_int              !< temperature
     1634    REAL(wp) :: afactor            !< curvature effects
     1635    REAL(wp) :: bfactor            !< solute effects
     1636    REAL(wp) :: dlogr              !< logarithmic width of radius bin
     1637    REAL(wp) :: e_a                !< vapor pressure
     1638    REAL(wp) :: e_s                !< saturation vapor pressure
     1639    REAL(wp) :: rmin = 0.005e-6_wp !< minimum aerosol radius
     1640    REAL(wp) :: rmax = 10.0e-6_wp  !< maximum aerosol radius
     1641    REAL(wp) :: r_mid              !< mean radius of bin
     1642    REAL(wp) :: r_l                !< left radius of bin
     1643    REAL(wp) :: r_r                !< right radius of bin
     1644    REAL(wp) :: sigma              !< surface tension
     1645    REAL(wp) :: t_int              !< temperature
    16431646
    16441647    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
    16451648
    1646     INTEGER(iwp)  :: n              !<
    1647     INTEGER(iwp)  :: ip             !<
    1648     INTEGER(iwp)  :: jp             !<
    1649     INTEGER(iwp)  :: kp             !<
     1649    INTEGER(iwp) :: n              !<
     1650    INTEGER(iwp) :: ip             !<
     1651    INTEGER(iwp) :: jp             !<
     1652    INTEGER(iwp) :: kp             !<
    16501653
    16511654!
    16521655!-- Set constants for different aerosol species
    1653     IF ( TRIM(aero_species) == 'nacl' )  THEN
     1656    IF ( TRIM( aero_species ) == 'nacl' )  THEN
    16541657       molecular_weight_of_solute = 0.05844_wp
    16551658       rho_s                      = 2165.0_wp
    16561659       vanthoff                   = 2.0_wp
    1657     ELSEIF ( TRIM(aero_species) == 'c3h4o4' )  THEN
     1660    ELSEIF ( TRIM( aero_species ) == 'c3h4o4' )  THEN
    16581661       molecular_weight_of_solute = 0.10406_wp
    16591662       rho_s                      = 1600.0_wp
    16601663       vanthoff                   = 1.37_wp
    1661     ELSEIF ( TRIM(aero_species) == 'nh4o3' )  THEN
     1664    ELSEIF ( TRIM( aero_species ) == 'nh4o3' )  THEN
    16621665       molecular_weight_of_solute = 0.08004_wp
    16631666       rho_s                      = 1720.0_wp
     
    16711674!-- The following typical aerosol spectra are taken from Jaenicke (1993):
    16721675!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
    1673     IF ( TRIM(aero_type) == 'polar' )  THEN
     1676    IF ( TRIM( aero_type ) == 'polar' )  THEN
    16741677       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp
    16751678       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp
    16761679       log_sigma = (/ 0.245, 0.300, 0.291 /)
    1677     ELSEIF ( TRIM(aero_type) == 'background' )  THEN
     1680    ELSEIF ( TRIM( aero_type ) == 'background' )  THEN
    16781681       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp
    16791682       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp
    16801683       log_sigma = (/ 0.645, 0.253, 0.425 /)
    1681     ELSEIF ( TRIM(aero_type) == 'maritime' )  THEN
     1684    ELSEIF ( TRIM( aero_type ) == 'maritime' )  THEN
    16821685       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp
    16831686       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp
    16841687       log_sigma = (/ 0.657, 0.210, 0.396 /)
    1685     ELSEIF ( TRIM(aero_type) == 'continental' )  THEN
     1688    ELSEIF ( TRIM( aero_type ) == 'continental' )  THEN
    16861689       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp
    16871690       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp
    16881691       log_sigma = (/ 0.161, 0.217, 0.380 /)
    1689     ELSEIF ( TRIM(aero_type) == 'desert' )  THEN
     1692    ELSEIF ( TRIM( aero_type ) == 'desert' )  THEN
    16901693       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp
    16911694       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp
    16921695       log_sigma = (/ 0.247, 0.770, 0.438 /)
    1693     ELSEIF ( TRIM(aero_type) == 'rural' )  THEN
     1696    ELSEIF ( TRIM( aero_type ) == 'rural' )  THEN
    16941697       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp
    16951698       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp
    16961699       log_sigma = (/ 0.225, 0.557, 0.266 /)
    1697     ELSEIF ( TRIM(aero_type) == 'urban' )  THEN
     1700    ELSEIF ( TRIM( aero_type ) == 'urban' )  THEN
    16981701       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp
    16991702       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp
    17001703       log_sigma = (/ 0.245, 0.666, 0.337 /)
    1701     ELSEIF ( TRIM(aero_type) == 'user' )  THEN
     1704    ELSEIF ( TRIM( aero_type ) == 'user' )  THEN
    17021705       CONTINUE
    17031706    ELSE
     
    21172120!
    21182121!--       If necessary, release new set of particles
    2119           IF ( ( simulated_time - last_particle_release_time ) >= dt_prel  .AND. end_time_prel > simulated_time ) &
    2120           THEN
     2122          IF ( ( simulated_time - last_particle_release_time ) >= dt_prel  .AND.     &
     2123                 end_time_prel > simulated_time )  THEN
    21212124             DO WHILE ( ( simulated_time - last_particle_release_time ) >= dt_prel )
    21222125                CALL lpm_create_particle( PHASE_RELEASE )
     
    22122215!
    22132216!--                   Particle advection
    2214                       CALL lpm_advec(i,j,k)
     2217                      CALL lpm_advec( i, j, k )
    22152218!
    22162219!--                   Particle reflection from walls. Only applied if the particles
    22172220!--                   are in the vertical range of the topography. (Here, some
    22182221!--                   optimization is still possible.)
    2219                       IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
     2222                      IF ( topography /= 'flat'  .AND. k < nzb_max + 2 )  THEN
    22202223                         CALL  lpm_boundary_conds( 'walls', i, j, k )
    22212224                      ENDIF
     
    22232226!--                   User-defined actions after the calculation of the new particle
    22242227!--                   position
    2225                       CALL user_lpm_advec(i,j,k)
     2228                      CALL user_lpm_advec( i, j, k )
    22262229!
    22272230!--                   Apply boundary conditions to those particles that have crossed
     
    22842287!
    22852288!--          IF .FALSE., lpm_sort_and_delete is done inside pcmp
    2286              IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run )   THEN   
     2289             IF ( .NOT. dt_3d_reached  .OR.  .NOT. nested_run )   THEN
    22872290!
    22882291!--             Pack particles (eliminate those marked for deletion),
     
    23862389 SUBROUTINE lpm_write_exchange_statistics
    23872390
    2388     INTEGER(iwp) :: ip         !<
    2389     INTEGER(iwp) :: jp         !<
    2390     INTEGER(iwp) :: kp         !<
    2391     INTEGER(iwp) :: tot_number_of_particles
     2391    INTEGER(iwp) ::  ip         !<
     2392    INTEGER(iwp) ::  jp         !<
     2393    INTEGER(iwp) ::  kp         !<
     2394    INTEGER(iwp) ::  tot_number_of_particles !<
    23922395
    23932396!
     
    24622465    WRITE ( 85 )  simulated_time
    24632466    WRITE ( 85 )  prt_count
    2464          
     2467
    24652468    DO  ip = nxl, nxr
    24662469       DO  jp = nys, nyn
     
    28982901    CHARACTER (LEN=10) ::  version_on_file            !<
    28992902
    2900     INTEGER(iwp) :: alloc_size !<
    2901     INTEGER(iwp) :: ip         !<
    2902     INTEGER(iwp) :: jp         !<
    2903     INTEGER(iwp) :: kp         !<
    2904 
    2905     TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles !<
     2903    INTEGER(iwp) ::  alloc_size !<
     2904    INTEGER(iwp) ::  ip         !<
     2905    INTEGER(iwp) ::  jp         !<
     2906    INTEGER(iwp) ::  kp         !<
     2907
     2908    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles !<
    29062909
    29072910!
     
    30163019    LOGICAL, INTENT(OUT)  ::  found
    30173020
    3018     REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
     3021    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::  tmp_3d   !<
    30193022
    30203023
     
    33593362!--    Predictor step
    33603363       kkw = kp - 1
    3361        DO n = 1, number_of_particles
     3364       DO  n = 1, number_of_particles
    33623365
    33633366          alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
     
    33743377!
    33753378!--    Corrector step
    3376        DO n = 1, number_of_particles
     3379       DO  n = 1, number_of_particles
    33773380
    33783381          IF ( .NOT. particles(n)%particle_mask )  CYCLE
    33793382
    3380           DO nn = 1, iteration_steps
     3383          DO  nn = 1, iteration_steps
    33813384
    33823385!
     
    34283431!--    The particle position for the w velociy is based on the value of kp and kp-1
    34293432       kkw = kp - 1
    3430        DO n = 1, number_of_particles
     3433       DO  n = 1, number_of_particles
    34313434          IF ( .NOT. particles(n)%particle_mask )  CYCLE
    34323435
     
    38303833       ENDDO
    38313834
    3832        DO nb = 0,7
     3835       DO  nb = 0,7
    38333836          i = ip + block_offset(nb)%i_off
    38343837          j = jp + block_offset(nb)%j_off
     
    39633966                ENDIF
    39643967
    3965                 CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),&
     3968                CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),    &
    39663969                                        de_dx_int(n), de_dt, diss_int(n),       &
    39673970                                        dt_particle(n), rg(n,1), term_1_2(n) )
    39683971
    3969                 CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),&
     3972                CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),    &
    39703973                                        de_dy_int(n), de_dt, diss_int(n),       &
    39713974                                        dt_particle(n), rg(n,2), term_1_2(n) )
    39723975
    3973                 CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),&
     3976                CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),    &
    39743977                                        de_dz_int(n), de_dt, diss_int(n),       &
    39753978                                        dt_particle(n), rg(n,3), term_1_2(n) )
     
    40114014                ENDIF
    40124015
    4013                 CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),&
     4016                CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n),    &
    40144017                                        de_dx_int(n), de_dt, diss_int(n),       &
    40154018                                        dt_particle(n), rg(n,1), term_1_2(n) )
    40164019
    4017                 CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),&
     4020                CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n),    &
    40184021                                        de_dy_int(n), de_dt, diss_int(n),       &
    40194022                                        dt_particle(n), rg(n,2), term_1_2(n) )
    40204023
    4021                 CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),&
     4024                CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n),    &
    40224025                                        de_dz_int(n), de_dt, diss_int(n),       &
    40234026                                        dt_particle(n), rg(n,3), term_1_2(n) )
     
    44204423    INTEGER(iwp) ::  tmp_z          !< dummy for sorting algorithm
    44214424
    4422     INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
    4423     INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
    4424     INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
     4425    INTEGER(iwp), DIMENSION(0:10) ::  x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
     4426    INTEGER(iwp), DIMENSION(0:10) ::  y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
     4427    INTEGER(iwp), DIMENSION(0:10) ::  z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
    44254428
    44264429    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
     
    45564559          i2 = particles(n)%x * ddx
    45574560          j2 = particles(n)%y * ddy
    4558           IF (zw(k)   < particles(n)%z ) k2 = k + 1
    4559           IF (zw(k)   > particles(n)%z .AND. zw(k-1) < particles(n)%z ) k2 = k
    4560           IF (zw(k-1) > particles(n)%z ) k2 = k - 1
     4561          IF ( zw(k)   < particles(n)%z ) k2 = k + 1
     4562          IF ( zw(k)   > particles(n)%z  .AND. zw(k-1) < particles(n)%z ) k2 = k
     4563          IF ( zw(k-1) > particles(n)%z ) k2 = k - 1
    45614564!
    45624565!--       Save current particle positions
     
    45904593!--       Wall to the north
    45914594          IF ( prt_y > pos_y_old )  THEN
    4592              ywall = ( j1 +1 ) * dy
     4595             ywall = ( j1 + 1 ) * dy
    45934596!--       Wall to the south
    45944597          ELSE
     
    46454648!--       Store these values only if particle really reaches any wall. t must
    46464649!--       be in a interval between [0:1].
    4647           IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     4650          IF ( t(t_index) <= 1.0_wp  .AND. t(t_index) >= 0.0_wp )  THEN
    46484651             t_index      = t_index + 1
    46494652             cross_wall_x = .TRUE.
     
    46614664          reach_y(t_index) = .TRUE.
    46624665          reach_z(t_index) = .FALSE.
    4663           IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     4666          IF ( t(t_index) <= 1.0_wp  .AND. t(t_index) >= 0.0_wp )  THEN
    46644667             t_index      = t_index + 1
    46654668             cross_wall_y = .TRUE.
     
    46784681          reach_y(t_index) = .FALSE.
    46794682          reach_z(t_index) = .TRUE.
    4680           IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp)  THEN
     4683          IF( t(t_index) <= 1.0_wp  .AND. t(t_index) >= 0.0_wp)  THEN
    46814684             t_index      = t_index + 1
    46824685             cross_wall_z = .TRUE.
     
    46864689!
    46874690!--       Carry out reflection only if particle reaches any wall
    4688           IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z )  THEN
     4691          IF ( cross_wall_x  .OR.  cross_wall_y  .OR. cross_wall_z )  THEN
    46894692!
    46904693!--          Sort fractional timesteps in ascending order. Also sort the
     
    48854888!--          If a particle was reflected, calculate final position from last
    48864889!--          intermediate position.
    4887              IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
     4890             IF ( reflect_x  .OR.  reflect_y  .OR. reflect_z )  THEN
    48884891
    48894892                particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
     
    49224925 SUBROUTINE lpm_droplet_condensation (i,j,k)
    49234926
    4924     INTEGER(iwp), INTENT(IN) :: i              !<
    4925     INTEGER(iwp), INTENT(IN) :: j              !<
    4926     INTEGER(iwp), INTENT(IN) :: k              !<
    4927     INTEGER(iwp) :: n                          !<
     4927    INTEGER(iwp), INTENT(IN) ::  i              !<
     4928    INTEGER(iwp), INTENT(IN) ::  j              !<
     4929    INTEGER(iwp), INTENT(IN) ::  k              !<
     4930    INTEGER(iwp) ::  n                          !<
    49284931
    49294932    REAL(wp) ::  afactor                       !< curvature effects
     
    49544957!
    49554958!-- Parameters for Rosenbrock method (see Verwer et al., 1999)
    4956     REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
    4957     REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
    4958     REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
    4959     REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
     4959    REAL(wp), PARAMETER ::  prec = 1.0E-3_wp     !< precision of Rosenbrock solution
     4960    REAL(wp), PARAMETER ::  q_increase = 1.5_wp  !< increase factor in timestep
     4961    REAL(wp), PARAMETER ::  q_decrease = 0.9_wp  !< decrease factor in timestep
     4962    REAL(wp), PARAMETER ::  gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
    49604963!
    49614964!-- Parameters for terminal velocity
     
    51525155!--    case, the model timestep might be too long.
    51535156       delta_r = new_r(n) - particles(n)%radius
    5154        IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
    5155           WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,             &
     5157       IF ( delta_r < 0.0_wp  .AND.  new_r(n) < 0.0_wp )  THEN
     5158          WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i,                &
    51565159                       ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,               &
    51575160                       '&delta_r=',delta_r,                                    &
     
    61586161 SUBROUTINE turb_enhance_eff
    61596162
    6160     INTEGER(iwp) :: i  !<
    6161     INTEGER(iwp) :: iq !<
    6162     INTEGER(iwp) :: ir !<
    6163     INTEGER(iwp) :: j  !<
    6164     INTEGER(iwp) :: k  !<
    6165     INTEGER(iwp) :: kk !<
     6163    INTEGER(iwp) ::  i  !<
     6164    INTEGER(iwp) ::  iq !<
     6165    INTEGER(iwp) ::  ir !<
     6166    INTEGER(iwp) ::  j  !<
     6167    INTEGER(iwp) ::  k  !<
     6168    INTEGER(iwp) ::  kk !<
    61666169
    61676170    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
     
    63486351    INTEGER(iwp) ::  np               !<
    63496352    INTEGER(iwp) ::  old_size         !< old particle array size
    6350    
     6353
    63516354    INTEGER(iwp), PARAMETER ::  n_max = 100 !< number of radii bin for splitting functions   
    6352    
     6355
    63536356    LOGICAL ::  first_loop_stride_sp = .TRUE. !< flag to calculate constants only once
    63546357
     
    63846387    REAL(wp), DIMENSION(0:n_max-1) ::  r_bin_mid  !< mass weighted mean radius of a bin
    63856388    REAL(wp), DIMENSION(0:n_max)   ::  r_bin      !< boundaries of a radius bin
    6386    
     6389
    63876390    TYPE(particle_type) ::  tmp_particle   !< temporary particle TYPE
    63886391
     
    64426445!--                   Reallocate particle array if necessary.
    64436446                      IF ( new_size > SIZE(particles) )  THEN
    6444                          CALL realloc_particles_array(i,j,k,new_size)
     6447                         CALL realloc_particles_array( i, j, k, new_size )
    64456448                      ENDIF
    64466449                      old_size = prt_count(k,j,i)
     
    66526655!--                      Reallocate particle array if necessary.
    66536656                         IF ( new_size > SIZE(particles) )  THEN
    6654                             CALL realloc_particles_array(i,j,k,new_size)
     6657                            CALL realloc_particles_array( i, j, k, new_size )
    66556658                         ENDIF
    66566659                         old_size  = prt_count(k,j,i)
     
    68196822!--                      Reallocate particle array if necessary.
    68206823                         IF ( new_size > SIZE(particles) )  THEN
    6821                             CALL realloc_particles_array(i,j,k,new_size)
     6824                            CALL realloc_particles_array( i, j, k, new_size )
    68226825                         ENDIF
    68236826!
     
    68296832!
    68306833!--                      Create splitting_factor-1 new particles.
    6831                          DO jpp = 1, splitting_factor-1
    6832                             grid_particles(k,j,i)%particles(jpp+old_size) =    &
    6833                                tmp_particle                 
     6834                         DO  jpp = 1, splitting_factor-1
     6835                            grid_particles(k,j,i)%particles( jpp + old_size ) = &
     6836                               tmp_particle
    68346837                         ENDDO
    68356838!
     
    68406843                                               splitting_factor - 1
    68416844                      ENDIF
    6842                    ENDDO 
     6845                   ENDDO
    68436846                ENDDO
    68446847             ENDDO
     
    68466849       ENDDO
    68476850    ENDIF
    6848        
     6851
    68496852    CALL cpu_log( log_point_s(80), 'lpm_splitting', 'stop' )
    68506853
     
    76107613    pack_done     = .FALSE.
    76117614
    7612     DO n = 1, SIZE(particle_array)
     7615    DO  n = 1, SIZE(particle_array)
    76137616
    76147617       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
     
    76337636          IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
    76347637             IF ( pack_done )  THEN
    7635                 CALL realloc_particles_array (ip,jp,kp)
     7638                CALL realloc_particles_array ( ip, jp, kp )
    76367639             ELSE
    76377640                CALL lpm_pack
     
    76397642                pindex = prt_count(kp,jp,ip)+1
    76407643                IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
    7641                    CALL realloc_particles_array (ip,jp,kp)
     7644                   CALL realloc_particles_array ( ip, jp, kp )
    76427645                ENDIF
    76437646                pack_done = .TRUE.
     
    78137816                   IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )     &
    78147817                   THEN
    7815                       CALL realloc_particles_array(i,j,k)
     7818                      CALL realloc_particles_array( i, j, k )
    78167819                   ENDIF
    78177820
     
    78557858             IF ( number_of_particles <= 0 )  CYCLE
    78567859             particles => grid_particles(k,j,i)%particles(1:number_of_particles)         
    7857              DO n = 1, number_of_particles
     7860             DO  n = 1, number_of_particles
    78587861!
    78597862!--             Note, check for CFL does not work at first particle timestep
     
    79507953
    79517954 
    7952     INTEGER(iwp) ::  i
    7953     INTEGER(iwp) ::  j
    7954     INTEGER(iwp) ::  k
    7955     INTEGER(iwp) :: old_size        !<
    7956     INTEGER(iwp) :: new_size        !<
    7957 
    7958     LOGICAL                                        :: dealloc 
    7959 
    7960     TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
    7961     TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !<
     7955    INTEGER(iwp) ::  i               !<
     7956    INTEGER(iwp) ::  j               !<
     7957    INTEGER(iwp) ::  k               !<
     7958    INTEGER(iwp) ::  old_size        !<
     7959    INTEGER(iwp) ::  new_size        !<
     7960
     7961    LOGICAL ::  dealloc
     7962
     7963    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles_d !<
     7964    TYPE(particle_type), DIMENSION(500)            ::  tmp_particles_s !<
    79627965
    79637966    DO  i = nxl, nxr
     
    80828085!--                the number_of_particles to the actual value.
    80838086                   nn = 0
    8084                    DO is = 0,7
     8087                   DO  is = 0,7
    80858088                      grid_particles(kp,jp,ip)%start_index(is) = nn + 1
    8086                       DO n = 1,sort_count(is)
     8089                      DO  n = 1,sort_count(is)
    80878090                         nn = nn + 1
    80888091                         particles(nn) = sort_particles(n,is)
     
    81758178    SUBROUTINE lpm_sort_timeloop_done
    81768179
    8177        INTEGER(iwp) :: end_index     !< particle end index for each sub-box
    8178        INTEGER(iwp) :: i             !< index of particle grid box in x-direction
    8179        INTEGER(iwp) :: j             !< index of particle grid box in y-direction
    8180        INTEGER(iwp) :: k             !< index of particle grid box in z-direction
    8181        INTEGER(iwp) :: n             !< running index for number of particles
    8182        INTEGER(iwp) :: nb            !< index of subgrid boux
    8183        INTEGER(iwp) :: nf            !< indices for particles in each sub-box that already finalized their substeps
    8184        INTEGER(iwp) :: nnf           !< indices for particles in each sub-box that need further treatment
    8185        INTEGER(iwp) :: num_finalized !< number of particles in each sub-box that already finalized their substeps
    8186        INTEGER(iwp) :: start_index   !< particle start index for each sub-box
    8187 
    8188        TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: sort_particles  !< temporary particle array
     8180       INTEGER(iwp) ::  end_index     !< particle end index for each sub-box
     8181       INTEGER(iwp) ::  i             !< index of particle grid box in x-direction
     8182       INTEGER(iwp) ::  j             !< index of particle grid box in y-direction
     8183       INTEGER(iwp) ::  k             !< index of particle grid box in z-direction
     8184       INTEGER(iwp) ::  n             !< running index for number of particles
     8185       INTEGER(iwp) ::  nb            !< index of subgrid boux
     8186       INTEGER(iwp) ::  nf            !< indices for particles in each sub-box that already finalized their substeps
     8187       INTEGER(iwp) ::  nnf           !< indices for particles in each sub-box that need further treatment
     8188       INTEGER(iwp) ::  num_finalized !< number of particles in each sub-box that already finalized their substeps
     8189       INTEGER(iwp) ::  start_index   !< particle start index for each sub-box
     8190
     8191       TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  sort_particles  !< temporary particle array
    81898192
    81908193       DO  i = nxl, nxr
Note: See TracChangeset for help on using the changeset viewer.