generate_particles_mod.f

Go to the documentation of this file.

      MODULE generate_particles

        DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: particle_count

      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                C
!                                                                      C
!  Purpose: Generate particle configuration based on maximum particle  C
!           radius and filling from top to bottom within specified     C
!           bounds                                                     C
!                                                                      C
!                                                                      C
!  Authors: Rahul Garg                                Date: 19-Mar-14  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE generate_particle_config

      use mfix_pic, only: mppic
      use discretelement, only: pip, particles
! Flag indicating that the IC region is defined.
      use ic, only: ic_defined
! Parameter for detecting unspecified values, zero, and one
      use param1, only: one
! Maximum number of initial conditions
      use param, only: dimension_ic
! IC Region gas volume fraction.
      use ic, only: ic_ep_g

      use mpi_utility

! Use the error manager for posting error messages.
!---------------------------------------------------------------------//
      use error_manager

      IMPLICIT NONE

      INTEGER :: ICV

! Initialize the error manager.
      CALL init_err_msg("Generate_Particle_Config")

      DO icv = 1, dimension_ic

         IF(.NOT.ic_defined(icv)) cycle
         IF(ic_ep_g(icv) == one) cycle

         IF(mppic) THEN
            CALL generate_particle_config_mppic(icv)
         ELSE
            CALL generate_particle_config_dem(icv)
         ENDIF

      ENDDO

      CALL global_sum(pip,particles)

      WRITE(err_msg, 1004) particles
 1004 FORMAT(/,'Total number of particles in the system: ',i15)

      CALL flush_err_msg(header=.false., footer=.false.)

      CALL finl_err_msg

      RETURN
      END SUBROUTINE generate_particle_config

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                !
!  Authors: Rahul Garg                               Date: 21-Mar-2014 !
!                                                                      !
!  Purpose: Generate particle configuration for DEM solids for each IC !
!           region. Now using the particle linked lists for initial    !
!           build                                                      !
!           This routine will ultimately supersede the older rouine    !
!           that has not been deleted yet                              !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE generate_particle_config_dem(ICV)

! Global Variables:
!---------------------------------------------------------------------//
! particle radius and density
      use discretelement, only: des_radius, ro_sol
! particle position new and old
      use discretelement, only: des_pos_new, des_pos_old
! particle velocity new and old
      use discretelement, only: des_vel_new, des_vel_old
! Simulation dimension (2D/3D)
      use discretelement, only: dimn
! Number of particles in the system (current)
      use discretelement, only: pip
! Number of DEM solids phases.
      use discretelement, only: des_mmax
! Flag to use _OLD variables
      use discretelement, only: do_old
! Angular velocity
      use discretelement, only: omega_old, omega_new, pijk
! solid phase diameters and densities.
      use physprop, only: d_p0, ro_s0, mmax
! IC Region solids volume fraction.
      use ic, only: ic_ep_s

! Constant: 3.14159...
      use constant, only: pi
! min and max physical co-ordinates of IC regions in each direction
      use ic, only: ic_x_w, ic_x_e, ic_y_s, ic_y_n, ic_z_b, ic_z_t
! initally specified velocity field and granular temperature
      use ic, only: ic_u_s, ic_v_s, ic_w_s, ic_theta_m
! Flag to extend the lattice distribution in a given IC to available area
      use ic, only: ic_des_fit_to_region
! Parameter for detecting unspecified values, zero, and one
      use param1, only: undefined, undefined_i, zero, one, half
! Parameter for small and large numbers
      use param1, only: small_number, large_number

! to access random number generator subroutines
      use randomno
      use mpi_utility
      use functions, only: set_normal

      use desgrid, only: dg_xstart, dg_ystart, dg_zstart
      use desgrid, only: dg_xend, dg_yend, dg_zend

! direction wise spans of the domain and grid spacing in each direction
      use geometry, only: xlength, ylength, zlength

      use cutcell, only : cartesian_grid
      use stl_functions_des, only: check_if_particle_overlaps_stl
      use run, only: solids_model
      use des_allocate, only: particle_grow

      use desgrid, only: iofpos, jofpos, kofpos
      use desgrid, only: dg_is_on_mype_owns
      use toleranc, only: compare

      use discretelement, only: max_pip, max_radius, xe, yn, zt
      use error_manager
      use functions
      use param, only: dim_m
      use param, only: dimension_i, dimension_j, dimension_k

      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      INTEGER, INTENT(IN) :: ICV

! Local variables
!---------------------------------------------------------------------//
! Starting positions in the axial directions
      DOUBLE PRECISION :: xINIT, yINIT, zINIT
! Fractor used to scale particle diameter
      DOUBLE PRECISION :: lFAC
! Particle position and velocity
      DOUBLE PRECISION :: POS(3), VEL(3)
! Number of particles in the lattice
      INTEGER :: SEED_X, SEED_Y, SEED_Z
! Loop indices phase/fluid cell
      INTEGER :: M, MM, I, J, K, IJK, LB, UB
! Loop indicies for seeding
      INTEGER :: II, JJ, KK
! Start and end bound for IC region.
      DOUBLE PRECISION :: IC_START(3), IC_END(3)
! Volume and lengths of the IC Region
      DOUBLE PRECISION :: DOM_VOL, DOML(3)
! Flag to skip the particle
      LOGICAL :: SKIP
! Diameter adjusted for space padding
      DOUBLE PRECISION :: ADJ_DIA
! Number of particles calculated from volume fracton
      INTEGER :: rPARTS(dim_m), tPARTS
! Spacing between particles.
      DOUBLE PRECISION :: lDEL, lDX, lDY, lDZ
! Flag that the setup failed to fit the particles to the IC region
      LOGICAL :: FIT_FAILED
! Number of seeded particles
      INTEGER :: pCOUNT(dim_m), tCOUNT

      DOUBLE PRECISION :: SOLIDS_DATA(0:dim_m)

      LOGICAL :: VEL_FLUCT
      DOUBLE PRECISION :: VEL_SIG
      DOUBLE PRECISION, ALLOCATABLE :: randVEL(:,:)

      logical :: report = .true.
      logical :: found

!......................................................................!

      CALL init_err_msg("GENERATE_PARTICLE_CONFIG_DEM")

      WRITE(err_msg,"(2/,'Generating initial particle configuration:')")
      CALL flush_err_msg(header=.false., footer=.false.)

      solids_data = zero
      CALL get_ic_volume(icv, solids_data(0))

! setting particle seed spacing grid to be slightly greater than
! the maximum particle diameter. seed at ~particle radii
      lfac = 1.05d0

! Setup local arrays with IC region bounds.
      ic_start(1)=ic_x_w(icv);   ic_end(1)=ic_x_e(icv)
      ic_start(2)=ic_y_s(icv);   ic_end(2)=ic_y_n(icv)
      ic_start(3)=ic_z_b(icv);   ic_end(3)=ic_z_t(icv)

      doml = ic_end-ic_start
      IF(no_k) doml(3)=dz(1)

! Volume of the IC region
      dom_vol = doml(1)*doml(2)*doml(3)

      rparts=0
      DO m=mmax+1,mmax+des_mmax
         IF(solids_model(m) == 'DEM') THEN
! Number of particles for phase M
            rparts(m) = &
               floor((6.0d0*ic_ep_s(icv,m)*dom_vol)/(pi*(d_p0(m)**3)))
         ENDIF
      ENDDO

! Total number of particles in this IC region.
      tparts = sum(rparts)
      IF(tparts == 0) RETURN

      adj_dia = 2.0d0*max_radius*lfac

! Attempt to seed particle throughout the IC region
      fit_failed=.false.
      IF(ic_des_fit_to_region(icv)) THEN
         IF(no_k) THEN
            ldel = (doml(1)-adj_dia)*(doml(2)-adj_dia)
            ldel = (ldel/dble(tparts))**(1.0/2.0)
            seed_x = max(1,ceiling((doml(1)-adj_dia)/ldel))
            seed_y = max(1,ceiling((doml(2)-adj_dia)/ldel))
            seed_z = 1
         ELSE
            ldel = (doml(1)-adj_dia)*(doml(2)-adj_dia)*(doml(3)-adj_dia)
            ldel = (ldel/dble(tparts))**(1.0/3.0)
            seed_x = max(1,ceiling((doml(1)-adj_dia)/ldel))
            seed_y = max(1,ceiling((doml(2)-adj_dia)/ldel))
            seed_z = max(1,ceiling((doml(3)-adj_dia)/ldel))
         ENDIF
         fit_failed=(dble(seed_x*seed_y*seed_z) < tparts)
      ENDIF

! Generic filling
      IF(.NOT.ic_des_fit_to_region(icv) .OR. fit_failed) THEN
         seed_x = max(1,floor((ic_end(1)-ic_start(1)-adj_dia)/adj_dia))
         seed_y = max(1,floor((ic_end(2)-ic_start(2)-adj_dia)/adj_dia))
         seed_z = max(1,floor((ic_end(3)-ic_start(3)-adj_dia)/adj_dia))
      ENDIF

      ldx = doml(1)/dble(seed_x)
      ldy = doml(2)/dble(seed_y)
      IF(do_k) THEN
         ldz = doml(3)/dble(seed_z)
      ELSE
         ldz = 0.0d0
      ENDIF

      xinit = ic_start(1)+half*ldx
      yinit = ic_start(2)+half*ldy
      zinit = ic_start(3)+half*ldz

      m=1
      pcount = 0
      tcount = 0

      vel_fluct = set_vel_fluct(icv,m)

      jj_lp: DO jj=1, seed_y
         pos(2) = yinit + (jj-1)*ldy
         IF(compare(pos(2),dg_ystart) .OR. compare(pos(2),dg_yend))    &
            pos(2) = pos(2) + small_number

      kk_lp: DO kk=1, seed_z
         pos(3) = zinit + (kk-1)*ldz
         IF(do_k) THEN
            IF(compare(pos(3),dg_zstart) .OR. compare(pos(3),dg_zend)) &
               pos(3) = pos(3) + small_number
         ENDIF

      ii_lp: DO ii=1, seed_x
         pos(1) = xinit + (ii-1)*ldx
         IF(compare(pos(1),dg_xstart) .OR. compare(pos(1),dg_xend))    &
            pos(1) = pos(1) + small_number

! Exit if all particles were seeded.
         IF(tcount > int(tparts)) THEN
            EXIT jj_lp
! Find the next phase that needs to be seeded
         ELSEIF(pcount(m) > int(rparts(m))) THEN
            mm_lp: DO mm=m+1,mmax+des_mmax
               IF(rparts(mm) > 0.0) THEN
                  m=mm
                  EXIT mm_lp
               ENDIF
            ENDDO mm_lp
            IF(m > mmax+des_mmax) EXIT jj_lp
            vel_fluct = set_vel_fluct(icv,m)
         ENDIF

         pcount(m) = pcount(m) + 1
         tcount = tcount + 1

! Keep only particles that belong to this process.
         IF(.NOT.dg_is_on_mype_owns(iofpos(pos(1)), &
            jofpos(pos(2)),kofpos(pos(3)))) cycle

! Bin the parcel to the fuild grid.
         k=1
         IF(do_k) CALL pic_search(k, pos(3), zt, dimension_k, kmin2, kmax2)
         CALL pic_search(j, pos(2), yn, dimension_j, jmin2, jmax2)
         CALL pic_search(i, pos(1), xe, dimension_i, imin2, imax2)

! Skip cells that return invalid IJKs.
         IF(dead_cell_at(i,j,k)) cycle

! Skip cells that are not part of the local fuild domain.
         ijk = funijk(i,j,k)
         IF(.NOT.fluid_at(ijk)) cycle

         IF(cartesian_grid) THEN
            CALL check_if_particle_overlaps_stl(pos, i, j, k, skip)
            IF(skip) cycle
         ENDIF

         pip = pip + 1
         CALL particle_grow(pip)
         max_pip = max(pip,max_pip)

         CALL set_normal(pip)

         IF(vel_fluct) THEN
            vel(1) = randvel(pcount(m),1)
            vel(2) = randvel(pcount(m),2)
            vel(3) = randvel(pcount(m),3)
         ELSE
            vel(1) = ic_u_s(icv,m)
            vel(2) = ic_v_s(icv,m)
            vel(3) = ic_w_s(icv,m)
         ENDIF
         IF(no_k) vel(3) = 0.0d0


         des_pos_new(pip,:) = pos(:)
         des_vel_new(pip,:) = vel(:)
         omega_new(pip,:) = 0.0d0

         des_radius(pip) = d_p0(m)*half
         ro_sol(pip) =  ro_s0(m)

         pijk(pip,1) = i
         pijk(pip,2) = j
         pijk(pip,3) = k
         pijk(pip,4) = ijk
         pijk(pip,5) = m

         IF(do_old) THEN
            des_vel_old(pip,:) = des_vel_new(pip,:)
            des_pos_old(pip,:) = des_pos_new(pip,:)
            omega_old(pip,:) = zero
         ENDIF

         solids_data(m) = solids_data(m) + 1.0

      ENDDO ii_lp
      ENDDO kk_lp
      ENDDO jj_lp

! Collect the data
      CALL global_all_sum(solids_data)

! Verify that the IC region volume is not zero.
      IF(solids_data(0) <= 0.0d0) THEN
         WRITE(err_msg,1000) icv, solids_data(0)
         CALL flush_err_msg(abort=.true.)
      ENDIF

1000 FORMAT('Error 1000: Invalid IC region volume: IC=',i3,' VOL=',&
         es15.4,/'Please correct the mfix.dat file.')

      WRITE(err_msg,2000) icv
      CALL flush_err_msg(header=.false., footer=.false.)

      DO m=mmax+1, mmax+des_mmax
         IF(solids_data(m) < small_number) cycle
         WRITE(err_msg,2010) m, int(solids_data(m)), ic_ep_s(icv,m),   &
            (dble(solids_data(m))*(pi/6.0d0)*d_p0(m)**3)/solids_data(0)
         CALL flush_err_msg(header=.false., footer=.false.)
      ENDDO

      IF(allocated(randvel)) deallocate(randvel)

      CALL finl_err_msg

      RETURN

 2000 FORMAT(/2x,'|',43('-'),'|',/2x,'| IC Region: ',i3,28x,'|',/2x,   &
         '|',43('-'),'|',/2x,'| Phase | Number of |    EPs    |    EP',&
         's    |',/2x,'|   ID  | Particles | Specified |   Actual  |', &
         /2x,'|-------|',3(11('-'),'|'))

 2010 FORMAT(2x,'|  ',i3,'  |',1x,i9,1x,'|',2(1x,es9.2,1x,'|'),/2x,    &
         '|-------|',3(11('-'),'|'))


      CONTAINS

!......................................................................!
! Function: SET_VEL_FLUCT                                              !
! Purpose: Set the flag for random velocity fluctuations. If needed    !
! the random velocities are calculated.                                !
!......................................................................!
      LOGICAL FUNCTION set_vel_fluct(lICV, lM)
      INTEGER, INTENT(IN) :: lICV, lM
      DOUBLE PRECISION :: VEL_SIG
      set_vel_fluct=(ic_theta_m(licv,lm) /= 0.0d0)
      IF(set_vel_fluct) THEN
         if(allocated(randvel)) deallocate(randvel)
         allocate(randvel(100+int(rparts(lm)),3))
         vel_sig = sqrt(ic_theta_m(licv,lm))
         CALL nor_rno(randvel(:,1), ic_u_s(licv,lm),vel_sig)
         CALL nor_rno(randvel(:,2), ic_v_s(licv,lm),vel_sig)
         IF(do_k) CALL nor_rno(randvel(:,3),ic_w_s(licv,lm),vel_sig)
      ENDIF

      END FUNCTION set_vel_fluct

      END SUBROUTINE generate_particle_config_dem




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
!  Author: Rahul Garg                                 Date: 3-May-2011 !
!                                                                      !
!  Purpose: Generates particle position distribution for MPPIC.        !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE generate_particle_config_mppic(ICV)

! Global variables
!---------------------------------------------------------------------//
! Number of DES solids phases.
      use discretelement, only: des_mmax
! Flag indicating that the IC region is defined.
      use ic, only: ic_defined
! IC Region bulk density (RO_s * EP_s)
      use ic, only: ic_rop_s
! IC Region gas volume fraction.
      use ic, only: ic_ep_g
! MPPIC specific IC region specification.
      use ic, only: ic_pic_const_npc, ic_pic_const_statwt

      use param1, only: undefined, undefined_i
      use param1, only: zero, one, half
! Maximum number of IC regions and solids phases
      use param, only: dimension_ic
      use param, only: dim_m
      use physprop, only: mmax

! The accumulated number of particles in each IJK.
      use mpi_utility, only: global_all_sum

      use error_manager

      use run, only: solids_model
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      INTEGER, INTENT(IN) :: ICV

! Local variables
!---------------------------------------------------------------------//
! Generic loop counters
      INTEGER :: M
! Actual volume of IC region
      DOUBLE PRECISION :: IC_VOL
! Solids data in IC Region by phase:
      DOUBLE PRECISION :: SOLIDS_DATA(0:4*dim_m)
!......................................................................!

      CALL init_err_msg("GENERATE_PARTICLE_CONFIG_MPPIC")

      WRITE(err_msg,"(2/,'Generating initial parcel configuration:')")
      CALL flush_err_msg(header=.false., footer=.false.)


      solids_data = zero
      CALL get_ic_volume(icv, solids_data(0))

      WRITE(err_msg,2000) icv
      CALL flush_err_msg(header=.false., footer=.false.)

! Set up the individual solids phases.
      DO m=mmax+1, des_mmax+mmax
         IF(solids_model(m) == 'PIC') THEN
            IF(ic_rop_s(icv,m) == zero) cycle
! Seed parcels with constant stastical weight.
            CALL gpc_mppic_const_npc(icv, m, solids_data(0), &
            solids_data((4*m-3):(4*m)))
         ENDIF
      ENDDO

! Collect the data
      CALL global_all_sum(solids_data)

! Verify that the IC region volume is not zero.
      IF(solids_data(0) <= 0.0d0) THEN
         WRITE(err_msg,1000) icv, solids_data(0)
         CALL flush_err_msg(abort=.true.)
      ENDIF

 1000 FORMAT('Error 1000: Invalid IC region volume: IC=',i3,' VOL=',&
         es15.4,/'Please correct the mfix.dat file.')

! Report solids information for the IC region.
      DO m=mmax+1, des_mmax+mmax
         IF(solids_model(m) == 'PIC') THEN
            WRITE(err_msg,2010) m, int(solids_data(4*m-3)),&
               int(solids_data(4*m-3)*solids_data(4*m-2)), &
               solids_data(4*m-2), solids_data(4*m-1),     &
               solids_data(4*m)/solids_data(0)
            CALL flush_err_msg(header=.false., footer=.false.)
         ENDIF
      ENDDO

      CALL finl_err_msg

      RETURN

 2000 FORMAT(/2x,'|',67('-'),'|',/2x,'| IC Region: ',i3,52x,'|',/2x,   &
         '|',67('-'),'|',/2x,'| Phase | Num. Comp | Num. Real ',       &
         '| Stastical |    EPs    |    EPs    |',/2x,'|   ID  |  ',    &
         'Parcels  | Particles |   Weight  | Specified |   Actual  |', &
         /2x,'|-------|',5(11('-'),'|'))

 2010 FORMAT(2x,'|  ',i3,'  |',2(1x,i9,1x,'|'),3(1x,es9.2,1x,'|'),/2x,&
         '|-------|',5(11('-'),'|'))

      END SUBROUTINE generate_particle_config_mppic



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
!  Author: Rahul Garg                                 Date: 3-May-2011 !
!                                                                      !
!  Purpose: generates particle position distribution for MPPIC         !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE gpc_mppic_const_npc(ICV, M, IC_VOL, sDATA)

! Global variables
!---------------------------------------------------------------------//
! Constant: 3.14159...
      use constant, only: pi
! Cut_cell identifier array
      use cutcell, only: cut_cell_at
      use cutcell, only : cartesian_grid
      use discretelement, only: xe, yn, zt

      use des_allocate, only: particle_grow
      use discretelement

! Flag indicating that the IC region is defined.
      use ic, only: ic_defined
! IC Region bulk density (RO_s * EP_s)
      use ic, only: ic_ep_s
      use ic, only: ic_i_w, ic_i_e, ic_j_s, ic_j_n, ic_k_b, ic_k_t

! initally specified velocity field and granular temperature
      use ic, only: ic_u_s, ic_v_s, ic_w_s, ic_theta_m
      use ic, only: ic_pic_const_npc
      use ic, only: ic_pic_const_statwt

      use mfix_pic, only: des_stat_wt
      use mpi_utility

! Maximum number of IC regions and solids phases
      use param, only: dimension_ic
      use param, only: dim_m
      use param1, only: undefined, undefined_i, zero, one, half

! solid phase diameters and densities.
      use physprop, only: d_p0, ro_s0
      use run, only: solids_model

      use randomno
      use error_manager
      use functions
      use run, only: solids_model
      use des_allocate, only: particle_grow

      IMPLICIT NONE

! Dummy Arguments
!----------------------------------------------------------------------//
! Index of IC region and solids phase
      INTEGER, INTENT(IN) :: ICV, M
! Specific volume of IC region (accounts for blocked cells)
      DOUBLE PRECISION, INTENT(IN) :: IC_VOL
! Data about solids in the IC region.
      DOUBLE PRECISION, INTENT(OUT) :: sDATA(4)

! Local variables
!----------------------------------------------------------------------//
      DOUBLE PRECISION :: EP_SM

! Number of real and comp. particles in a cell.
      DOUBLE PRECISION ::  rPARTS
      INTEGER :: maxPARTS
      DOUBLE PRECISION :: DOML(3), IC_START(3)
! Parcel position with random
      DOUBLE PRECISION :: POS(3)
! Average velocity and standard deivation
      DOUBLE PRECISION :: IC_VEL(3), VEL_SIG
! Flag to not keep parcel.
      LOGICAL :: SKIP
! Arrasy for assigning random position and velocities
      DOUBLE PRECISION, ALLOCATABLE :: randVEL(:,:)
      DOUBLE PRECISION :: RAND(3)
! Statistical weights
      DOUBLE PRECISION :: STAT_WT, SUM_STAT_WT
! Volume of a parcel and total solids volume
      DOUBLE PRECISION :: sVOL, sVOL_TOT
! Counter for seeded parcles.
      INTEGER :: SEEDED
! Generic loop indices and loop counters
      INTEGER :: I, J, K, IJK, LC, LC_MAX
!......................................................................!

      CALL init_err_msg("GPC_MPPIC_CONST_NPC")

      maxparts=25
      allocate(randvel(maxparts,3))

      ic_vel(1) = ic_u_s(icv,m)
      ic_vel(2) = ic_v_s(icv,m)
      ic_vel(3) = merge(ic_w_s(icv,m),0.0d0,do_k)

      vel_sig = sqrt(ic_theta_m(icv,m))

! Volume occupied by one particle
      svol = (pi/6.0d0)*(d_p0(m)**3.d0)

      seeded = 0
      svol_tot = 0.0d0
      sum_stat_wt = 0.0d0

      DO k = ic_k_b(icv), ic_k_t(icv)
      DO j = ic_j_s(icv), ic_j_n(icv)
      DO i = ic_i_w(icv), ic_i_e(icv)

         IF(.not.is_on_mype_wobnd(i,j,k)) cycle
         IF(dead_cell_at(i,j,k)) cycle

         ijk = funijk(i,j,k)
         IF(.not.fluid_at(ijk)) cycle

         rparts = ic_ep_s(icv,m)*vol(ijk)/svol

! Seed parcels with a constant stastical weight
         IF(ic_pic_const_statwt(icv,m) /= zero) THEN
            stat_wt = ic_pic_const_statwt(icv,m)
            lc_max = int(rparts/stat_wt)

! Seed parcels with a constant number per cell
         ELSEIF(ic_pic_const_npc(icv,m) /= 0) THEN
            lc_max = ic_pic_const_npc(icv,m)
            stat_wt = rparts/dble(lc_max)
            IF(cut_cell_at(ijk)) lc_max = max(1,int(vol(ijk)*dble(&
               ic_pic_const_npc(icv,m))/(dx(i)*dy(j)*dz(k))))
         ENDIF

! Increase particle buffer
         IF(lc_max > maxparts) THEN
            maxparts = 2*lc_max
            if(allocated(randvel)) deallocate(randvel)
            allocate(randvel(maxparts,3))
         ENDIF

         DO lc=1, merge(2,3,no_k)
            IF(vel_sig > zero) THEN
               CALL nor_rno(randvel(1:lc_max,lc), ic_vel(lc), vel_sig)
            ELSE
               randvel(1:lc_max,lc) = ic_vel(lc)
            ENDIF
         ENDDO
         IF(no_k) randvel(1:lc_max,3) = 0.0d0

         ic_start(1) = xe(i-1)
         ic_start(2) = yn(j-1)
         ic_start(3) = zero;  IF(do_k) ic_start(3) = zt(k-1)

         doml(1) = dx(i)
         doml(2) = dy(j)
         doml(3) = zero;  IF(do_k) doml(3) = dz(k)

         DO lc=1,lc_max

            CALL random_number(rand)
            pos(:) = ic_start + doml*rand

            IF(cartesian_grid) THEN
               CALL check_if_parcel_overlaps_stl(pos, skip)
               DO WHILE(skip)
                  CALL random_number(rand)
                  pos(:) = ic_start + doml*rand
                  CALL check_if_parcel_overlaps_stl(pos, skip)
               ENDDO
            ENDIF

            pip = pip + 1
            CALL particle_grow(pip)
            max_pip = max(pip,max_pip)

            des_pos_new(pip,:) = pos(:)
            des_vel_new(pip,:) = randvel(lc,:)

            des_radius(pip) = d_p0(m)*half
            ro_sol(pip) =  ro_s0(m)

            pijk(pip,1) = i
            pijk(pip,2) = j
            pijk(pip,3) = k
            pijk(pip,4) = ijk
            pijk(pip,5) = m

            sum_stat_wt = sum_stat_wt + stat_wt

            des_stat_wt(pip) = stat_wt
            svol_tot = svol_tot + svol*stat_wt

            CALL set_normal(pip)

            seeded = seeded + 1

         ENDDO

      ENDDO
      ENDDO
      ENDDO

      IF(allocated(randvel)) deallocate(randvel)

      sdata(1) = dble(seeded)
      sdata(2) = sum_stat_wt/dble(seeded)
      sdata(3) = ic_ep_s(icv,m)
      sdata(4) = svol_tot

      CALL finl_err_msg

      RETURN
      END SUBROUTINE gpc_mppic_const_npc

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GET_IC_VOLUME                                           !
!  Author: J.Musser                                 Date: 26-Aug-2015  !
!                                                                      !
!  Purpose: Calculate the actual volume of the IC region.              !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE get_ic_volume(ICV, IC_VOL)

! IC region index bounds
      use ic, only: ic_i_w, ic_i_e
      use ic, only: ic_j_s, ic_j_n
      use ic, only: ic_k_b, ic_k_t

! Volume of computational cells.
      use geometry, only: vol

      use functions
      use compar, only: dead_cell_at

      IMPLICIT NONE

! Dummy Arguments
!---------------------------------------------------------------------//
! Index of IC region
      INTEGER, INTENT(IN) :: ICV
! Total calculated volume of IC region
      DOUBLE PRECISION, INTENT(OUT) :: IC_VOL

! Local variables
!---------------------------------------------------------------------//
      INTEGER :: I, J, K, IJK
!......................................................................!


      ic_vol = 0.0d0
      DO k = ic_k_b(icv), ic_k_t(icv)
      DO j = ic_j_s(icv), ic_j_n(icv)
      DO i = ic_i_w(icv), ic_i_e(icv)

         IF(.NOT.is_on_mype_wobnd(i,j,k)) cycle
         IF(dead_cell_at(i,j,k)) cycle

         ijk = funijk(i,j,k)
         IF(fluid_at(ijk)) ic_vol = ic_vol + vol(ijk)

      ENDDO
      ENDDO
      ENDDO

      RETURN
      END SUBROUTINE get_ic_volume

      END MODULE generate_particles