PLAG_ModInjection.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
! ******************************************************************************
!
! Purpose: Suite of routines for stochastic injection algorithm.
!
! Description: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: PLAG_ModInjection.F90,v 1.11 2008/12/06 08:44:34 mtcampbe Exp $
!
! Copyright: (c) 2004-2005 by the University of Illinois
!
! ******************************************************************************

MODULE plag_modinjection
  
  USE modglobal, ONLY: t_global
  USE modparameters
  USE moddatatypes
  USE modglobal,     ONLY: t_global
  USE modgrid,       ONLY: t_grid
  USE modpartlag,    ONLY: t_plag, t_plag_input, t_tile_plag
  USE modbndpatch,   ONLY: t_patch
  USE moddatastruct, ONLY: t_region
  USE moderror
  USE modrandom,     ONLY: rand1uniform

  USE plag_modparameters
  USE inrt_modparameters

  USE plag_modinterfaces, ONLY : plag_injcmakeparticle, &
                                 plag_injcsetinjection

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: plag_ejectparticle, &
            plag_invokeejecmodel1, &
            plag_invokeconsrandejec
        
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
   
  CHARACTER(CHRLEN) :: RCSIdentString = & 
    '$RCSfile: PLAG_ModInjection.F90,v $ $Revision: 1.11 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS



!******************************************************************************
!
! Purpose: Fills PLAG datastructure for ejected particle.
!
! Description: none.
!
! Input:
!   pGlobal      Pointer to global data
!   pPlag        Pointer to particle data
!   pTilePlag    Pointer to tile particle data
!   iTile        Tile index
!   icg          Cell index
!   burnStat     Particle burning status
!   nVert        Number of vertices on tile element surface
!   normDirFlag  Normal direction flag
!   xyzVertex    Coordinates of element vertices
!   fn           Face normal coordinates
!   fc           Face centroids coordinates
!   cellHeight   Cell height
!   plagVolRatio Volume ratio particles get from tiles
!
! Output:
!   pPlag      Plag values for cv, aiv, arv of current region.
!   pTilePlag  Tile values for cv, dv of current region.
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE plag_ejectparticle( pRegion,pPlag,pTilePlag,      &
                               itile,icg,burnstat,           &
                               nvert,normdirflag,xyzvertex,  &
                               fn,fc,cellheight,plagvolratio )

#ifdef RFLU
  USE rflu_modincelltest, ONLY: rflu_ict_testincell
#endif

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

  INTEGER, INTENT(IN) :: burnstat,icg,itile,normdirflag,nvert

  REAL(RFREAL),                      INTENT(IN) :: cellheight,plagvolratio
  REAL(RFREAL), DIMENSION(ZCOORD),   INTENT(IN) :: fc,fn
  REAL(RFREAL), DIMENSION(ZCOORD,4), INTENT(IN) :: xyzvertex     
   
  TYPE(t_plag),      POINTER :: pplag
  TYPE(t_tile_plag), POINTER :: ptileplag
  TYPE(t_region), POINTER  :: pregion

! ==============================================================================
! Locals
! ==============================================================================

  INTEGER, PARAMETER :: iter_cell_locate_max = 20
  INTEGER :: icont,icvplagmass,icvtilemass,ipcl,ireg,itercelllocate, &
             ncont,nextidnumber,npcls,npclsmax
  INTEGER, POINTER, DIMENSION(:)   :: pcvplagmass, pcvtilemass
  INTEGER, POINTER, DIMENSION(:,:) :: paiv

  LOGICAL :: celllocate
   
  REAL(RFREAL) :: plagmomenrm,sxn,syn,szn,xloc,yloc,zloc
  REAL(RFREAL), DIMENSION(3) :: poolxyz
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: parv,pcvplag,pcvtile,pdvtile
  
  TYPE(t_global),    POINTER :: global 
! ******************************************************************************
! Start
! ******************************************************************************

  rcsidentstring = '$RCSfile: PLAG_ModInjection.F90,v $ $Revision: 1.11 $'

  global => pregion%global

  CALL registerfunction( global, 'PLAG_EjectParticle',&
  'PLAG_ModInjection.F90' )

! ******************************************************************************
! Set pointers and variables
! ******************************************************************************

  pcvplag     => pplag%cv
  pcvplagmass => pplag%cvPlagMass
  paiv        => pplag%aiv
  parv        => pplag%arv

  pcvtile     => ptileplag%cv
  pcvtilemass => ptileplag%cvTileMass
  pdvtile     => ptileplag%dv

  ireg     = pregion%iRegionGlobal
  ncont    = pregion%plagInput%nCont

! ******************************************************************************
! Load normals to point inwards
!   Note: RFLU defines normals pointing outwards.
! ******************************************************************************

  sxn = fn(xcoord) *normdirflag
  syn = fn(ycoord) *normdirflag
  szn = fn(zcoord) *normdirflag

! ******************************************************************************
! Increment PLAG datastructure for new injected particle
! ******************************************************************************

  pplag%nPcls = pplag%nPcls + 1
  ipcl        = pplag%nPcls

  ptileplag%nPclsInjc(itile) = ptileplag%nPclsInjc(itile) + 1

! ******************************************************************************
! Trap error if nPcls exceeds maximum datastructure dimension, nPclsMax
! ******************************************************************************

  npcls    = pplag%nPcls
  npclsmax = pregion%plag%nPclsMax 

  IF ( npcls >= npclsmax ) THEN
    WRITE(stdout,*) ' PLAG_EjectParticle: Datastructure Dimension Exceeded ',&
                      npcls,npclsmax
    CALL errorstop( global,err_plag_memoverflow,__line__ )
  END IF ! nPcls

! ******************************************************************************
! Keep track of global count in particle ejected in region
! ******************************************************************************

  pplag%nextIdNumber = pplag%nextIdNumber + 1
  nextidnumber       = pplag%nextIdNumber

! ******************************************************************************
! Compute PLAG cv datastructure for new injected particle
! ******************************************************************************

  DO icont = 1, ncont
    icvplagmass = pcvplagmass(icont)
    icvtilemass = pcvtilemass(icont)
    pcvplag(icvplagmass,ipcl) = pcvtile(icvtilemass,itile)    *plagvolratio
  END DO ! iCont

  pcvplag(cv_plag_ener,ipcl)  = pcvtile(cv_tile_ener,  itile) *plagvolratio

  plagmomenrm                 = pcvtile(cv_tile_momnrm,itile) *plagvolratio
  pcvplag(cv_plag_xmom,ipcl)  = plagmomenrm*sxn
  pcvplag(cv_plag_ymom,ipcl)  = plagmomenrm*syn
  pcvplag(cv_plag_zmom,ipcl)  = plagmomenrm*szn

  pcvplag(cv_plag_enervapor,ipcl) = 0._rfreal

! ******************************************************************************
! Determine particle position
! include a test to insure particle is in tile cell
!  Note: RFLU defines normals pointing outwards.
! ******************************************************************************

  itercelllocate = 1

199 CONTINUE

  CALL plag_injcsetpoolpos( pregion,nvert,normdirflag,xyzvertex, &
                            fc(xcoord:zcoord),fn(xcoord:zcoord), &
                            cellheight,poolxyz )

  xloc = poolxyz(xcoord)
  yloc = poolxyz(ycoord)
  zloc = poolxyz(zcoord)

  celllocate = rflu_ict_testincell(pregion,xloc,yloc,zloc,icg)

! TEMPORARY
!   PRINT*,' nVert       = ', nVert
!   PRINT*,' icg         = ' ,icg
!   PRINT*,' iTile       = ' ,iTile
!   PRINT*,' normDirFlag = ', normDirFlag
!   PRINT*,' fc          = ', fc(XCOORD:ZCOORD)
!   PRINT*,' fn          = ', fn(XCOORD:ZCOORD)
!   PRINT*,' cellHeight  = ', cellHeight
!   PRINT*,' poolxyz     = ', poolxyz
!   PRINT*,' cellLocate  = ', cellLocate
! END TEMPORARY

  IF (celllocate .EQV. .false.) THEN
    itercelllocate = itercelllocate+1
    IF ( itercelllocate > iter_cell_locate_max ) THEN
       WRITE(stdout,*) ' PLAG_EjectParticle: Unable to create a particle after ',&
                         iter_cell_locate_max, ' iterations'

       CALL errorstop( global,err_plag_cellindex,__line__ )
    END IF ! iterCellLocate
    goto 199
  ENDIF ! cellLocate

  pcvplag(cv_plag_xpos,ipcl) = poolxyz(xcoord)
  pcvplag(cv_plag_ypos,ipcl) = poolxyz(ycoord)
  pcvplag(cv_plag_zpos,ipcl) = poolxyz(zcoord)

! ******************************************************************************
! Set aiv and arv
! ******************************************************************************

  paiv(aiv_plag_pidini,ipcl) = nextidnumber
  paiv(aiv_plag_regini,ipcl) = ireg
  paiv(aiv_plag_icells,ipcl) = icg
  paiv(aiv_plag_status,ipcl) = plag_status_keep
  paiv(aiv_plag_burnstat,ipcl) = burnstat

  parv(arv_plag_spload,ipcl) = pdvtile(dv_tile_spload,itile)

! ******************************************************************************
! Decrease pool variables
! ******************************************************************************

  DO icont = 1, ncont
    icvplagmass = pcvplagmass(icont)
    icvtilemass = pcvtilemass(icont)
    pcvtile(icvtilemass,itile) = pcvtile(icvtilemass,   itile) -  &
                                 pdvtile(dv_tile_spload,itile) *  &
                                 pcvplag(icvplagmass,ipcl)
  END DO ! iCont

  pcvtile(cv_tile_momnrm,itile) = pcvtile(cv_tile_momnrm,itile) - &
                                  pdvtile(dv_tile_spload,itile) * &
                                  plagmomenrm

  pcvtile(cv_tile_ener  ,itile) = pcvtile(cv_tile_ener  ,itile) - &
                                  pdvtile(dv_tile_spload,itile) * &
                                  pcvplag(cv_plag_ener,ipcl)

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction( global )

END SUBROUTINE plag_ejectparticle




! ******************************************************************************
!
! Purpose: Determine coordinates for injected particle on patch surface
!   based on random number generator.
!
! Description: None.
!
! Input:
!   pRegion       Pointer to region data
!   nVert         Number of vertices
!   normDirFlag   Flag setting normal direction
!   xyzVertex     Coordinates of patch
!   faceCentroid  Coordinates of face centroid
!   sNormal       Vector of face normal
!   cellHeight    Height of cell to which patch is attached
!
! Output:
!   posPlag       Coordinates of injected particle
!
! Notes:
!   1. Normals need to be set pointing inward
!   2. normDirFlag is set to -1 for RFLU and +1 for RFLO
!   3. RFLU always set normals pointing outwards of the computational domain
!   4. RFLU numbers vertices in the counterclockwise direction as 1-2-3-4
!
! ******************************************************************************

  SUBROUTINE plag_injcsetpoolpos( pRegion,nVert,normDirFlag,xyzVertex,    &
                                  facecentroid,snormal,cellheight,posplag )
    
    IMPLICIT NONE

! ******************************************************************************
!   Definitions and declarations
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================
    
    INTEGER,      INTENT(IN) :: nvert,normdirflag
       
    REAL(RFREAL), INTENT(IN) :: cellheight
    REAL(RFREAL), DIMENSION(ZCOORD),   INTENT(IN)  :: facecentroid,snormal
    REAL(RFREAL), DIMENSION(ZCOORD,4), INTENT(IN)  :: xyzvertex
    REAL(RFREAL), DIMENSION(ZCOORD),   INTENT(OUT) :: posplag
     
    TYPE(t_region), POINTER  :: pregion

! ==============================================================================
!   Locals
! ==============================================================================

    LOGICAL :: usetriangle1
     
    REAL(RFREAL), PARAMETER :: height_fraction = 1.0e-3_rfreal
    REAL(RFREAL) :: areatriangle1,areatriangle2,xrand,xtrirand,yrand,zrand
    REAL(RFREAL), DIMENSION(ZCOORD) :: snormalinward, v1,v2  

    TYPE(t_global), POINTER :: global
      
! ******************************************************************************
!   Start
! ******************************************************************************

    rcsidentstring = '$RCSfile: PLAG_ModInjection.F90,v $ $Revision: 1.11 $'

    global => pregion%global

    CALL registerfunction(global,'PLAG_InjcSetPoolPos',&
  'PLAG_ModInjection.F90' )

! ****************************************************************************** 
!   Set normals to point inwards 
! ****************************************************************************** 

    snormalinward(xcoord:zcoord) = normdirflag  *snormal(xcoord:zcoord) 

! ****************************************************************************** 
!   Select a random point within the appropriate triangle 
! ****************************************************************************** 

    xtrirand = rand1uniform(pregion%randData)   
    xrand    = rand1uniform(pregion%randData) 
    yrand    = rand1uniform(pregion%randData) 

! ******************************************************************************     
!   Reflect back into triangle
! ******************************************************************************

    IF ( xrand+yrand  > 1.0_rfreal ) THEN    
      xrand = 1.0_rfreal-xrand
      yrand = 1.0_rfreal-yrand
    ENDIF ! xRand
     
    zrand = 1.0_rfreal -(xrand+yrand)

! ******************************************************************************  
!   Setup triangles depending on element type
! ******************************************************************************  

    SELECT CASE(nvert)     

! ==============================================================================
!     For triangular elements, always select triangle1
! ==============================================================================

      CASE (3)
        usetriangle1 = .true.        

! ==============================================================================
!     Partition the face into two triangles and compute the areas 
!     of their projections normal to sNormal
! ==============================================================================

      CASE (4)
        v1(xcoord:zcoord) = xyzvertex(xcoord:zcoord,2) -xyzvertex(xcoord:zcoord,1)
        v2(xcoord:zcoord) = xyzvertex(xcoord:zcoord,3) -xyzvertex(xcoord:zcoord,1)
    
        areatriangle1 = 0.5_rfreal*abs( snormal(xcoord)*(v1(2)*v2(3)-v1(3)*v2(2)) &
                                      + snormal(ycoord)*(v1(3)*v2(1)-v1(1)*v2(3)) &
                                      + snormal(zcoord)*(v1(1)*v2(2)-v1(2)*v2(1)) )  

        v1(xcoord:zcoord) = xyzvertex(xcoord:zcoord,3) -xyzvertex(xcoord:zcoord,1)
        v2(xcoord:zcoord) = xyzvertex(xcoord:zcoord,4) -xyzvertex(xcoord:zcoord,1)

        areatriangle2 = 0.5_rfreal*abs( snormal(xcoord)*(v1(2)*v2(3)-v1(3)*v2(2)) &
                                      + snormal(ycoord)*(v1(3)*v2(1)-v1(1)*v2(3)) &
                                      + snormal(zcoord)*(v1(1)*v2(2)-v1(2)*v2(1)) )

! ------------------------------------------------------------------------------
!     Select which triangle to select a point in 
! ------------------------------------------------------------------------------

        usetriangle1 = ( (areatriangle1+areatriangle2)*xtrirand < areatriangle1 ) 
    END SELECT ! nVert 

! ******************************************************************************  
!   Set particle position on triangular tile
! ****************************************************************************** 
     
    IF ( usetriangle1 ) THEN
      posplag(xcoord:zcoord) = xrand*xyzvertex(xcoord:zcoord,1) &
                             + yrand*xyzvertex(xcoord:zcoord,2) &
                             + zrand*xyzvertex(xcoord:zcoord,3) 
    ELSE
      posplag(xcoord:zcoord) = xrand*xyzvertex(xcoord:zcoord,1) &
                             + yrand*xyzvertex(xcoord:zcoord,3) &
                             + zrand*xyzvertex(xcoord:zcoord,4) 
    ENDIF ! useTriangle1

! ****************************************************************************** 
!   Adjust posPlag to be on tile surface      
! ******************************************************************************     

    posplag(xcoord:zcoord) = posplag(xcoord:zcoord)                   &
                           - dot_product(     posplag(xcoord:zcoord)- &
                                         facecentroid(xcoord:zcoord), &
                                        snormalinward(xcoord:zcoord)) &
                           * snormalinward(xcoord:zcoord) 

! ******************************************************************************     
!   Add tiny offset to put posPlag inside cell
! ******************************************************************************      

    posplag(xcoord:zcoord) = posplag(xcoord:zcoord)      &
                           + height_fraction*cellheight  &
                           *snormalinward(xcoord:zcoord) 

! *****************************************************************************
!   End
! *****************************************************************************

    CALL deregisterfunction(global)

  END SUBROUTINE plag_injcsetpoolpos



!******************************************************************************
!
! Purpose: Invoke ejection model1 to generate particles on tile surfaces.
!
! Description: none.
!
! Input:
!   pRegion      Pointer to region data
!   pPlag        Pointer to particle data
!   pTilePlag    Pointer to tile particle data
!   iTile        Tile identifier
!   icg          Cell index
!   burnStat     Particle burning status
!   nVert        Number of vertices on tile element surface
!   normDirFlag  Normal direction flag
!   xyzVertex    Coordinates of element vertices
!   volMeanPart  Mean particle volume
!   fn           Face normal coordinates
!   fc           Face centroids coordinates
!   cellHeight   Cell height
!
! Output:
!   pPlag      Plag values for cv, aiv, arv of current region.
!   pTilePlag  Tile values for cv, dv of current region.
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE plag_invokeejecmodel1( pRegion,pPlag,pTilePlag,      &
                                  itile,icg,burnstat,nvert,     &
                                  normdirflag,xyzvertex,        &
                                  volmeanpart,fn,fc,cellheight  )

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

  INTEGER, INTENT(IN) :: burnstat,icg,itile,normdirflag,nvert

  REAL(RFREAL),                      INTENT(IN) :: cellheight,volmeanpart
  REAL(RFREAL), DIMENSION(ZCOORD),   INTENT(IN) :: fc,fn
  REAL(RFREAL), DIMENSION(ZCOORD,4), INTENT(IN) :: xyzvertex     

  TYPE(t_region),    POINTER :: pregion  
  TYPE(t_plag),      POINTER :: pplag
  TYPE(t_tile_plag), POINTER :: ptileplag

! ==============================================================================
! Locals
! ==============================================================================

  INTEGER :: injcdiamdist
  INTEGER, POINTER, DIMENSION(:)   :: pcvplagmass, pcvtilemass

  LOGICAL :: injectq

  REAL(RFREAL) :: injcbeta,pi,plagvolratio,poolvolsum,tinjccoeff,tinjcsum
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcvtile,pdvtile

  TYPE(t_global),    POINTER :: global 

! ******************************************************************************
! Start
! ******************************************************************************

  rcsidentstring = '$RCSfile: PLAG_ModInjection.F90,v $ $Revision: 1.11 $'

  global => pregion%global

  CALL registerfunction( global, 'PLAG_InvokeEjecModel1',&
  'PLAG_ModInjection.F90' )

! ******************************************************************************
! Set pointers and variables
! ******************************************************************************

  pcvtile     => ptileplag%cv
  pcvtilemass => ptileplag%cvTileMass
  pdvtile     => ptileplag%dv

  injcdiamdist = pregion%plagInput%injcDiamDist
  injcbeta     = pregion%plagInput%injcBeta

! ******************************************************************************
! Begin injection algorithm
! ******************************************************************************
  
  poolvolsum = sum( pcvtile( pcvtilemass(:),itile ) / &
                     pregion%plagInput%dens(:) )

  pdvtile(dv_tile_poolvold,itile) = poolvolsum

  tinjccoeff = injcbeta *volmeanpart

  tinjcsum   = 0.0_rfreal

  CALL plag_injcsetinjection( pregion,ptileplag,itile,        &
                              tinjccoeff,tinjcsum,poolvolsum, &
                              injectq,plagvolratio            )

! ******************************************************************************
! If injectQ set to True, PLAG datastructure for new injected particle
! ******************************************************************************

  DO WHILE ( injectq )

! TEMPORARY
!   PRINT*,'Ejection Model =', pRegion%plagInput%ejecModel
!   PRINT*,' IFC = ' ,iTile
!   PRINT*,'injectQ     = ',injectQ
!   PRINT*,' plagVolRatio = ',plagVolRatio
!   PRINT*,' tInjcSum   = ',tInjcSum
!   PRINT*,' poolVolSum         = ',poolVolSum
!   PRINT*,' pDvTile(DV_TILE_DIAM,iTile)        = ',pDvTile(DV_TILE_DIAM,iTile)
!   PRINT*,' pDvTile(DV_TILE_COUNTDOWN,iTile)   = ',pDvTile(DV_TILE_COUNTDOWN,iTile)
! END TEMPORARY

! ==============================================================================
!  Eject particle
! ==============================================================================

    CALL plag_ejectparticle( pregion,pplag,ptileplag,       &
                             itile,icg,burnstat,            &
                             nvert,normdirflag,xyzvertex,   &
                             fn,fc,cellheight,plagvolratio  )

! ==============================================================================
!   Update tile diameter and superparticle loading factor
! ==============================================================================

    CALL plag_injcmakeparticle( pregion, injcdiamdist,         &
                                pdvtile(dv_tile_diam,  itile), &
                                pdvtile(dv_tile_spload,itile)  )

! ==============================================================================
!   Check if another particle could be injected
! ==============================================================================

    CALL plag_injcsetinjection( pregion,ptileplag,itile,        &
                                tinjccoeff,tinjcsum,poolvolsum, &
                                injectq,plagvolratio            )

  END DO ! injectQ

! *****************************************************************************
! End
! *****************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_invokeejecmodel1



!******************************************************************************
!
! Purpose: Invoke conservative random ejection model
!          to generate particles on tile surfaces.
!
! Description: none.
!
! Input:
!   pRegion      Pointer to region data
!   pPlag        Pointer to particle data
!   pTilePlag    Pointer to tile particle data
!   iTile        Tile identifier
!   icg          Cell index
!   burnStat     Particle burning status
!   nVert        Number of vertices on tile element surface
!   normDirFlag  Normal direction flag
!   xyzVertex    Coordinates of element vertices
!   volMeanPart  Mean particle volume
!   fn           Face normal coordinates
!   fc           Face centroids coordinates
!   cellHeight   Cell height
!
! Output:
!   pPlag      Plag values for cv, aiv, arv of current region.
!   pTilePlag  Tile values for cv, dv of current region.
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE plag_invokeconsrandejec( pRegion,pPlag,pTilePlag,      &
                                    itile,icg,burnstat,nvert,     &
                                    normdirflag,xyzvertex,        &
                                    volmeanpart,fn,fc,cellheight  )

! ******************************************************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
!   Arguments
! ==============================================================================

  INTEGER, INTENT(IN) :: burnstat,icg,itile,normdirflag,nvert

  REAL(RFREAL),                      INTENT(IN) :: cellheight,volmeanpart
  REAL(RFREAL), DIMENSION(ZCOORD),   INTENT(IN) :: fc,fn
  REAL(RFREAL), DIMENSION(ZCOORD,4), INTENT(IN) :: xyzvertex     

  TYPE(t_region),    POINTER :: pregion  
  TYPE(t_plag),      POINTER :: pplag
  TYPE(t_tile_plag), POINTER :: ptileplag

! ==============================================================================
! Locals
! ==============================================================================

  INTEGER :: ejecmodel,injcdiamdist
  INTEGER, POINTER, DIMENSION(:)   :: pcvplagmass, pcvtilemass

  REAL(RFREAL) :: injcbeta,pi,plagvolratio,poolvolsum,tinjccoeff,tinjcsum
  REAL(RFREAL) :: meansuperparticlevolume, spload
  REAL(RFREAL) :: injcbetafac, injcbetafacinv
  REAL(RFREAL) :: poolvololdsum, poolvolfinal, remainingvolume
  REAL(RFREAL) :: currentsuperparticlevolume, poolvolume
  REAL(RFREAL) :: poolexcess, pexcesss, possibleexcess
  REAL(RFREAL) :: countdown, countdownnext, deltavolume, randunif
  REAL(RFREAL) :: currentparticlevolume
  REAL(RFREAL) :: poolvolcurr
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcvtile,pcvoldtile,pdvtile

  TYPE(t_global),    POINTER :: global 

! ******************************************************************************
! Start
! ******************************************************************************

  rcsidentstring = '$RCSfile: PLAG_ModInjection.F90,v $ $Revision: 1.11 $'

  global => pregion%global

  CALL registerfunction( global, 'PLAG_InvokeConsRandEjec',&
  'PLAG_ModInjection.F90' )

! ******************************************************************************
! Set pointers and variables
! ******************************************************************************

  pcvtile     => ptileplag%cv
  pcvoldtile  => ptileplag%cvOld
  pcvtilemass => ptileplag%cvTileMass
  pdvtile     => ptileplag%dv

 
  injcdiamdist = pregion%plagInput%injcDiamDist
  injcbeta     = pregion%plagInput%injcBeta
  
  pi        = global%pi
  spload       = pregion%plagInput%spload
  meansuperparticlevolume = volmeanpart *spload

  injcbetafac = 2.0_rfreal *injcbeta *meansuperparticlevolume**2
  injcbetafacinv = 1.0_rfreal/injcbetafac

! ******************************************************************************
! Begin injection algorithm
! ******************************************************************************

! ==============================================================================
! poolVolSum is the pool volume after the current timestep is finished
! ==============================================================================

  poolvolsum    = sum( pcvtile( pcvtilemass(:),itile ) / &
                        pregion%plagInput%dens(:) )

! ==============================================================================
! poolVolOldSum is the pool volume before the current timestep began
! ==============================================================================

  poolvololdsum = sum( pcvoldtile( pcvtilemass(:),itile ) / &
                        pregion%plagInput%dens(:) )

! ==============================================================================
! poolVolume is built up incrementally with volume injection,
!  and is decremented when (super)particles are ejected
! ==============================================================================

  poolvolume   = poolvololdsum

! ==============================================================================
! poolVolFinal starts with all volume injection having occurred,
!  and is decremented when (super)particles are ejected
! ==============================================================================

  poolvolfinal = poolvolsum

! ==============================================================================
! poolVolCurr starts with all volume injection having occurred,
!  and is decremented when (super)particles are ejected
! ==============================================================================

  poolvolcurr = poolvolsum

! ==============================================================================
! remainingVolume is the volume to be added to the pool over the
!  current time step: it decreases to 0 as all volume is added
! ==============================================================================

  remainingvolume = poolvolsum -poolvololdsum

! TEMPORARY
!   PRINT*,'Ejection Model =',  pRegion%plagInput%ejecModel
!   PRINT*,' IFC = ' ,iTile
!   PRINT*,' plagVolRatio = ',plagVolRatio
!   PRINT*,' poolVolSum         = ',poolVolSum
!   PRINT*,' remainingVolume    = ',remainingVolume
!   PRINT*,' pDvTile(DV_TILE_DIAM,iTile)        = ',pDvTile(DV_TILE_DIAM,iTile)
!   PRINT*,' pDvTile(DV_TILE_COUNTDOWN,iTile)   = ',pDvTile(DV_TILE_COUNTDOWN,iTile)
! END TEMPORARY

! ==============================================================================
! Start the ejection loop for CRE model
! ==============================================================================

  creejectloop: DO
    currentsuperparticlevolume = pi/6.0_rfreal *spload &
                               * pdvtile(dv_tile_diam,itile)**3.0_rfreal

    currentparticlevolume = pi/6.0_rfreal &
                          * pdvtile(dv_tile_diam,itile)**3.0_rfreal

! ------------------------------------------------------------------------------
!   poolExcess is the volume that would remain in the pool after ejection
! ------------------------------------------------------------------------------

!    poolExcess = poolVolume -currentSuperParticleVolume

! TEMPORARY
    poolexcess = poolvolcurr -currentsuperparticlevolume
! END TEMPORARY

! ------------------------------------------------------------------------------
!   possibleExcess is the volume that would remain in the pool at the end of
!     the timestep after the current particle is (and no others are) ejected
! ------------------------------------------------------------------------------

    possibleexcess = poolexcess + remainingvolume

! ------------------------------------------------------------------------------
!   pExcessS is an auxiliary variable that occurs in a few formulas
! ------------------------------------------------------------------------------

    IF ( poolexcess > 0.0_rfreal ) THEN
      pexcesss = poolexcess**2
    ELSE
      pexcesss = 0.0_rfreal
    ENDIF ! poolExcess

! ------------------------------------------------------------------------------
!   countdown is the internal ejection clock for the current particle:  it
!    starts to tick down when pool volume > current (super)particle volume
! ------------------------------------------------------------------------------
 
    countdown = pdvtile(dv_tile_countdown,itile)

! ------------------------------------------------------------------------------
!   countdownNext is what the clock could tick down to within remainingVolume
! ------------------------------------------------------------------------------

    IF ( possibleexcess > 0.0_rfreal ) THEN
      countdownnext = countdown + injcbetafacinv &
                    * (pexcesss -possibleexcess**2.0_rfreal)
    ELSE
      countdownnext = countdown
    ENDIF ! possibleExcess
 
! ------------------------------------------------------------------------------
!   Pool volume too small to eject particle
! ------------------------------------------------------------------------------

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     if the clock cannot tick down to zero within remainingVolume, no more
!       particles will be ejected in this timestep
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IF ( countdownnext> 0.0_rfreal ) THEN
      poolvolume = poolvolume + remainingvolume ! not needed except as check
      pdvtile(dv_tile_countdown,itile) = countdownnext
      EXIT creejectloop

!-------------------------------------------------------------------------------
!   Pool volume big enough to eject particle
!-------------------------------------------------------------------------------

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     if the clock can tick down below zero, then we need to find when it
!       reaches zero, and then eject the particle
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ELSE

! --- deltaVolume is the volume that must be added to the pool (out of
!     remainingVolume) in order to bring countdown exactly to zero

      deltavolume = sqrt(injcbetafac*countdown +pexcesss) -poolexcess

! --- therefore deltaVolume is added to the pool, the countdown hits zero,
!      and the particle is ejected
      
      poolvolume  = poolvolume +deltavolume -currentsuperparticlevolume

! --- the deltaVolume added to the pool is taken from remainingVolume
      
      remainingvolume = remainingvolume -deltavolume

! --- plagVolRatio must be in terms of poolVolFinal:  the pool cv variables
!     are expressed in terms of values at the end of the time step, not
!     in terms of the intermediate values like poolVolume

      plagvolratio = currentparticlevolume / poolvolfinal

! --- poolVolFinal, like the cv variables, is affected by particle ejection
      
      poolvolfinal = poolvolfinal -currentsuperparticlevolume

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     Eject particle
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      CALL plag_ejectparticle( pregion,pplag,ptileplag,       &
                               itile,icg,burnstat,            &
                               nvert,normdirflag,xyzvertex,   &
                               fn,fc,cellheight,plagvolratio  )

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     Update tile diameter and superparticle loading factor
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      CALL plag_injcmakeparticle( pregion, injcdiamdist,         &
                                  pdvtile(dv_tile_diam,  itile), &
                                  pdvtile(dv_tile_spload,itile)  )

! TEMPORARY
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     Compute current pool volume after being decreased in PLAG_EjectParticle
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
       poolvolcurr  = sum( pcvtile( pcvtilemass(:),itile ) / &
                            pregion%plagInput%dens(:) )

! END TEMPORARY

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     Set countdown for new particle
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      randunif = rand1uniform(pregion%randData)
      
      IF ( randunif <= 0.0_rfreal) THEN
! ----- treats randUnif as 1.9e-22
        pdvtile(dv_tile_countdown,itile) = 50.0_rfreal
      ELSE
        pdvtile(dv_tile_countdown,itile) = -log(randunif)
      END IF ! randUnif

    END IF ! countdownNext
  END DO creejectloop

!------------------------------------------------------------------------------
! Temporary check: poolVolume and poolVolFinal should now agree
!------------------------------------------------------------------------------

  IF ( abs(poolvolume - poolvolfinal) / poolvolfinal > 1.0e-10_rfreal) &
  print*,'CRE check: iTile ', itile,(poolvolume - poolvolfinal) / poolvolfinal  

! *****************************************************************************
! End
! *****************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE plag_invokeconsrandejec

! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE plag_modinjection

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_ModInjection.F90,v $
! Revision 1.11  2008/12/06 08:44:34  mtcampbe
! Updated license.
!
! Revision 1.10  2008/11/19 22:17:46  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.9  2007/03/06 23:15:32  fnajjar
! Renamed nPclsTot to nPclsMax
!
! Revision 1.8  2006/04/07 15:19:23  haselbac
! Removed tabs
!
! Revision 1.7  2005/12/30 16:26:25  fnajjar
! Bug fix for vertex and triangle defs in posPlag
!
! Revision 1.6  2005/12/24 21:35:49  haselbac
! Adapted to changes in ICT, cosmetics
!
! Revision 1.5  2005/01/20 15:36:19  fnajjar
! Bug fix to properly bypass negative pool volume
!
! Revision 1.4  2004/07/28 18:56:54  fnajjar
! Included error trap for memory overflow
!
! Revision 1.3  2004/07/12 15:52:38  fnajjar
! Included interface and routine for conservative random ejection model
!
! Revision 1.2  2004/07/12 14:43:17  fnajjar
! Added pertinent routines to prepare for CRE ejection model
!
! Revision 1.1  2004/03/08 22:35:37  fnajjar
! Initial import of Injection Module
!
! ******************************************************************************