particles.f90

Go to the documentation of this file.

!********************************************************************************
!
!********************************************************************************
MODULE particles_module
  !
  USE moments_module, ONLY : n_mom , n_nodes , n_sections

  USE variables, ONLY : aggregation_flag , verbose_level , indent_space , fmt

  USE variables, ONLY : k_b
  USE variables, ONLY : wp

  IMPLICIT NONE
  
  INTEGER :: n_part

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_partial_mass_fraction

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_partial_mass_fraction0

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_partial_volume_fraction

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_mass_fraction

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_mass_fraction0

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: solid_volume_fraction

  REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: bin_partial_mass_fraction

  REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: particle_loss_rate
  
  REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cum_particle_loss_rate

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: mom

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: mom0
  
  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: set_mom

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: set_cp_mom

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: birth_mom

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: death_mom

  REAL(wp), ALLOCATABLE :: shape_factor(:)

  REAL(wp), ALLOCATABLE :: phi1(:)

  REAL(wp), ALLOCATABLE :: rho1(:)

  REAL(wp), ALLOCATABLE :: shape1(:)
  
  REAL(wp), ALLOCATABLE :: phi2(:)

  REAL(wp), ALLOCATABLE :: rho2(:)

  REAL(wp), ALLOCATABLE :: shape2(:)

  REAL(wp), ALLOCATABLE :: cp_part(:)

  REAL(wp) :: cpsolid

  
  REAL(wp) :: particles_beta0 

  CHARACTER*10 :: settling_model

  CHARACTER(LEN=20) :: distribution

  LOGICAL, ALLOCATABLE :: aggregation_array(:)

  REAL(wp), ALLOCATABLE :: aggregate_porosity(:)
  
  CHARACTER(LEN=20) :: aggregation_model

  INTEGER, ALLOCATABLE :: aggr_idx(:)

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: m_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: w_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: phi_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: diam_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: vol_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: rho_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: shape_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: set_vel_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: cp_quad

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: f_quad
  
  REAL(wp) :: t_part

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: phil 

  REAL(wp), ALLOCATABLE, DIMENSION(:) :: phir 

  REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: m

  LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: q_flag

  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:,:) :: kernel_aggr
  
  REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:,:,:,:) :: a

  REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:,:,:,:) :: wij

  SAVE

CONTAINS

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE allocate_particles

    use, intrinsic :: iso_fortran_env
    use, intrinsic :: ieee_arithmetic

    IMPLICIT NONE

    REAL(wp) :: notSet

    notset = ieee_value(0.0_wp, ieee_quiet_nan)

    ALLOCATE ( solid_partial_mass_fraction(1:n_part) )
    ALLOCATE ( solid_partial_mass_fraction0(1:n_part) )
    ALLOCATE ( solid_partial_volume_fraction(1:n_part) )
    ALLOCATE ( solid_mass_fraction(1:n_part) )
    ALLOCATE ( solid_mass_fraction0(1:n_part) )
    ALLOCATE ( solid_volume_fraction(1:n_part) )

    ALLOCATE ( bin_partial_mass_fraction(1:n_sections,1:n_part) )

    ALLOCATE ( particle_loss_rate(1:n_part,1:n_sections) )
    ALLOCATE ( cum_particle_loss_rate(1:n_part,1:n_sections) )
    
    ! Allocation of the arrays for the moments
    ALLOCATE ( mom0(0:n_mom-1,1:n_sections,1:n_part) )
    ALLOCATE ( mom(0:n_mom-1,1:n_sections,1:n_part) )
    ALLOCATE ( set_mom(0:n_mom-1,1:n_sections,1:n_part) )
    ALLOCATE ( set_cp_mom(0:n_mom-1,1:n_sections,1:n_part) )
    ALLOCATE ( birth_mom(0:n_mom-1,1:n_sections,1:n_part) )
    ALLOCATE ( death_mom(0:n_mom-1,1:n_sections,1:n_part) )

    ! Allocation of the parameters for the variable density
    ALLOCATE ( shape_factor(n_part) )
    ALLOCATE ( phi1(n_part) )
    ALLOCATE ( rho1(n_part) )
    ALLOCATE ( shape1(n_part) )
    ALLOCATE ( phi2(n_part) )
    ALLOCATE ( rho2(n_part) )
    ALLOCATE ( shape2(n_part) )

    ALLOCATE ( cp_part(n_part) )

    phi1(1:n_part) = notset
    phi2(1:n_part) = notset
    
    shape_factor(1:n_part) = notset
    shape1(1:n_part) = notset
    shape2(1:n_part) = notset

    !Allocation of arrays for quadrature variables
    ALLOCATE ( m_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( w_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( f_quad(n_nodes,n_sections,n_part) )

    ALLOCATE ( phi_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( diam_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( vol_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( rho_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( shape_quad(n_nodes,n_sections,n_part) )
    ALLOCATE ( set_vel_quad(n_nodes,n_sections,n_part) )
    
    ! Allocation of arrays for aggregation
    ALLOCATE ( aggregation_array(n_part) )
    ALLOCATE ( aggregate_porosity(n_part) )
    ALLOCATE ( aggr_idx(n_part) )

    ALLOCATE ( cp_quad(n_nodes,n_sections,n_part) )

    ALLOCATE ( phil(1:n_sections) , phir(1:n_sections) )

    ALLOCATE ( m(1:n_sections+1,1:n_part) )

    ALLOCATE( q_flag(n_sections,n_sections,n_sections,1:n_part,1:n_part) )

    ALLOCATE( kernel_aggr(n_nodes,n_nodes,n_sections,n_part,n_sections,n_part) )

    ALLOCATE( a(n_nodes,n_nodes,n_sections,n_sections,n_part,n_sections,n_part) )

    ALLOCATE( wij(n_nodes,n_nodes,0:n_mom-1,n_sections,n_part,n_sections,       &
         n_part) )
  

  END SUBROUTINE allocate_particles

  SUBROUTINE deallocate_particles

    IMPLICIT NONE

    DEALLOCATE ( solid_partial_mass_fraction )
    DEALLOCATE ( solid_partial_mass_fraction0 )
    DEALLOCATE ( solid_partial_volume_fraction )
    DEALLOCATE ( solid_mass_fraction )
    DEALLOCATE ( solid_mass_fraction0 )
    DEALLOCATE ( solid_volume_fraction )

    DEALLOCATE ( bin_partial_mass_fraction )

    DEALLOCATE ( particle_loss_rate , cum_particle_loss_rate )
    
    ! Allocation of the arrays for the moments
    DEALLOCATE ( mom0 )
    DEALLOCATE ( mom )
    DEALLOCATE ( set_mom )
    DEALLOCATE ( set_cp_mom )
    DEALLOCATE ( birth_mom )
    DEALLOCATE ( death_mom )

    ! Allocation of the parameters for the variable density
    DEALLOCATE ( shape_factor )
    DEALLOCATE ( phi1 )
    DEALLOCATE ( rho1 )
    DEALLOCATE ( shape1 )
    DEALLOCATE ( phi2 )
    DEALLOCATE ( rho2 )
    DEALLOCATE ( shape2 )

    DEALLOCATE ( cp_part )

    DEALLOCATE ( m_quad )
    DEALLOCATE ( w_quad )
    DEALLOCATE ( f_quad )
    DEALLOCATE ( cp_quad )

    DEALLOCATE ( set_vel_quad )

    ! Allocation of arrays for aggregation
    DEALLOCATE ( aggregation_array )
    DEALLOCATE ( aggregate_porosity )
    DEALLOCATE ( aggr_idx )

    DEALLOCATE ( phil , phir )

    DEALLOCATE ( m)

    DEALLOCATE ( q_flag )
    DEALLOCATE ( kernel_aggr )

    DEALLOCATE ( a )
    DEALLOCATE ( wij )
    
  END SUBROUTINE deallocate_particles

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION particles_settling_velocity(diam,rhop,shape_fact)
    !
    USE meteo_module, ONLY : rho_atm , rho_atm0 , visc_atm

    USE variables, ONLY : gi , pi_g , verbose_level

    IMPLICIT NONE

    REAL(wp) :: particles_settling_velocity

    REAL(wp), INTENT(IN) :: diam
    REAL(wp), INTENT(IN) :: rhop
    REAL(wp), INTENT(IN) :: shape_fact

    REAL(wp) :: k1 , k2 , k3
    REAL(wp) :: CD , CD1 , CD2

    REAL(wp) :: Reynolds , Reynoldsk1k2
    REAL(wp) :: Vinit , Vg_Ganser

    INTEGER :: i

    REAL(wp) :: A_cs

    REAL(wp) :: Cd_100 , Cd_1000

    REAL(wp) :: Cd_interp

    REAL(wp) :: Us_100 , Us_1000

    REAL(wp) :: mass

    REAL(wp) :: Us , Us_1 ,Us_2

    REAL(wp) :: Rey1 , Rey2

    REAL(wp) :: c0 , c1 , c2

    REAL(wp) :: sqrt_delta

    IF ( settling_model .EQ. 'textor' ) THEN

       ! Textor et al. 2006

       IF ( diam .LE. 1.e-4_wp ) THEN

          k1 = 1.19d5   ! (m^2 kg^-1 s^-1 )
          
          particles_settling_velocity = k1 * rhop * sqrt( rho_atm0 / rho_atm ) &
               * ( 0.5_wp * diam )**2

       ELSEIF ( diam .LE. 1.e-3_wp ) THEN

          k2 = 8.0_wp    ! (m^3 kg^-1 s^-1 )

          particles_settling_velocity = k2 * rhop * sqrt( rho_atm0 / rho_atm ) &
               * ( 0.5_wp * diam )

       ELSE 

          k3 = 4.833_wp ! (m^2 kg^-0.5 s^-1 )
          cd = 0.75_wp

          particles_settling_velocity = k3 * sqrt( rhop / cd )                 &
               * sqrt(  rho_atm0 / rho_atm ) * sqrt( 0.5_wp * diam )

       END IF

    ELSEIF ( settling_model .EQ. 'ganser' ) THEN 

       vinit = diam**2 * gi * ( rhop - rho_atm ) / (18.0_wp*visc_atm)

       DO i=1,10

          IF (i.EQ.1) reynolds = rho_atm * vinit * diam / visc_atm

          k1 = 3.0_wp / (1.0_wp + 2.0_wp * ( shape_fact**(-0.5_wp)))

          k2 = 10.0_wp**(1.8148_wp * ((-1.0_wp*log10(shape_fact))**0.5743_wp))

          reynoldsk1k2 = reynolds * k1 * k2

          cd1 = k2 * 24.0_wp / reynoldsk1k2  *                                     &
               ( 1.0_wp + 0.1118_wp * reynoldsk1k2**0.6567_wp )

          cd2 = 0.4305_wp * k2 / ( 1.0_wp + 3305.0_wp / reynoldsk1k2 )

          cd = cd1 + cd2

          vg_ganser = ( ( 4.0_wp * gi * diam * ( rhop - rho_atm ) )                &
               / ( 3.0_wp * cd * rho_atm) )**0.5_wp

          reynolds = rho_atm * vg_ganser * diam / visc_atm

       ENDDO

       particles_settling_velocity = vg_ganser

       IF ( vg_ganser .LE. 0.0_wp ) THEN

          WRITE(*,*) 'diam',diam
          WRITE(*,*) 'NEGATIVE VALUE', vinit,vg_ganser
          READ(*,*)
          
       END IF

    ELSEIF ( settling_model .EQ. 'pfeiffer' ) THEN

       k1 = shape_fact**(-0.828_wp)
       k2 = 2.0_wp * sqrt( 1.07_wp - shape_fact )

       mass = rhop * 4.0_wp/3.0_wp * pi_g * ( 0.5_wp * diam )**3

       a_cs = pi_g * ( 0.5_wp * diam )**2

       c0 = -2.0_wp * diam * mass * gi
       c1 = 24.0_wp * visc_atm * k1 * a_cs
       c2 = rho_atm * diam * k2 * a_cs

       sqrt_delta = sqrt( c1**2 - 4.0_wp * c0*c2 )

       us_1 = ( - c1 + sqrt_delta ) / ( 2.0_wp * c2 )
       us_2 = ( - c1 - sqrt_delta ) / ( 2.0_wp * c2 )


       cd_100 = 24.0_wp/100.0_wp * k1 + k2
       us_100 = sqrt( 2.0_wp * mass * gi / ( cd_100*rho_atm * a_cs ) )

       cd_1000 = 1.0_wp
       us_1000 = sqrt( 2.0_wp * mass * gi / ( cd_1000*rho_atm * a_cs ) )

       rey1 = rho_atm * diam * us_1 / visc_atm
       rey2 = rho_atm * diam * us_2 / visc_atm

       IF ( verbose_level .GE. 4 ) THEN

          WRITE(*,*) 'rho_atm , diam , Us_1 , visc_atm',rho_atm , diam , us_1 , &
               visc_atm
          WRITE(*,*) 'Rey1,Rey2',rey1,rey2
          READ(*,*)

       END IF

       ! Initialization only
       us = us_1000

       IF ( ( rey1 .GT. 0.0_wp ) .AND. ( rey1 .LE. 100.0_wp ) ) THEN

          ! For small Reynolds numbers the drag coefficient is given by Eq.8
          ! of Pfeiffer et al. 2005 and the settling velocity is Us_1

          us = us_1  

       ELSEIF ( ( rey1 .GT. 100.0_wp ) .AND. ( rey1 .LE. 1000.0_wp ) ) THEN

          ! For intermediate Reyonlds numbers, 100<Re<1000, the drag coefficient 
          ! is linearly interpolated between Cd_100 and Cd_1000

          cd_interp = cd_100 + ( rey1 - 100.0_wp ) / ( 1000.0_wp - 100.0_wp ) *                &
               ( cd_1000 - cd_100)
          us = sqrt( 2.0_wp * mass * gi / ( cd_interp * rho_atm * a_cs ) )

       ELSEIF ( rey1 .GT. 1000.0_wp ) THEN

          ! For large Reynolds numbers the drag coefficient is taken as Cd=1,
          ! as in Pfeiffer et al. 2005 with the settling velocity is Us_1000

          us = us_1000

       END IF

       IF ( ( rey2 .GT. 0.0_wp ) .AND. ( rey2 .LE. 100.0_wp ) ) THEN 

          us = us_2

       ELSEIF ( ( rey2 .GT. 100.0_wp ) .AND. ( rey2 .LE. 1000.0_wp ) ) THEN 

          cd_interp = cd_100 + ( rey2 - 100.0_wp ) / ( 1000.0_wp - 100.0_wp )                  &
               * ( cd_1000 - cd_100 )

          us = sqrt( 2.0_wp * mass * gi / ( cd_interp * rho_atm * a_cs ) )

       ELSEIF ( rey2 .GT. 1000.0_wp ) THEN

          us = us_1000

       END IF

       particles_settling_velocity = us

    ELSE

       WRITE(*,*) 'wrong settling model'
       stop

    END IF

    !WRITE(*,*) 'diam',diam
    !WRITE(*,*) 'rhop',rhop
    !WRITE(*,*) 'shape factor',shape_fact
    !WRITE(*,*) 'particles_settling_velocit',particles_settling_velocity
    !READ(*,*)

    RETURN

  END FUNCTION particles_settling_velocity

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION particles_heat_capacity(i_part,diam)
    !
    IMPLICIT NONE

    REAL(wp) :: particles_heat_capacity
    INTEGER, INTENT(IN) :: i_part
    REAL(wp), INTENT(IN) :: diam

    particles_heat_capacity = cp_part(i_part)

  END FUNCTION particles_heat_capacity

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION particles_density(i_part,phi)
    !
    IMPLICIT NONE

    REAL(wp) :: particles_density

    INTEGER, INTENT(IN) :: i_part
    REAL(wp), INTENT(IN) :: phi

    REAL(wp) :: phi_temp,rho_temp
    
    IF ( phi1(i_part) .GT. phi2(i_part) ) THEN

       phi_temp = phi1(i_part)
       rho_temp = rho1(i_part)

       phi1(i_part) = phi2(i_part)
       rho1(i_part) = rho2(i_part)

       phi2(i_part) = phi_temp
       rho2(i_part) = rho_temp

    END IF
       
    IF ( phi .LE. phi1(i_part) ) THEN

       particles_density = rho1(i_part)

    ELSEIF ( phi .LE. phi2(i_part) ) THEN


       particles_density = rho1(i_part) + ( phi - phi1(i_part) ) /              &
            ( phi2(i_part) - phi1(i_part) ) * ( rho2(i_part) - rho1(i_part) )

    ELSE

       particles_density = rho2(i_part)

    END IF
    
    RETURN

  END FUNCTION particles_density

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION particles_shape(i_part,phi)
    !
    IMPLICIT NONE

    REAL(wp) :: particles_shape

    INTEGER, INTENT(IN) :: i_part
    REAL(wp), INTENT(IN) :: phi

    REAL(wp) :: phi_temp,shape_temp
    
    IF ( phi1(i_part) .GT. phi2(i_part) ) THEN

       phi_temp = phi1(i_part)
       shape_temp = shape1(i_part)

       phi1(i_part) = phi2(i_part)
       shape1(i_part) = shape2(i_part)

       phi2(i_part) = phi_temp
       shape2(i_part) = shape_temp

    END IF
       
    IF ( phi .LE. phi1(i_part) ) THEN

       particles_shape = shape1(i_part)

    ELSEIF ( phi .LE. phi2(i_part) ) THEN


       particles_shape = shape1(i_part) + ( phi - phi1(i_part) ) /              &
            ( phi2(i_part) - phi1(i_part) ) * ( shape2(i_part) - shape1(i_part) )

    ELSE

       particles_shape = shape2(i_part)

    END IF
    
    RETURN

  END FUNCTION particles_shape

  
  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION particles_beta(temp,visc,diam_i,diam_j,rho_i,rho_j,Vs_i,Vs_j,lw_vf,  &
       ice_vf,solid_mf)

    IMPLICIT NONE

    REAL(wp) :: particles_beta

    REAL(wp), INTENT(IN) :: temp
    REAL(wp), INTENT(IN) :: visc
    REAL(wp), INTENT(IN) :: diam_i
    REAL(wp), INTENT(IN) :: diam_j
    REAL(wp), INTENT(IN), OPTIONAL :: rho_i
    REAL(wp), INTENT(IN), OPTIONAL :: rho_j
    REAL(wp), INTENT(IN), OPTIONAL :: Vs_i
    REAL(wp), INTENT(IN), OPTIONAL :: Vs_j
    REAL(wp), INTENT(IN), OPTIONAL :: lw_vf 
    REAL(wp), INTENT(IN), OPTIONAL :: ice_vf 
    REAL(wp), INTENT(IN), OPTIONAL :: solid_mf 

    SELECT CASE ( aggregation_model )
       
    CASE DEFAULT

       particles_beta = 0.0_wp

    CASE ( 'constant' )

       particles_beta = particles_beta0

    CASE ( 'brownian' )

       particles_beta = ( 2.0_wp * k_b * temp ) / ( 3.0_wp * visc ) *               &
            ( diam_i + diam_j ) ** 2 / ( diam_i * diam_j ) 

    CASE ( 'sum' )

       particles_beta = diam_i**3 + diam_j**3

    CASE ( 'costa')

       IF ( max(lw_vf,ice_vf) .LT. 1.0e-10_wp ) THEN

          particles_beta = 0.0_wp

       ELSE
       
          particles_beta = aggregation_kernel(diam_i,rho_i,vs_i,diam_j,rho_j,vs_j, &
               lw_vf,ice_vf,solid_mf,temp,visc)

       END IF
          
    END SELECT

    IF ( verbose_level .GE. 2 ) THEN

       WRITE(*,*) 'beta =',particles_beta
       WRITE(*,*) 

       WRITE(*,fmt) ' ','END particles_beta'
       indent_space = indent_space - 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"

    END IF

    RETURN

  END FUNCTION particles_beta

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION aggregation_kernel( diam_i , rho_i , Vs_i , diam_j , rho_j , Vs_j ,  &
       lw_vf , ice_vf , solid_mf , temp , visc )

    IMPLICIT NONE

    REAL(wp) :: aggregation_kernel

    REAL(wp), INTENT(IN) :: diam_i
    REAL(wp), INTENT(IN) :: rho_i
    REAL(wp), INTENT(IN) :: Vs_i
    REAL(wp), INTENT(IN) :: diam_j
    REAL(wp), INTENT(IN) :: rho_j
    REAL(wp), INTENT(IN) :: Vs_j
    REAL(wp), INTENT(IN) :: lw_vf 
    REAL(wp), INTENT(IN) :: ice_vf 
    REAL(wp), INTENT(IN) :: solid_mf 
    REAL(wp), INTENT(IN) :: temp
    REAL(wp), INTENT(IN) :: visc

    REAL(wp) :: beta
    REAL(wp) :: alfa

    IF ( verbose_level .GE. 2 ) THEN

       indent_space = indent_space + 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"
       WRITE(*,fmt) ' ','BEGINNING aggregation_kernel'
       READ(*,*)

    END IF

    !WRITE(*,*) 'aggregation_kernel: Vs_i,Vs_j',Vs_i,Vs_j

    beta = collision_kernel(diam_i,rho_i,vs_i,diam_j,rho_j,vs_j,temp,visc)

    alfa = coalescence_efficiency(diam_i,rho_i,vs_i,diam_j,rho_j,vs_j,lw_vf,    &
         ice_vf,solid_mf)

    aggregation_kernel = beta * alfa

    IF ( aggregation_kernel .GT. 0.0_wp ) THEN

       !WRITE(*,*) 'lw_vf,ice_vf',lw_vf,ice_vf
       !WRITE(*,*) 'aggregation_kernel, beta, alfa',aggregation_kernel, beta, alfa
       !READ(*,*)

    END IF
       
    IF ( verbose_level .GE. 2 ) THEN

       WRITE(*,fmt) ' ','END aggregation_kernel'
       indent_space = indent_space - 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"

    END IF

    RETURN

  END FUNCTION aggregation_kernel

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION collision_kernel(diam_i,rho_i,Vs_i,diam_j,rho_j,Vs_j,temp,visc)

    USE meteo_module, ONLY : visc_atm

    USE variables, ONLY: pi_g

    IMPLICIT NONE

    REAL(wp) :: collision_kernel

    REAL(wp),INTENT(IN) :: diam_i
    REAL(wp),INTENT(IN) :: rho_i
    REAL(wp),INTENT(IN) :: Vs_i
    REAL(wp),INTENT(IN) :: diam_j
    REAL(wp),INTENT(IN) :: rho_j
    REAL(wp),INTENT(IN) :: Vs_j
    REAL(wp),INTENT(IN) :: temp
    REAL(wp),INTENT(IN) :: visc

    REAL(wp) :: beta_B   

    REAL(wp) :: beta_S

    REAL(wp) :: beta_DS

    REAL(wp) :: E_coll

    REAL(wp) :: epsilon

    REAL(wp) :: Gamma_s

    !WRITE(*,*) 'collision_kernel'

    IF ( verbose_level .GE. 2 ) THEN

       indent_space = indent_space + 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"
       WRITE(*,fmt) ' ','BEGINNING collision_kernel'
       READ(*,*)

    END IF

    ! Eq. 3, first term Costa et al. JGR 2010
    beta_b = 2.0_wp / 3.0_wp * k_b * temp / visc * ( diam_i + diam_j )**2           &
         / ( diam_i*diam_j ) 

    ! Value from Table 1 (Costa 2010)
    gamma_s = 0.0045_wp 

    ! Eq. 3, second term Costa et al. JGR 2010
    beta_s = 1.0_wp / 6.0_wp * gamma_s * ( diam_i + diam_j )**3

    ! Eq. 3, third term Costa et al. JGR 2010
    beta_ds = pi_g / 4.0_wp * ( diam_i + diam_j )**2 * abs( vs_j - vs_i )

    ! WRITE(*,*) 'beta_B,beta_S,beta_DS',beta_B,beta_S, beta_DS

    collision_kernel = beta_b + beta_s + beta_ds

    IF ( verbose_level .GE. 2 ) THEN

       WRITE(*,fmt) ' ','END collision_kernel'
       indent_space = indent_space - 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"

    END IF

    RETURN

  END FUNCTION collision_kernel

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION coalescence_efficiency(diam_i,rho_i,Vs_i,diam_j,rho_j,Vs_j,lw_vf,    &
       ice_vf,solid_mf)

    USE variables, ONLY: gi

    IMPLICIT NONE

    REAL(wp) :: coalescence_efficiency

    REAL(wp), INTENT(IN) :: diam_i
    REAL(wp), INTENT(IN) :: rho_i
    REAL(wp), INTENT(IN) :: Vs_i
    REAL(wp), INTENT(IN) :: diam_j
    REAL(wp), INTENT(IN) :: rho_j
    REAL(wp), INTENT(IN) :: Vs_j
    REAL(wp), INTENT(IN) :: lw_vf 
    REAL(wp), INTENT(IN) :: ice_vf 
    REAL(wp), INTENT(IN) :: solid_mf 
    
    REAL(wp) :: coalescence_efficiency_ice , coalescence_efficiency_water
    
    REAL(wp) :: Stokes

    REAL(wp) :: Stokes_cr

    REAL(wp) :: q

    REAL(wp) :: mu_liq

    IF ( verbose_level .GE. 2 ) THEN

       indent_space = indent_space + 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"
       WRITE(*,fmt) ' ','BEGINNING coalescence_efficiency'

    END IF

    ! Modified from Eq. 5 Costa et al. JGR 2010
    coalescence_efficiency_ice = 0.09_wp * ice_vf
            
    mu_liq = 5.43e-4_wp
    
    ! Eq. 6 Costa et al. JGR 2010 (CHECK DENSITY!)
    stokes = 8.0_wp * ( 0.5_wp * ( rho_i + rho_j ) ) * abs( vs_i - vs_j ) /     &
         ( 9.0_wp * mu_liq ) * diam_i * diam_j / ( diam_i + diam_j )
    
    stokes_cr = 1.3_wp
    
    q = 0.8_wp
    
    ! Modified from Eq. 8 Costa et al. JGR 2010
    coalescence_efficiency_water = 1.0_wp / ( 1.0_wp+( stokes/stokes_cr ) )**q  &
         * lw_vf

    IF ( lw_vf .GT. 0.0_wp ) THEN

       IF ( ice_vf .GT. 0.0_wp ) THEN

          coalescence_efficiency = coalescence_efficiency_water                 &
               + coalescence_efficiency_ice

       ELSE
       
          coalescence_efficiency = coalescence_efficiency_water

       END IF
          
    ELSEIF ( ice_vf .GT. 0.0_wp ) THEN

       coalescence_efficiency = coalescence_efficiency_ice
          
    ELSE

       coalescence_efficiency = 0.0_wp

    END IF
          
    
    IF ( verbose_level .GE. 2 ) THEN

       WRITE(*,fmt) ' ','END coalescence_efficiency'
       indent_space = indent_space - 2
       WRITE(fmt,*) indent_space
       fmt = "(A" // trim(fmt) // ",A)"

    END IF

    RETURN

  END FUNCTION coalescence_efficiency


  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE eval_particles_moments

    ! external variables
    USE variables, ONLY : verbose_level

    IMPLICIT NONE

    INTEGER :: i_part , j_part
    INTEGER :: i_sect , j_sect , k_sect
    
    INTEGER :: i_mom
    INTEGER :: i_node , j_node

    DO i_part=1,n_part
       
       DO i_sect=1,n_sections

          DO i_node=1,n_nodes

             ! Settling velocities at the quadrature points
             set_vel_quad(i_node,i_sect,i_part) =                               &
                  particles_settling_velocity( diam_quad(i_node,i_sect,i_part), &
                  rho_quad(i_node,i_sect,i_part) ,                              &
                  shape_quad(i_node,i_sect,i_part) )

          END DO
          
          DO i_mom=0,n_mom-1

             IF ( mom(1,i_sect,i_part) .GT. 0.0_wp ) THEN
             
                set_mom(i_mom,i_sect,i_part) = sum( set_vel_quad(:,i_sect,i_part)  &
                     * f_quad(:,i_sect,i_part) * w_quad(:,i_sect,i_part)           &
                     * m_quad(:,i_sect,i_part)**i_mom ) / mom(i_mom,i_sect,i_part)
                
                set_cp_mom(i_mom,i_sect,i_part) =                                  &
                     sum( set_vel_quad(:,i_sect,i_part )                           &
                     * cp_quad(:,i_sect,i_part) * f_quad(:,i_sect,i_part)          &
                     * w_quad(:,i_sect,i_part) * m_quad(:,i_sect,i_part)**i_mom )  &
                     / mom(i_mom,i_sect,i_part) 
                
             ELSE

                set_mom(i_mom,i_sect,i_part) = 0.0_wp
                
                set_cp_mom(i_mom,i_sect,i_part) = 0.0_wp

             END IF
             
             IF ( verbose_level .GE. 2 ) THEN
                
                WRITE(*,*) 'i_part,i_mom',i_part,i_mom
                WRITE(*,*) 'abscissas', m_quad(1:n_nodes,i_sect,i_part)
                WRITE(*,*) 'set_mom(i_mom,i_sect,i_part) = ' ,                  &
                     set_mom(i_mom,i_sect,i_part)
                
             END IF
             
          END DO
          
       END DO

    END DO
       
    RETURN

  END SUBROUTINE eval_particles_moments

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE eval_quad_values

    USE moments_module, ONLY : linear_reconstruction
    
    IMPLICIT NONE

    INTEGER :: i_part , i_sect
    REAL(wp) :: coeff_lin(4)
    REAL(wp) :: Mai , Mbi
    REAL(wp) :: alfai , betai,gamma1,gamma2

    REAL(wp) :: Ml , Mr

    INTEGER :: condition1(n_nodes) , condition2(n_nodes)
       
    DO i_part=1,n_part
       
       DO i_sect = 1,n_sections

          ml = m(i_sect,i_part)
          mr = m(i_sect+1,i_part)
          
          CALL linear_reconstruction( ml,mr,mom(:,i_sect,i_part), mai , mbi ,   &
               alfai , betai , gamma1 , gamma2 )

          f_quad(:,i_sect,i_part) = alfai * ( ( mr - m_quad(:,i_sect,i_part) )  &
               / ( mr - ml ) )**gamma1 + ( betai - alfai )                      &
               * ( ( m_quad(:,i_sect,i_part) - ml ) / ( mr - ml ) )**gamma2
          
       END DO
       
    END DO

    RETURN

  END SUBROUTINE eval_quad_values

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE init_quadrature_points

    USE moments_module, ONLY : gaulegf
    
    IMPLICIT NONE

    REAL(wp) :: x(n_nodes) , w(n_nodes)

    INTEGER :: i_part , i_sect

    REAL(wp) :: Ml , Mr

    ! Compute the quadrature abscissas and weights on the interval [-1;1]
    CALL gaulegf(-1.0_wp, 1.0_wp, x, w, n_nodes)

    IF ( verbose_level .GE. 1 ) THEN

       WRITE(*,*) 'original quadrature points in [-1;1]'
       WRITE(*,*) x
       WRITE(*,*) w
       READ(*,*)

    END IF
    
    DO i_part=1,n_part
    
       DO i_sect=1,n_sections

          ml = m(i_sect,i_part)
          mr = m(i_sect+1,i_part)

          ! https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval
          m_quad(:,i_sect,i_part) = 0.5_wp * ( ( mr - ml ) * x + ( mr + ml ) )
          w_quad(:,i_sect,i_part) = 0.5_wp * ( mr - ml ) * w

          IF ( verbose_level .GE. 1 ) THEN

             WRITE(*,*) 'i_part,i_sect',i_part,i_sect
             WRITE(*,*) 'Ml,Mr',ml,mr
             WRITE(*,"(100ES12.4)") m_quad(:,i_sect,i_part)
             WRITE(*,"(100ES12.3)") w_quad(:,i_sect,i_part)
             READ(*,*)
             
          END IF
        
       END DO
    
    END DO

  END SUBROUTINE init_quadrature_points

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE phifromm

    IMPLICIT NONE

    REAL(wp) :: diam1 , diam2
    REAL(wp) :: Vol1 , Vol2
    REAL(wp) :: M1 , M2
    REAL(wp) :: Vol(n_nodes,n_sections)
    REAL(wp) :: f1(n_nodes,n_sections) , f2(n_nodes,n_sections)
    REAL(wp) :: f(n_nodes,n_sections)
    REAL(wp) :: phiL(n_nodes,n_sections) , phiR(n_nodes,n_sections)
    REAL(wp) :: phi(n_nodes,n_sections)
    REAL(wp) :: diam(n_nodes,n_sections)
    REAL(wp) :: Mi(n_nodes,n_sections)
    REAL(wp) :: rho_p(1:n_nodes,1:n_sections)
    REAL(wp) :: shape_p(1:n_nodes,1:n_sections)
    REAL(wp) :: cond1(n_nodes,n_sections), cond2(n_nodes,n_sections)
    REAL(wp) :: cond3(n_nodes,n_sections)
    
    INTEGER :: i_part
    INTEGER :: i_sect
    INTEGER :: iter
    INTEGER :: j
    
    DO i_part=1,n_part

       ! convert from Krumbein scale to meters
       diam1 = 1.e-3_wp * 2.0_wp**( - phi2(i_part) )
       ! compute the volume from diameter and shape factor
       ! Vol1 = shape_factor(i_part) * diam1**3
       vol1 = shape2(i_part) * diam1**3
       ! compute the mass from volume and density
       m1 = vol1 * rho2(i_part)

       ! convert from Krumbein scale to meters
       diam2 = 1.e-3_wp * 2.0_wp**( - phi1(i_part) )
       ! compute the volume from diameter and shape factor
       ! Vol2 = shape_factor(i_part) * diam2**3
       vol2 = shape1(i_part) * diam2**3
       ! compute the mass from volume and density
       m2 = vol2 * rho1(i_part)
       
       ! Initialize the functions for the bisection method
       f1 = m_quad(1:n_nodes,1:n_sections,i_part) - m1
       f2 = m_quad(1:n_nodes,1:n_sections,i_part) - m2

       IF ( verbose_level .GE. 2 ) THEN
       
          WRITE(*,*) 'm_quad'
          DO i_sect=1,n_sections
             
             WRITE(*,*) m_quad(1:n_nodes,i_sect,i_part)
             
          END DO
          
          WRITE(*,*) 'f1'
          DO i_sect=1,n_sections
             
             WRITE(*,*) f1(1:n_nodes,i_sect)
             
          END DO
          
          WRITE(*,*) 'f2'
          DO i_sect=1,n_sections
             
             WRITE(*,*) f2(1:n_nodes,i_sect)
             
          END DO

          WRITE(*,*) 'f1*f2'
          DO i_sect=1,n_sections
             
             WRITE(*,*) f1(1:n_nodes,i_sect)*f2(1:n_nodes,i_sect)
             
          END DO
          
          READ(*,*)

       END IF

       ! Initialize the arrays of the left and right guess for the solution
       phil(1:n_nodes,1:n_sections) = phi2(i_part)
       phir(1:n_nodes,1:n_sections) = phi1(i_part)
       
       DO iter=1,30
          
          phi  = 0.5_wp * (phil+phir)
          diam = 1.e-3_wp * 2.0_wp**( - phi )
          ! Vol = shape_factor(i_part) * diam**3

          shape_p = shape1(i_part) + ( phi - phi1(i_part) ) / ( phi2(i_part) -  &
               phi1(i_part) ) * ( shape2(i_part) - shape1(i_part) )
          shape_p = 0.6_wp
          vol = shape_p * diam**3

          rho_p = rho1(i_part) + ( phi - phi1(i_part) ) / ( phi2(i_part) -      &
               phi1(i_part) ) * ( rho2(i_part) - rho1(i_part) )
    
          mi = vol * rho_p
          f = m_quad(1:n_nodes,1:n_sections,i_part) - mi
          
          ! The phi-rho relationship is linear-piecewise, so we need some
          ! condition defining in which part we are
          cond1 = merge( 1.0_wp , 0.0_wp , f1*f2 .LT. 0.0_wp )
          cond2 = merge( 1.0_wp , 0.0_wp , f*f1 .LT. 0.0_wp )
          cond3 = merge( 1.0_wp , 0.0_wp , f1 .GT. 0.0_wp )
          
          phir = cond1 * ( cond2 * phi + (1.0_wp-cond2) * phir ) +                &
               (1.0_wp-cond1) * ( cond3*phir + (1.0_wp-cond3)*phi )

          f2 = cond1 * ( cond2 * f + (1.0_wp-cond2) * f2 ) +                      &
               (1.0_wp-cond1) * ( cond3 * f2 + (1.0_wp-cond3) * f )
          
          phil = cond1 * ( (1.0_wp-cond2) * phi + cond2 * phil ) +                &
               (1.0_wp-cond1) * ( (1.0_wp-cond3) * phil + cond3 * phi )

          f1 = cond1 * ( (1.0_wp-cond2) * f + cond2 * f1 ) +                      &
               (1.0_wp-cond1) * ( (1.0_wp-cond3) * f1 + cond3 * f )
          
       END DO

       ! the output of the bilinear iterations is the density
       rho_quad(1:n_nodes,1:n_sections,i_part) = rho_p

       shape_quad(1:n_nodes,1:n_sections,i_part) = shape_p
       
       ! we compute the volume from mass and density
       vol_quad(1:n_nodes,1:n_sections,i_part) = m_quad(:,:,i_part) / rho_p

       ! the diameter in m is computed from volume and shapefactor
       diam_quad(1:n_nodes,1:n_sections,i_part) = ( vol_quad(:,:,i_part)        &
            / shape_quad(1:n_nodes,1:n_sections,i_part) )**( 1.0_wp/3.0_wp )

       ! the diameter in m is converted to phi
       phi_quad(1:n_nodes,1:n_sections,i_part) =                                &
            - log(1.e3_wp * diam_quad(1:n_nodes,1:n_sections,i_part)) / log(2.0_wp)

       DO i_sect=1,n_sections
       
          DO j=1,n_nodes
             
             cp_quad(j,i_sect,i_part) = particles_heat_capacity( i_part,        &
                  diam_quad(j,i_sect,i_part))  
             
          END DO

       END DO
       
       IF ( verbose_level .GE. 1 ) THEN

          WRITE(*,*) 'i_part',i_part
          WRITE(*,*) 'phi'
          WRITE(*,"(100ES8.1)") phi_quad(1:n_nodes,1:n_sections,i_part)
          WRITE(*,*) 'diam'
          WRITE(*,"(100ES8.1)") diam_quad(1:n_nodes,1:n_sections,i_part)
          WRITE(*,*) 'vol'
          WRITE(*,"(100ES8.1)") vol_quad(1:n_nodes,1:n_sections,i_part)
          WRITE(*,*) 'rho'
          WRITE(*,"(100ES8.1)") rho_quad(1:n_nodes,1:n_sections,i_part)
          READ(*,*)

          END IF
       
    END DO
    
    RETURN
    
  END SUBROUTINE phifromm

  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE update_aggregation(temp,visc,lw_vf,ice_vf,solid_mf)

    IMPLICIT NONE

    REAL(wp),INTENT(IN) :: temp
    REAL(wp),INTENT(IN) :: visc
    REAL(wp),INTENT(IN) :: lw_vf
    REAL(wp),INTENT(IN) :: ice_vf
    REAL(wp),INTENT(IN) :: solid_mf
  
    INTEGER :: i_part , j_part
    INTEGER :: i_sect , j_sect , k_sect
    
    INTEGER :: i_node , j_node
    
    REAL(wp) :: kernel_ji(n_nodes,n_nodes)

    REAL(wp) :: f1i(n_nodes)
    REAL(wp) :: f2j(n_nodes)
    REAL(wp) :: integrand_ijkm
    REAL(wp) :: f1i_f2j(n_nodes,n_nodes)


    DO i_part=1,n_part 

       DO i_sect=1,n_sections

          DO j_part=1,n_part 

             DO j_sect=1,n_sections

                DO i_node=1,n_nodes

                   DO j_node=1,n_nodes

                      ! If the kernel is a function of particle sizes only, then
                      ! it can be initialized before the integration of the
                      ! plume equation
                      kernel_ji(j_node,i_node) = particles_beta( temp , visc ,  &
                           diam_quad(j_node,j_sect,j_part) ,                    &
                           diam_quad(i_node,i_sect,i_part) ,                    &
                           rho_quad(j_node,j_sect,j_part) ,                     &
                           rho_quad(i_node,i_sect,i_part) ,                     &
                           set_vel_quad(j_node,j_sect,j_part) ,                 &
                           set_vel_quad(i_node,i_sect,i_part) ,                 &
                           lw_vf , ice_vf , solid_mf )

                   END DO

                END DO
                
                kernel_aggr(:,:,j_sect,j_part,i_sect,i_part) = kernel_ji

             END DO

          END DO

       END DO

    END DO

    !--- Birth and death terms due to aggregation
    birth_mom(0:n_mom-1,1:n_sections,1:n_part) = 0.0_wp
    death_mom(0:n_mom-1,1:n_sections,1:n_part) = 0.0_wp

    ! loop over all particles which aggregate i_part
    DO i_part=1,n_part 
       
       ! loop over sections of i_part-particles
       DO i_sect=1,n_sections
          
          ! linear reconstruction at nodes of section i_sect of part(i_part)
          f1i = f_quad(:,i_sect,i_part)
          
          ! loop over all particles which aggregate part(jp)
          DO j_part=1,n_part 

             ! loop over sections of j_part-particles
             DO j_sect=1,n_sections 
         
                f2j = f_quad(:,j_sect,j_part)
                       
                ! interval of particles of family jp
                DO k_sect=1,n_sections             
                   
                   ! check for combination of sections aggregation
                   IF ( q_flag(k_sect,j_sect,i_sect,j_part,i_part) ) THEN

                      DO i_node=1,n_nodes

                         ! term accounting for the linear recontruction
                         ! values at the quadrature nodes
                         f1i_f2j(i_node,1:n_nodes) = f1i(1:n_nodes) * f2j(i_node)

                      END DO

                      integrand_ijkm = sum(                                     &
                           a(:,:,k_sect,j_sect,j_part,i_sect,i_part)            &
                           * wij(:,:,0,j_sect,j_part,i_sect,i_part)             &
                           * kernel_aggr(:,:,j_sect,j_part,i_sect,i_part)       &
                           * f1i_f2j )
                      
                      ! 0-moment of loss rate of i_part-particles in section 
                      ! i_sect becauseof aggregation with j_part-particles in 
                      ! section j_sect to form n_part-particles in section k_sect
                      death_mom(0,i_sect,i_part) = death_mom(0,i_sect,i_part)   &
                           + integrand_ijkm

                      ! 0-moment of birth rate of n_part-particles in section
                      ! k_sect because of aggregation of i_part-particles in
                      ! section i_sect with j_part-particles in section j_sect
                      birth_mom(0,k_sect,n_part) = birth_mom(0,k_sect,n_part)   &
                           + 0.5_wp * integrand_ijkm
                      
                      integrand_ijkm = sum(                                     &
                           a(:,:,k_sect,j_sect,j_part,i_sect,i_part)            &
                           * wij(:,:,1,j_sect,j_part,i_sect,i_part)             &
                           * kernel_aggr(:,:,j_sect,j_part,i_sect,i_part)       &
                           * f1i_f2j )
                      
                      ! 1-moment of loss rate of i_part-particles in section 
                      ! i_sect becauseof aggregation with j_part-particles in 
                      ! section j_sect to form n_part-particles in section k_sect
                      death_mom(1,i_sect,i_part) = death_mom(1,i_sect,i_part)   &
                           + integrand_ijkm

                      ! 1-moment of birth rate of n_part-particles in section
                      ! k_sect because of aggregation of i_part-particles in
                      ! section i_sect with j_part-particles in section j_sect
                      birth_mom(1,k_sect,n_part) = birth_mom(1,k_sect,n_part)   &
                           + integrand_ijkm

                   END IF
                        
                END DO
                    
             END DO
                
          END DO
            
       END DO

    END DO

    RETURN
    
  END SUBROUTINE update_aggregation


  !******************************************************************************
  !
  !******************************************************************************
  
  SUBROUTINE init_aggregation

    IMPLICIT NONE

    REAL(wp) :: mi1(n_nodes) , mj2(n_nodes)
    REAL(wp) :: wi1(n_nodes) , wj2(n_nodes)

    REAL(wp) :: mij12(n_nodes,n_nodes)
    REAL(wp) :: wij12(n_nodes,n_nodes)
    INTEGER :: Aijk(n_nodes,n_nodes)
    REAL(wp) :: Wijm(n_nodes,n_nodes)

    REAL(wp) :: kernel_ji(n_nodes,n_nodes)
    
    INTEGER :: i_part , j_part
    INTEGER :: i_sect , j_sect , k_sect
    INTEGER :: i_mom
    INTEGER :: i_node , j_node
    
    DO i_part=1,n_part 

       DO i_sect=1,n_sections

          mi1(1:n_nodes) = m_quad(:,i_sect,i_part)
          wi1(1:n_nodes) = w_quad(:,i_sect,i_part)

          DO j_part=1,n_part 

             DO j_sect=1,n_sections

                mj2(1:n_nodes) = m_quad(:,j_sect,j_part)
                wj2(1:n_nodes) = w_quad(:,j_sect,j_part)

                DO i_node=1,n_nodes

                   mij12(i_node,1:n_nodes) = mi1(1:n_nodes)+mj2(i_node)
                   wij12(i_node,1:n_nodes) = wi1(1:n_nodes)*wj2(i_node)

                END DO
               
                DO k_sect=1,n_sections

                   ! The array Aijk depends only on the sections i, j and k.
                   ! The value is 1 when the sum of the mass of the particles
                   ! i_part,i_sect,i_node and j_part,j_sect,j_node falls in
                   ! the k_section of the class n_part
                   aijk = merge(1 , 0 , ( mij12 .GE. m(k_sect,n_part) ) .AND.   &
                        ( mij12 .LE. m(k_sect+1,n_part) ) )

                   ! this logical variable is true only when there are particles 
                   ! i_part/i_sect and particles j_part/j_sect aggregating
                   ! to particles n_part/k_sect
                   q_flag(k_sect,j_sect,i_sect,j_part,i_part) =                 &
                        ( sum(aijk) .GT. 0 )
 
                   a(:,:,k_sect,j_sect,j_part,i_sect,i_part) = aijk
                   
                END DO

                DO i_mom=0,n_mom-1

                   DO i_node=1,n_nodes

                      wijm(i_node,1:n_nodes) = wij12(i_node,1:n_nodes)          &
                           * mi1(1:n_nodes)**i_mom

                   END DO

                   wij(:,:,i_mom,j_sect,j_part,i_sect,i_part) = wijm
                      
                END DO

             END DO

          END DO

       END DO

    END DO

    WRITE(*,*) 'Aggregation initialization completed'
    
    RETURN

  END SUBROUTINE init_aggregation
  
END MODULE particles_module