part.f90 - fds-smv - Project Hosting on Google Code

‹r5379

r5420

Source path: svn/ trunk/ FDS/ trunk/ FDS_Source/ part.f90

MODULE PART
 
USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS
USE TRAN
IMPLICIT NONE
PRIVATE
PUBLIC INSERT_DROPLETS_AND_PARTICLES,UPDATE_PARTICLES,REMOVE_DROPLETS,INITIALIZE_DROPLETS,GET_REV_part
CHARACTER(255), PARAMETER :: partid='$Id$'
CHARACTER(255), PARAMETER :: partrev='$Revision$'
CHARACTER(255), PARAMETER :: partdate='$Date$'
REAL(EB) :: INSERT_RATE = 0._EB
INTEGER :: INSERT_COUNT = 0
 
CONTAINS
 

SUBROUTINE INITIALIZE_DROPLETS(NM)

! Insert droplets into the domain at the start of calculation

USE COMP_FUNCTIONS, ONLY : SECOND 
USE PHYSICAL_FUNCTIONS, ONLY : DROPLET_SIZE_DISTRIBUTION 
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS
 
REAL(EB) :: LL,UL,BIN_SIZE,TNOW
INTEGER  :: I,J,IL,IU,IPC
INTEGER, INTENT(IN) :: NM
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC

IF (N_PART==0) RETURN            ! Don't waste time if no particles
IF (EVACUATION_ONLY(NM)) RETURN  ! Don't waste time if an evac mesh
 
TNOW=SECOND()
CALL POINT_TO_MESH(NM)

PART_CLASS_LOOP: DO IPC=1,N_PART

   PC=>PARTICLE_CLASS(IPC)
   
   ! If particles/droplets have a size distribution, initialize here
 
   IF_SIZE_DISTRIBUTION: IF (.NOT.PC%MONODISPERSE .AND. (PC%DIAMETER > 0._EB)) THEN
      CALL DROPLET_SIZE_DISTRIBUTION(PC%DIAMETER,PC%R_CDF(:),PC%CDF(:),NDC,PC%GAMMA,PC%SIGMA)
      BIN_SIZE = PC%R_CDF(NDC)/REAL(NSTRATA,EB)
      STRATIFY: DO I=1,NSTRATA
         LL = (I-1)*BIN_SIZE
         UL =  I   *BIN_SIZE
         LL_LOOP: DO J=1,NDC
            IF (PC%R_CDF(J)>LL) THEN
               IL = J-1 
               PC%IL_CDF(I) = J-1
               EXIT LL_LOOP
            ENDIF
         ENDDO LL_LOOP
         UL_LOOP: DO J=NDC,1,-1
            IF (PC%R_CDF(J)<=UL) THEN
               IU = J 
               PC%IU_CDF(I) = J
               EXIT UL_LOOP
            ENDIF
         ENDDO UL_LOOP
         PC%W_CDF(I) = PC%CDF(IU) - PC%CDF(IL)
      ENDDO STRATIFY
   ENDIF IF_SIZE_DISTRIBUTION

   ! If pacticles/droplets can break up, compute normalized (median = 1) size distribution for child droplets

   IF (PC%BREAKUP .AND. .NOT.PC%MONODISPERSE) THEN
      CALL DROPLET_SIZE_DISTRIBUTION(1._EB,PC%CHILD_R_CDF(:),PC%CHILD_CDF(:),NDC, &
                                     PC%BREAKUP_CHILD_GAMMA,PC%BREAKUP_CHILD_SIGMA)
   ENDIF

ENDDO PART_CLASS_LOOP

TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE INITIALIZE_DROPLETS
 


SUBROUTINE INSERT_DROPLETS_AND_PARTICLES(T,NM)

! Insert sprinkler droplets and lagrangian particles into the domain every DT_INSERT seconds

USE COMP_FUNCTIONS, ONLY : SECOND 
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE GEOMETRY_FUNCTIONS, ONLY: RANDOM_RECTANGLE,RANDOM_CONE
USE DEVICE_VARIABLES
USE CONTROL_VARIABLES
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL(EB) :: PHI_RN,RN,FLOW_RATE,THETA_RN,SPHI,CPHI,MASS_SUM,D_PRES_FACTOR, &
            STHETA,CTHETA,PWT0,DROPLET_SPEED,XI,YJ,ZK,SHIFT1,SHIFT2,XTMP,YTMP,ZTMP,VLEN, &
            TRIGT1,TRIGT2,TNOW,TSI,PIPE_PRESSURE,MASS_PER_TIME,MASS_PER_VOLUME,X1,X2,Y1,Y2,Z1,Z2,BLOCK_VOLUME
REAL(EB), PARAMETER :: VENT_OFFSET=0.1
INTEGER :: I,KS,II,JJ,KK,IC,IL,IU,IPC,DROP_SUM,IIG,JJG,KKG,IW,IBC,IOR,STRATUM,IB
INTEGER :: N_OPEN_NOZZLES,N_INITIAL,N
LOGICAL :: INSERT_ANOTHER_BATCH
TYPE (PROPERTY_TYPE), POINTER :: PY
TYPE (TABLES_TYPE), POINTER :: TA
TYPE (DEVICE_TYPE), POINTER :: DV
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC
TYPE (INITIALIZATION_TYPE), POINTER :: IN
 
IF (EVACUATION_ONLY(NM)) RETURN  ! Don't waste time if an evac mesh
IF (N_PART==0) RETURN            ! Don't waste time if no particles

TNOW=SECOND()
CALL POINT_TO_MESH(NM)

! Count active sprinklers and nozzles

N_OPEN_NOZZLES = 0
N_ACTUATED_SPRINKLERS = 0
COUNT_OPEN_NOZZLES_LOOP: DO KS=1,N_DEVC ! Loop over all devices, but look for sprinklers or nozzles
   DV => DEVICE(KS)
   PY => PROPERTY(DV%PROP_INDEX)
   IF (.NOT. DV%CURRENT_STATE) CYCLE COUNT_OPEN_NOZZLES_LOOP
   IF (PY%PART_ID == 'null')   CYCLE COUNT_OPEN_NOZZLES_LOOP
   N_OPEN_NOZZLES = N_OPEN_NOZZLES + 1
   IF (PY%QUANTITY=='SPRINKLER LINK TEMPERATURE') N_ACTUATED_SPRINKLERS = N_ACTUATED_SPRINKLERS + 1
ENDDO COUNT_OPEN_NOZZLES_LOOP

! Count the number of droplets/particles inserted

INSERT_RATE = 0._EB
INSERT_COUNT = 0

! Add more than one batch of particles/droplets if FDS time step is large enough

OVERALL_INSERT_LOOP: DO  
   
   INSERT_ANOTHER_BATCH = .FALSE.

   ! Loop over all devices, but look for sprinklers or nozzles
      
   SPRINKLER_INSERT_LOOP: DO KS=1,N_DEVC  
      DV => DEVICE(KS)
      PY => PROPERTY(DV%PROP_INDEX)
      IF (PY%PART_ID == 'null')   CYCLE SPRINKLER_INSERT_LOOP
      IF (DV%MESH/=NM)            CYCLE SPRINKLER_INSERT_LOOP
      IF (.NOT. DV%CURRENT_STATE) CYCLE SPRINKLER_INSERT_LOOP
      IPC = PY%PART_INDEX
      PC=>PARTICLE_CLASS(IPC)
      INSERT_RATE = INSERT_RATE + PY%N_INSERT/PY%DT_INSERT
      INSERT_COUNT = INSERT_COUNT + PY%N_INSERT
      IF (DV%T_CHANGE == T) THEN 
         DV%T = T
         CYCLE SPRINKLER_INSERT_LOOP
      ENDIF
      IF (T < PY%PARTICLE_INSERT_CLOCK(NM)) CYCLE SPRINKLER_INSERT_LOOP
   
      ! Determine sprinkler/nozzle flow rate
   
      IF (T_BEGIN == DV%T_CHANGE .AND. PY%FLOW_RAMP_INDEX>=1) THEN
         TSI = T
      ELSE
         TSI = T - DV%T_CHANGE
      ENDIF
   
      IF (PY%PRESSURE_RAMP_INDEX>0) THEN
         PIPE_PRESSURE = EVALUATE_RAMP(REAL(N_OPEN_NOZZLES,EB),0._EB,PY%PRESSURE_RAMP_INDEX)
         D_PRES_FACTOR = (PY%OPERATING_PRESSURE/PIPE_PRESSURE)**(1._EB/3._EB)
         FLOW_RATE = PY%K_FACTOR*SQRT(PIPE_PRESSURE)
      ELSE
         PIPE_PRESSURE = PY%OPERATING_PRESSURE
         D_PRES_FACTOR = 1.0_EB
         FLOW_RATE = PY%FLOW_RATE
      ENDIF      
   
      FLOW_RATE = EVALUATE_RAMP(TSI,PY%FLOW_TAU,PY%FLOW_RAMP_INDEX)*FLOW_RATE*(PC%DENSITY/1000._EB)/60._EB  ! kg/s
   
      IF (FLOW_RATE <= 0._EB) THEN
         DV%T = T
         CYCLE SPRINKLER_INSERT_LOOP
      ENDIF
   
      ! Direction initialization stuff
   
      TRIGT1 = ACOS(-DV%ORIENTATION(3))
      IF (DV%ORIENTATION(2)==0._EB) THEN
         TRIGT2 = ACOS(1._EB)
      ELSE
         TRIGT2 = ACOS(ABS(DV%ORIENTATION(1))/SQRT(DV%ORIENTATION(1)**2+DV%ORIENTATION(2)**2))
      ENDIF
    
      ! Droplet insertion loop
    
      MASS_SUM = 0._EB
      DROP_SUM = 0
   
      DROPLET_INSERT_LOOP: DO I=1,PY%N_INSERT
         IF (NLP+1>MAXIMUM_DROPLETS) EXIT DROPLET_INSERT_LOOP
         NLP = NLP+1
         PARTICLE_TAG = PARTICLE_TAG + NMESHES
         IF (NLP>NLPDIM) THEN
            CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
            DROPLET=>MESHES(NM)%DROPLET
         ENDIF
         DR=>DROPLET(NLP)
    
         ! Set droplet properties
    
         DR%TMP    = PC%TMP_INITIAL         ! Initial temperature
         DR%T      = T                      ! Time of insertion
         DR%IOR    = 0                      ! Orientation index (0 means it is not attached to a solid)
         DR%A_X    = 0._EB                  ! x component of drag on the gas in units of m/s2
         DR%A_Y    = 0._EB                  ! y component of drag on the gas in units of m/s2
         DR%A_Z    = 0._EB                  ! z component of drag on the gas in units of m/s2
         DR%CLASS = IPC                     ! The class of particles
         DR%TAG    = PARTICLE_TAG           ! A unique identifying integer
         IF (MOD(NLP,PC%SAMPLING)==0) THEN  ! Decide whether to show the droplet in Smokeview
            DR%SHOW = .TRUE.    
         ELSE
            DR%SHOW = .FALSE.
         ENDIF
         DR%SPLAT   = .FALSE.
         DR%WALL_INDEX = 0
   
         ! Randomly choose theta and phi
    
         CHOOSE_COORDS: DO
            PICK_PATTERN: IF(PY%SPRAY_PATTERN_INDEX>0) THEN !Use spray pattern table
               TA => TABLES(PY%SPRAY_PATTERN_INDEX)
               CALL RANDOM_NUMBER(RN)
               FIND_ROW: DO II=1,TA%NUMBER_ROWS
                  IF (RN>PY%TABLE_ROW(II)) CYCLE FIND_ROW
                  EXIT FIND_ROW
               END DO FIND_ROW
               CALL RANDOM_NUMBER(RN)
               THETA_RN = TA%TABLE_DATA(II,1) + RN*(TA%TABLE_DATA(II,2)-TA%TABLE_DATA(II,1))
               CALL RANDOM_NUMBER(RN)
               PHI_RN = TA%TABLE_DATA(II,3) + RN*(TA%TABLE_DATA(II,4)-TA%TABLE_DATA(II,3))
               IF (PY%PRESSURE_RAMP_INDEX>0) THEN
                  DROPLET_SPEED = PY%V_FACTOR(II)*SQRT(PIPE_PRESSURE)
               ELSE 
                  DROPLET_SPEED = TA%TABLE_DATA(II,5)
               ENDIF
            ELSE PICK_PATTERN !Use conical spray
               CALL RANDOM_NUMBER(RN)
               THETA_RN = PY%SPRAY_ANGLE(1) + RN*(PY%SPRAY_ANGLE(2)-PY%SPRAY_ANGLE(1))
               CALL RANDOM_NUMBER(RN)
               PHI_RN = RN*TWOPI
               IF (PY%PRESSURE_RAMP_INDEX>0) THEN
                  DROPLET_SPEED = PY%V_FACTOR(1)*SQRT(PIPE_PRESSURE)
               ELSE
                  DROPLET_SPEED = PY%DROPLET_VELOCITY
               ENDIF
            ENDIF PICK_PATTERN
            PHI_RN = PHI_RN + DV%ROTATION  ! Adjust for rotation of head by rotating about z-axis
   
            !  Adjust for tilt of sprinkler pipe
   
            SPHI   = SIN(PHI_RN)
            CPHI   = COS(PHI_RN)
            STHETA = SIN(THETA_RN)
            CTHETA = COS(THETA_RN)
            XTMP   = STHETA*CPHI
            YTMP   = STHETA*SPHI
            ZTMP   = -CTHETA
    
            ! First rotate about y-axis away from x-axis
    
            VLEN   = SQRT(XTMP**2+ZTMP**2)
            SHIFT1 = ACOS(ABS(XTMP)/VLEN)
               SELECT CASE(INT(SIGN(1._EB,ZTMP)))
               CASE (-1)
                  IF (XTMP<0) SHIFT1 = PI-SHIFT1
               CASE ( 1)
               SELECT CASE(INT(SIGN(1._EB,XTMP)))
                  CASE (-1)
                     SHIFT1 = SHIFT1+PI
                  CASE ( 1)
                     SHIFT1 = TWOPI - SHIFT1
               END SELECT
            END SELECT
       
            SHIFT1 = SHIFT1 + TRIGT1
            XTMP = VLEN * COS(SHIFT1)
            ZTMP = -VLEN * SIN(SHIFT1)
    
            ! Second rotate about z-axis away from x-axis
    
            VLEN   = SQRT(XTMP**2+YTMP**2)
            SHIFT1 = ACOS(ABS(XTMP)/VLEN)
            SELECT CASE(INT(SIGN(1._EB,YTMP)))
               CASE ( 1)
                  IF (XTMP<0) SHIFT1 = PI-SHIFT1
               CASE (-1)
               SELECT CASE(INT(SIGN(1._EB,XTMP)))
                  CASE (-1)
                     SHIFT1 = SHIFT1+PI
                  CASE ( 1) 
                     SHIFT1 = TWOPI - SHIFT1
               END SELECT
            END SELECT
    
            SHIFT2 = TRIGT2
            SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(1))))
               CASE (-1)
                  IF (DV%ORIENTATION(2)>0) SHIFT2 = TWOPI - SHIFT2
               CASE ( 1)
               SELECT CASE(INT(SIGN(1._EB,DV%ORIENTATION(2))))
                  CASE (-1) 
                     SHIFT2 = PI-SHIFT2
                  CASE ( 1)
                     SHIFT2 = SHIFT2+ PI
               END SELECT
            END SELECT
            SHIFT1=SHIFT1+SHIFT2
            XTMP = VLEN * COS(SHIFT1)
            YTMP = VLEN * SIN(SHIFT1)
    
            ! Compute initial position and velocity of droplets
    
            DR%U = DROPLET_SPEED*XTMP
            DR%V = DROPLET_SPEED*YTMP
            DR%W = DROPLET_SPEED*ZTMP
            DR%X = DV%X + PY%OFFSET*XTMP
            DR%Y = DV%Y + PY%OFFSET*YTMP
            DR%Z = DV%Z + PY%OFFSET*ZTMP
            IF (TWO_D) THEN
               DR%V = 0._EB
               DR%Y = DV%Y
            ENDIF
            IF (DR%X<=XS .OR. DR%X>=XF) CYCLE CHOOSE_COORDS
            IF (DR%Y<=YS .OR. DR%Y>=YF) CYCLE CHOOSE_COORDS
            IF (DR%Z<=ZS .OR. DR%Z>=ZF) CYCLE CHOOSE_COORDS
            XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
            YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
            ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
            II = FLOOR(XI+1._EB)
            JJ = FLOOR(YJ+1._EB)
            KK = FLOOR(ZK+1._EB)
            IC = CELL_INDEX(II,JJ,KK)
            IF (.NOT.SOLID(IC)) EXIT CHOOSE_COORDS
    
         ENDDO CHOOSE_COORDS
   
         ! Randomly choose droplet size according to Cumulative Distribution Function (CDF)
    
         IF (PC%MONODISPERSE) THEN
            DR%R   = 0.5_EB*PC%DIAMETER*D_PRES_FACTOR
            DR%PWT = 1._EB
         ELSE
            CALL RANDOM_NUMBER(RN)            
            STRATUM = NINT(REAL(NSTRATA,EB)*RN+0.5_EB)
            IL = PC%IL_CDF(STRATUM)
            IU = PC%IU_CDF(STRATUM)
            CALL RANDOM_CHOICE(PC%CDF(IL:IU), PC%R_CDF(IL:IU), IU-IL,DR%R)
            DR%PWT = PC%W_CDF(STRATUM)
            DR%R = D_PRES_FACTOR*DR%R
            IF (2._EB*DR%R > PC%MAXIMUM_DIAMETER) THEN
               DR%PWT = DR%PWT*DR%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
               DR%R = 0.5_EB*PC%MAXIMUM_DIAMETER
            ENDIF
         ENDIF
    
         ! Sum up mass of liquid being introduced
   
         MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
         DROP_SUM = DROP_SUM + 1
         
      ENDDO DROPLET_INSERT_LOOP
    
      ! Compute weighting factor for the droplets just inserted
      
      IF (DROP_SUM > 0) THEN
         PWT0 = FLOW_RATE*PY%DT_INSERT/MASS_SUM
         DROPLET(NLP-DROP_SUM+1:NLP)%PWT = DROPLET(NLP-DROP_SUM+1:NLP)%PWT*PWT0
      ENDIF
      
      ! Indicate that droplets from this device have been inserted at this time T
   
      DV%T = T 
   
   ENDDO SPRINKLER_INSERT_LOOP
   
   ! Loop through all boundary cells and insert particles if appropriate
    
   WALL_INSERT_LOOP: DO IW=1,NWC
    
      IF (T < TW(IW))                        CYCLE WALL_INSERT_LOOP
      IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) CYCLE WALL_INSERT_LOOP
      IBC =  IJKW(5,IW)
      SF  => SURFACE(IBC)
      IPC =  SF%PART_INDEX
      IF (IPC < 1)    CYCLE WALL_INSERT_LOOP
      PC  => PARTICLE_CLASS(IPC)
      IF (PC%DEVC_INDEX>0) THEN
         IF (.NOT.DEVICE(PC%DEVC_INDEX)%CURRENT_STATE) CYCLE WALL_INSERT_LOOP
      ENDIF
      IF (PC%CTRL_INDEX>0) THEN
         IF (.NOT.CONTROL(PC%CTRL_INDEX)%CURRENT_STATE) CYCLE WALL_INSERT_LOOP
      ENDIF
      IF (T < SF%PARTICLE_INSERT_CLOCK(NM)) CYCLE WALL_INSERT_LOOP

      INSERT_RATE = INSERT_RATE + NPPCW(IW)/SF%DT_INSERT
      INSERT_COUNT = INSERT_COUNT + NPPCW(IW)
      IF (UW(IW) >= -0.0001_EB) CYCLE WALL_INSERT_LOOP
    
      II = IJKW(1,IW)
      JJ = IJKW(2,IW)
      KK = IJKW(3,IW)
      IC = CELL_INDEX(II,JJ,KK)
      IF (.NOT.SOLID(IC)) CYCLE WALL_INSERT_LOOP
    
      IF (NM > 1) THEN
         IIG = IJKW(6,IW)
         JJG = IJKW(7,IW)
         KKG = IJKW(8,IW)
         IF (INTERPOLATED_MESH(IIG,JJG,KKG) > 0) CYCLE WALL_INSERT_LOOP
      ENDIF
    
      ! Loop over all particles for the IW-th cell
   
      IOR = IJKW(4,IW)
      MASS_SUM = 0._EB
      PARTICLE_INSERT_LOOP: DO I=1,NPPCW(IW) 
    
         IF (NLP+1 > MAXIMUM_DROPLETS) EXIT PARTICLE_INSERT_LOOP
    
         NLP = NLP+1
         PARTICLE_TAG = PARTICLE_TAG + NMESHES
         IF (NLP > NLPDIM) THEN
            CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
            DROPLET => MESHES(NM)%DROPLET
         ENDIF
         DR=>DROPLET(NLP)
    
         IF (ABS(IOR)==1) THEN
            IF (IOR== 1) DR%X = X(II)   + VENT_OFFSET*DX(II+1)
            IF (IOR==-1) DR%X = X(II-1) - VENT_OFFSET*DX(II-1)
            CALL RANDOM_NUMBER(RN)
            DR%Y = Y(JJ-1) + DY(JJ)*RN
            CALL RANDOM_NUMBER(RN)
            DR%Z = Z(KK-1) + DZ(KK)*RN
         ENDIF
         IF (ABS(IOR)==2) THEN
            IF (IOR== 2) DR%Y = Y(JJ)   + VENT_OFFSET*DY(JJ+1)
            IF (IOR==-2) DR%Y = Y(JJ-1) - VENT_OFFSET*DY(JJ-1)
            CALL RANDOM_NUMBER(RN)
            DR%X = X(II-1) + DX(II)*RN
            CALL RANDOM_NUMBER(RN)
            DR%Z = Z(KK-1) + DZ(KK)*RN
         ENDIF
         IF (ABS(IOR)==3) THEN
            IF (IOR== 3) DR%Z = Z(KK)   + VENT_OFFSET*DZ(KK+1)
            IF (IOR==-3) DR%Z = Z(KK-1) - VENT_OFFSET*DZ(KK-1)
            CALL RANDOM_NUMBER(RN)
            DR%X = X(II-1) + DX(II)*RN
            CALL RANDOM_NUMBER(RN)
            DR%Y = Y(JJ-1) + DY(JJ)*RN
         ENDIF
    
         SELECT CASE(IOR)  ! Give particles an initial velocity (if desired)
            CASE( 1) 
               DR%U = -UW(IW)
            CASE(-1) 
               DR%U =  UW(IW)
            CASE( 2) 
               DR%V = -UW(IW)
            CASE(-2) 
               DR%V =  UW(IW)
            CASE( 3)
               DR%W = -UW(IW)
            CASE(-3) 
               DR%W =  UW(IW)
         END SELECT
    
         ! Save the insertion time (TP) and scalar property (SP) for the particle
    
         DR%A_X    = 0._EB
         DR%A_Y    = 0._EB
         DR%A_Z    = 0._EB
         DR%CLASS  = IPC
         DR%TAG    = PARTICLE_TAG
         IF (MOD(NLP,PC%SAMPLING)==0) THEN
            DR%SHOW = .TRUE.
         ELSE
            DR%SHOW = .FALSE.
         ENDIF
         DR%SPLAT   = .FALSE.
         DR%WALL_INDEX = 0
         DR%R   = 0.0_EB
         DR%TMP = PC%TMP_INITIAL
         DR%T   = T
         DR%PWT = 1._EB
         DR%IOR = 0
         IF (PC%DIAMETER > 0._EB) THEN
            IF (PC%MONODISPERSE) THEN
               DR%R   = 0.5_EB*PC%DIAMETER
               DR%PWT = 1._EB
            ELSE
               CALL RANDOM_NUMBER(RN)            
               STRATUM = NINT(REAL(NSTRATA,EB)*RN+0.5_EB)
               IL = PC%IL_CDF(STRATUM)
               IU = PC%IU_CDF(STRATUM)
               CALL RANDOM_CHOICE(PC%CDF(IL:IU),PC%R_CDF(IL:IU),IU-IL,DR%R)
               DR%PWT = PC%W_CDF(STRATUM)
               IF (2._EB*DR%R > PC%MAXIMUM_DIAMETER) THEN
                  DR%PWT = DR%PWT*DR%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
                  DR%R = 0.5_EB*PC%MAXIMUM_DIAMETER
               ENDIF
            ENDIF
            MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
         ENDIF
   
      ENDDO PARTICLE_INSERT_LOOP

      IF (MASS_SUM > 0._EB) THEN
         IF (SF%PARTICLE_MASS_FLUX > 0._EB) THEN
            DROPLET(NLP-NPPCW(IW)+1:NLP)%PWT = &
            DROPLET(NLP-NPPCW(IW)+1:NLP)%PWT*SF%PARTICLE_MASS_FLUX*AREA_ADJUST(IW)*AW(IW)*SF%DT_INSERT/MASS_SUM
         ENDIF
      ENDIF
   
   ENDDO WALL_INSERT_LOOP
   

   ! Loop over all INIT lines and look for particles inserted within a specified volume

   VOLUME_INSERT_LOOP: DO IB=1,N_INIT

      IN => INITIALIZATION(IB)

      IPC = IN%PART_INDEX
      IF (IPC<1) CYCLE VOLUME_INSERT_LOOP
      PC => PARTICLE_CLASS(IPC)
      IF (T < IN%PARTICLE_INSERT_CLOCK(NM)) CYCLE VOLUME_INSERT_LOOP
 
      N_INITIAL = IN%NUMBER_INITIAL_DROPLETS
      IF (N_INITIAL==0) CYCLE VOLUME_INSERT_LOOP

      MASS_PER_VOLUME = IN%MASS_PER_VOLUME
      MASS_PER_TIME   = IN%MASS_PER_TIME

      SELECT CASE(IN%SHAPE)
         CASE('BLOCK')
            IF (IN%X1>XF .OR. IN%X2<XS .OR. IN%Y1>YF .OR. IN%Y2<YS .OR. IN%Z1>ZF .OR. IN%Z2<ZS) CYCLE VOLUME_INSERT_LOOP
            X1 = MAX(IN%X1,XS) 
            X2 = MIN(IN%X2,XF)
            Y1 = MAX(IN%Y1,YS) 
            Y2 = MIN(IN%Y2,YF)
            Z1 = MAX(IN%Z1,ZS) 
            Z2 = MIN(IN%Z2,ZF)
            BLOCK_VOLUME = (X2-X1)*(Y2-Y1)*(Z2-Z1)
         CASE ('CONE')
      END SELECT
   
      ! Assign properties to the initial droplets/particles
   
      MASS_SUM = 0._EB
      INSERT_PARTICLE_LOOP: DO I=1,N_INITIAL
         NLP = NLP + 1
         PARTICLE_TAG = PARTICLE_TAG + NMESHES
         IF (NLP>NLPDIM) THEN
            CALL RE_ALLOCATE_DROPLETS(1,NM,0,1000)
            DROPLET=>MESHES(NM)%DROPLET
         ENDIF
         DR=>DROPLET(NLP)     
         BLOCK_OUT_LOOP:  DO
        !!  CALL RANDOM_CONE(DR%X,DR%Y,DR%Z,0.5_EB*(X2+X1),0.5_EB*(Y2+Y1),Z1,0.25_EB,0.5_EB)
            CALL RANDOM_RECTANGLE(DR%X,DR%Y,DR%Z,X1,X2,Y1,Y2,Z1,Z2)
            II = CELLSI(FLOOR((DR%X-XS)*RDXINT)) + 1._EB
            JJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT)) + 1._EB
            KK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT)) + 1._EB
            IF (.NOT.SOLID(CELL_INDEX(II,JJ,KK))) EXIT BLOCK_OUT_LOOP
         ENDDO BLOCK_OUT_LOOP
         DR%U   = 0._EB                     ! No initial velocity
         DR%V   = 0._EB
         DR%W   = 0._EB
         DR%TMP = PC%TMP_INITIAL            ! Initial temperature
         DR%T   = T                         ! Insertion time is current time
         DR%R   = 0._EB                     ! Radius is zero unless DIAMETER has been specified
         DR%PWT = 1._EB                     ! Weighting factor is one unless changed below
         DR%IOR = 0                         ! Orientation of solid surface (0 means the droplet/particle is not attached)
         DR%CLASS = IPC                     ! Class identifier
         DR%TAG   = PARTICLE_TAG            ! Unique integer tag
         IF (MOD(NLP,PC%SAMPLING)==0) THEN  ! Decide whether to show or not show in Smokeview
            DR%SHOW = .TRUE.    
         ELSE
            DR%SHOW = .FALSE.    
         ENDIF
         DR%SPLAT   = .FALSE.
         DR%WALL_INDEX = 0
    
         IF (PC%DIAMETER>0._EB) THEN
            IF (PC%MONODISPERSE) THEN
               DR%R   = 0.5_EB*PC%DIAMETER
               DR%PWT = 1._EB
            ELSE
               CALL RANDOM_NUMBER(RN)            
               STRATUM = NINT(REAL(NSTRATA,EB)*RN+0.5_EB)
               IL = PC%IL_CDF(STRATUM)
               IU = PC%IU_CDF(STRATUM)
               CALL RANDOM_CHOICE(PC%CDF(IL:IU),PC%R_CDF(IL:IU),IU-IL,DR%R)
               DR%PWT = PC%W_CDF(STRATUM)
               IF (2._EB*DR%R > PC%MAXIMUM_DIAMETER) THEN
                  DR%PWT = DR%PWT*DR%R**3/(0.5_EB*PC%MAXIMUM_DIAMETER)**3
                  DR%R = 0.5_EB*PC%MAXIMUM_DIAMETER
               ENDIF
            ENDIF
            MASS_SUM = MASS_SUM + DR%PWT*PC%FTPR*DR%R**3
         ENDIF
    
         ! Process special particles that are associated with a particular SURFace type
   
         IF (PC%SURF_INDEX>0) THEN
            SF => SURFACE(PC%SURF_INDEX)
            DR%WALL_INDEX = IN%WALL_INDEX_START + I - 1
            IF (SF%THICKNESS>0._EB) THEN
               DR%R = SF%THICKNESS
            ELSEIF (SF%RADIUS>0._EB) THEN
               DR%R = SF%RADIUS
            ENDIF
            DR%PWT = 1._EB
            XW(DR%WALL_INDEX) = DR%X
            YW(DR%WALL_INDEX) = DR%Y
            ZW(DR%WALL_INDEX) = DR%Z
            IJKW(1,DR%WALL_INDEX) = II
            IJKW(2,DR%WALL_INDEX) = JJ
            IJKW(3,DR%WALL_INDEX) = KK
            IJKW(4,DR%WALL_INDEX) = 0
            IJKW(5,DR%WALL_INDEX) = PC%SURF_INDEX
            IJKW(6,DR%WALL_INDEX) = II
            IJKW(7,DR%WALL_INDEX) = JJ
            IJKW(8,DR%WALL_INDEX) = KK
         ENDIF
   
      ENDDO INSERT_PARTICLE_LOOP
    
      ! Adjust particle weighting factor PWT so that desired MASS_PER_VOLUME is achieved
      
      IF (MASS_PER_TIME>0._EB) MASS_PER_VOLUME = MASS_PER_TIME*IN%DT_INSERT/BLOCK_VOLUME


      IF (MASS_PER_VOLUME>0._EB) DROPLET(NLP-N_INITIAL+1:NLP)%PWT = & 
                                 DROPLET(NLP-N_INITIAL+1:NLP)%PWT*MASS_PER_VOLUME*BLOCK_VOLUME/MASS_SUM

   ENDDO VOLUME_INSERT_LOOP


   ! Reset particle/droplet insertion clocks
 
   DO N=1,N_INIT
      IN => INITIALIZATION(N)
      IF (T >= IN%PARTICLE_INSERT_CLOCK(NM)) IN%PARTICLE_INSERT_CLOCK(NM) = IN%PARTICLE_INSERT_CLOCK(NM) + IN%DT_INSERT
      IF (T >= IN%PARTICLE_INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
   ENDDO

   DO N=1,N_SURF
      SF => SURFACE(N)
      IF (T >= SF%PARTICLE_INSERT_CLOCK(NM)) SF%PARTICLE_INSERT_CLOCK(NM) = SF%PARTICLE_INSERT_CLOCK(NM) + SF%DT_INSERT
      IF (T >= SF%PARTICLE_INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
   ENDDO
   
   DO N=1,N_PROP
      PY => PROPERTY(N)
      IF (T >= PY%PARTICLE_INSERT_CLOCK(NM)) PY%PARTICLE_INSERT_CLOCK(NM) = PY%PARTICLE_INSERT_CLOCK(NM) + PY%DT_INSERT
      IF (T >= PY%PARTICLE_INSERT_CLOCK(NM)) INSERT_ANOTHER_BATCH = .TRUE.
   ENDDO

   IF (.NOT.INSERT_ANOTHER_BATCH) EXIT OVERALL_INSERT_LOOP

ENDDO OVERALL_INSERT_LOOP
 
TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE INSERT_DROPLETS_AND_PARTICLES
 
 

SUBROUTINE RANDOM_CHOICE(CDF,VAR,NPTS,CHOICE)
 
INTEGER,  INTENT(IN)  :: NPTS
REAL(EB), INTENT(IN)  :: CDF(0:NPTS),VAR(0:NPTS)
REAL(EB), INTENT(OUT) :: CHOICE
INTEGER  :: IT
REAL(EB) :: CFRAC,RN,A,B
 
CALL RANDOM_NUMBER(RN)
A = MINVAL(CDF)
B = MAXVAL(CDF)
RN = A + (B-A)*RN
 
CDF_LOOP: DO IT=1,NPTS
   IF (CDF(IT) > RN) THEN
      CFRAC  = (RN-CDF(IT-1))/(CDF(IT)-CDF(IT-1))
      CHOICE = VAR(IT-1) + (VAR(IT)-VAR(IT-1))*CFRAC  
      EXIT CDF_LOOP
   ENDIF
ENDDO CDF_LOOP
 
END SUBROUTINE RANDOM_CHOICE
 
 

SUBROUTINE UPDATE_PARTICLES(T,NM)

USE COMP_FUNCTIONS, ONLY : SECOND  
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
REAL(EB) :: TNOW

! Return if this is an evacuation mesh
IF (EVACUATION_ONLY(NM)) RETURN  

! Zero out the number of the droplets in the "orphanage"; that is, the place to hold droplets transferring from mesh to mesh

MESHES(NM)%OMESH(1:NMESHES)%N_DROP_ORPHANS = 0

! Return if there are no particles in this mesh

IF (MESHES(NM)%NLP==0)   RETURN

! Set the CPU timer and point to the current mesh variables

TNOW=SECOND()
CALL POINT_TO_MESH(NM)

! Zero out the contribution by lagrangian particles to divergence

IF (N_EVAP_INDICES>0 .AND. .NOT.EVACUATION_ONLY(NM) .AND. CORRECTOR) D_LAGRANGIAN = 0._EB

! Move the droplets/particles, then compute mass and energy transfer, then add droplet momentum to gas

IF (CORRECTOR) CALL MOVE_PARTICLES(T,NM)
IF (CORRECTOR) CALL PARTICLE_MASS_ENERGY_TRANSFER(T,NM)
CALL PARTICLE_MOMENTUM_TRANSFER

TUSED(8,NM)=TUSED(8,NM)+SECOND()-TNOW
END SUBROUTINE UPDATE_PARTICLES



SUBROUTINE MOVE_PARTICLES(T,NM)

! Momentum transfer from all particles and droplets
 
USE COMP_FUNCTIONS, ONLY : SECOND  
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP, AFILL2
USE PHYSICAL_FUNCTIONS, ONLY : DRAG
 
REAL(EB) :: RHO_G,RVC,RDS,RDC,QREL,SFAC,UREL,VREL,WREL,TMP_G,RN,THETA_RN, &
            RD,T,C_DRAG,XI,YJ,ZK,MU_AIR,WE_G,DROP_VOL_FRAC,DROP_DEN, &
            DTSP,DTMIN,UBAR,VBAR,WBAR,BFAC,GRVT1,GRVT2,GRVT3,AUREL,AVREL,AWREL,CONST, &
            UVW,DUMMY=0._EB,X_OLD,Y_OLD,Z_OLD,STEP_FRACTION(-3:3),SURFACE_DROPLET_DIAMETER, &
            T_BU_BAG,T_BU_STRIP,B_1,THROHALF,P_UVWMAX,UVWMAX
LOGICAL :: HIT_SOLID
INTEGER :: ICN,I,IIN,JJN,KKN,II,JJ,KK,IIX,JJY,KKZ,IW,N,NITER,IWP1,IWM1,IWP2,IWM2,IWP3,IWM3,IOR_OLD,IC,IOR_FIRST,IML
INTEGER, INTENT(IN) :: NM
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC

SURFACE_DROPLET_DIAMETER = 0.001_EB  ! All droplets adjusted to this size when on solid (m)
THROHALF = (0.5_EB)**(1./3.)
B_1 = 10._EB

GRVT1 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GX)*GVEC(1) 
GRVT2 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GY)*GVEC(2) 
GRVT3 = -EVALUATE_RAMP(T,DUMMY,I_RAMP_GZ)*GVEC(3) 

UVWMAX = 0._EB
P_UVWMAX = UVWMAX

! Loop through all Lagrangian particles and move them one time step

!$OMP PARALLEL SHARED(UVWMAX) FIRSTPRIVATE(GRVT1,GRVT2,GRVT3,P_UVWMAX)
!$OMP DO PRIVATE(I,DR,PC,XI,YJ,ZK,II,JJ,KK,IC,IIX,JJY,KKZ,UBAR,VBAR,WBAR,RVC,RD,RDS,RDC,UREL,VREL,WREL,QREL) &
!$OMP PRIVATE(TMP_G,RHO_G,MU_AIR,C_DRAG,WE_G,T_BU_BAG,T_BU_STRIP,DROP_VOL_FRAC,DROP_DEN) &
!$OMP PRIVATE(SFAC,CONST,BFAC,AUREL,AVREL,AWREL,DTMIN,UVW,NITER,DTSP,N) &
!$OMP PRIVATE(X_OLD,Y_OLD,Z_OLD,IW,RN,THETA_RN,IIN,JJN,KKN,ICN,IOR_OLD,HIT_SOLID) &
!$OMP PRIVATE(IWP1,IWM1,IWP2,IWM2,IWP3,IWM3,STEP_FRACTION,IML,IOR_FIRST) 


DROPLET_LOOP: DO I=1,NLP  

   DR => DROPLET(I)
   PC => PARTICLE_CLASS(DR%CLASS)

   RD  = DR%R
   IF (.NOT. PC%MASSLESS .AND. RD<=0._EB) CYCLE DROPLET_LOOP
   RDS = RD*RD
   RDC = RD*RDS

   ! Determine the current coordinates of the particle

   XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
   YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
   ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
   II  = FLOOR(XI+1._EB)
   JJ  = FLOOR(YJ+1._EB)
   KK  = FLOOR(ZK+1._EB)
   IC  = CELL_INDEX(II,JJ,KK)

   ! Throw out particles that are inside a solid obstruction

   IF (SOLID(IC)) THEN
      DR%X = 1.E6_EB
      CYCLE DROPLET_LOOP
   ENDIF

   DR%A_X    = 0._EB
   DR%A_Y    = 0._EB
   DR%A_Z    = 0._EB

   ! Decide how many time steps to use in tracking particle
   ! Now using initial droplet velocity -> more iterations
    
   DTMIN = DT
   UVW = MAX( ABS(DR%U)*RDX(II),ABS(DR%V)*RDY(JJ),ABS(DR%W)*RDZ(KK) )
   P_UVWMAX = MAX(P_UVWMAX,UVW)
   IF (UVW/=0._EB) DTMIN = MIN(DTMIN,1._EB/UVW)
   NITER = MAX(1,CEILING(DT/DTMIN))
   DTSP  = DT/REAL(NITER,EB)

    
   ! Iterate over a single time step
    
   SUB_TIME_STEP_ITERATIONS: DO N=1,NITER

      ! Get current particle coordinates

      IF (N>1) THEN
         XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
         YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
         ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
         II = FLOOR(XI+1._EB)
         JJ = FLOOR(YJ+1._EB)
         KK = FLOOR(ZK+1._EB)
         IC = CELL_INDEX(II,JJ,KK)
      ENDIF

   ! Interpolate the nearest velocity components of the gas

   IIX  = FLOOR(XI+.5_EB)
   JJY  = FLOOR(YJ+.5_EB)
   KKZ  = FLOOR(ZK+.5_EB)
   UBAR = AFILL2(U,II-1,JJY,KKZ,(DR%X-X(II-1))*RDX(II),YJ-JJY+.5_EB          ,ZK-KKZ+.5_EB)
   VBAR = AFILL2(V,IIX,JJ-1,KKZ,XI-IIX+.5_EB          ,(DR%Y-Y(JJ-1))*RDY(JJ),ZK-KKZ+.5_EB)
   WBAR = AFILL2(W,IIX,JJY,KK-1,XI-IIX+.5_EB          ,YJ-JJY+.5_EB          ,(DR%Z-Z(KK-1))*RDZ(KK))
    
   ! If the particle is just a massless tracer, just move it and go on to the next particle


   IF (PC%MASSLESS) THEN
      DR%U = UBAR
      DR%V = VBAR
      DR%W = WBAR
      DR%X = DR%X + DR%U*DTSP
      DR%Y = DR%Y + DR%V*DTSP
      DR%Z = DR%Z + DR%W*DTSP
      CYCLE DROPLET_LOOP
   ENDIF

   ! Calculate the particle velocity components and the amount of momentum to transfer to the gas

   RVC = RDX(II)*RDY(JJ)*RDZ(KK)

   UREL  = DR%U - UBAR
   VREL  = DR%V - VBAR
   WREL  = DR%W - WBAR
   QREL  = SQRT(UREL*UREL + VREL*VREL + WREL*WREL)
   TMP_G = MAX(TMPMIN,TMP(II,JJ,KK))
   RHO_G = RHO(II,JJ,KK)
   MU_AIR = Y2MU_C(MIN(5000,NINT(TMP_G)))*SPECIES(0)%MW
   DR%RE  = RHO_G*QREL*2._EB*RD/MU_AIR
   DRAG_CALC: IF (DR%IOR==0 .AND. DR%RE>0) THEN
      IF (PC%DRAG_LAW==USER_DRAG) THEN
         C_DRAG = PC%USER_DRAG_COEFFICIENT
      ELSE
         C_DRAG = DRAG(DR%RE,PC%DRAG_LAW)
      ENDIF

      !Secondary break-up
      BREAKUP: IF (PC%BREAKUP) THEN
         WE_G=RHO(II,JJ,KK)*QREL**2*2._EB*RD/PC%SURFACE_TENSION
         T_BU_BAG    = T_END-T_BEGIN
         T_BU_STRIP  = T_END-T_BEGIN
         ! Breakup conditions according to WAVE model by Reitz (1987)
         IF (WE_G >= 12.0_EB)             T_BU_BAG   = 1.72_EB*B_1*SQRT(PC%DENSITY*RDC/(2*PC%SURFACE_TENSION))
         IF (WE_G/SQRT(DR%RE) >= 1.0_EB)  T_BU_STRIP = B_1*(RD/QREL)*SQRT(PC%DENSITY/RHO_G)
         ! droplet age is larger than smallest characteristic breakup time
         IF ((T-DR%T) > MIN(T_BU_BAG,T_BU_STRIP)) THEN
            IF (PC%MONODISPERSE) THEN
               RD       = THROHALF*RD
            ELSE
               DO WHILE (RD >= DR%R)
                  CALL RANDOM_CHOICE(PC%CHILD_CDF(:),PC%BREAKUP_CHILD_DIAMETER*DR%R*PC%CHILD_R_CDF(:),NDC,RD)
               END DO
               RD = MAX(RD,PC%MINIMUM_DIAMETER/2._EB)
            ENDIF
            DR%RE    = RHO_G*QREL*2._EB*RD/MU_AIR
            IF (PC%DRAG_LAW==USER_DRAG) THEN
               C_DRAG = PC%USER_DRAG_COEFFICIENT
            ELSE
               C_DRAG = DRAG(DR%RE,PC%DRAG_LAW)
            ENDIF

            DR%PWT   = DR%PWT*RDC/RD**3
            DR%T     = T
            DR%R     = RD
            RDS      = RD*RD
            RDC      = RD*RDS
         ENDIF
      ENDIF BREAKUP

      ! Drag reduction

      DROP_DEN      = AVG_DROP_DEN(II,JJ,KK,PC%EVAP_INDEX) 
      DROP_VOL_FRAC = DROP_DEN/PC%DENSITY 
      IF (DROP_VOL_FRAC > 1.0E-5_EB) C_DRAG = WAKE_REDUCTION(DROP_VOL_FRAC,DR%RE,C_DRAG)

      IF (.NOT. PC%TREE) THEN
         SFAC    = DR%PWT*RDS*PIO2*QREL*C_DRAG
         DR%A_X  = (REAL(N-1,EB)*DR%A_X + SFAC*UREL*RVC)/REAL(N,EB)
         DR%A_Y  = (REAL(N-1,EB)*DR%A_Y + SFAC*VREL*RVC)/REAL(N,EB)
         DR%A_Z  = (REAL(N-1,EB)*DR%A_Z + SFAC*WREL*RVC)/REAL(N,EB)
      ELSE
         SFAC    = PC%VEG_DRAG_COEFFICIENT*PC%VEG_SV*DR%VEG_PACKING_RATIO*QREL*C_DRAG
         DR%A_X  = SFAC*UREL
         DR%A_Y  = SFAC*VREL
         DR%A_Z  = SFAC*WREL
      ENDIF

      IF (.NOT.PC%STATIC) THEN
         CONST   = 8._EB*PC%DENSITY*RD/(3._EB*RHO_G*C_DRAG*QREL)
         BFAC    = EXP(-DTSP/CONST)
         IF (SPATIAL_GRAVITY_VARIATION) THEN
            GRVT1 = -EVALUATE_RAMP(DR%X,DUMMY,I_RAMP_GX)*GVEC(1) 
            GRVT2 = -EVALUATE_RAMP(DR%X,DUMMY,I_RAMP_GY)*GVEC(2) 
            GRVT3 = -EVALUATE_RAMP(DR%X,DUMMY,I_RAMP_GZ)*GVEC(3) 
         ENDIF
         AUREL   = CONST*GRVT1
         AVREL   = CONST*GRVT2
         AWREL   = CONST*GRVT3
         DR%U    = UBAR + (UREL+AUREL)*BFAC - AUREL
         DR%V    = VBAR + (VREL+AVREL)*BFAC - AVREL
         DR%W    = WBAR + (WREL+AWREL)*BFAC - AWREL
      ENDIF
   ELSE DRAG_CALC ! No drag
      DR%A_X  = 0._EB
      DR%A_Y  = 0._EB
      DR%A_Z  = 0._EB
   ENDIF DRAG_CALC

   ! If the particle does not move, but does drag, go on to the next particle

   IF (PC%STATIC) CYCLE DROPLET_LOOP
    
      ! Update droplet position

      X_OLD = DR%X
      Y_OLD = DR%Y
      Z_OLD = DR%Z
      ! pc = predictor corrector update of X
      DR%X  = DR%X + DR%U*DTSP
      DR%Y  = DR%Y + DR%V*DTSP
      DR%Z  = DR%Z + DR%W*DTSP
    
      ! Droplet hits the floor
    
      IF (POROUS_FLOOR .AND. DR%Z<ZS .AND. PC%EVAP_INDEX>0) THEN
         IC = CELL_INDEX(II,JJ,1)
         IW = WALL_INDEX(IC,-3)
         IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY .AND. ACCUMULATE_WATER .AND. .NOT.DR%SPLAT) THEN
            AWMPUA(IW,PC%EVAP_INDEX) = AWMPUA(IW,PC%EVAP_INDEX) + DR%PWT*PC%FTPR*RDC*RAW(IW)
            DR%SPLAT = .TRUE.
         ENDIF
         CYCLE DROPLET_LOOP
      ENDIF
    
      IF (.NOT.POROUS_FLOOR .AND. DR%Z<ZS) THEN
         IC = CELL_INDEX(II,JJ,1)
         IW = WALL_INDEX(IC,-3)
         IF (BOUNDARY_TYPE(IW)==SOLID_BOUNDARY) THEN
            IF (ACCUMULATE_WATER .AND. .NOT.DR%SPLAT) THEN
               AWMPUA(IW,PC%EVAP_INDEX) = AWMPUA(IW,PC%EVAP_INDEX) + DR%PWT*PC%FTPR*RDC*RAW(IW)
               DR%SPLAT = .TRUE.
            ENDIF
            DR%Z = ZS + 0.05_EB*DZ(1)
            DR%IOR = 3
            CALL RANDOM_NUMBER(RN)
            THETA_RN = TWOPI*RN
            DR%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
            DR%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
            DR%W = 0._EB
         ENDIF
      ENDIF
    
      ! Where is the droplet now?
      XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
      YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
      ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
      IIN = FLOOR(XI+1._EB)
      JJN = FLOOR(YJ+1._EB)
      KKN = FLOOR(ZK+1._EB)
      IF (IIN<0 .OR. IIN>IBP1) CYCLE DROPLET_LOOP
      IF (JJN<0 .OR. JJN>JBP1) CYCLE DROPLET_LOOP
      IF (KKN<0 .OR. KKN>KBP1) CYCLE DROPLET_LOOP
      ICN = CELL_INDEX(IIN,JJN,KKN)
      IF (IC==0 .OR. ICN==0) CYCLE SUB_TIME_STEP_ITERATIONS

      IF (DR%X<XS .AND. BOUNDARY_TYPE(WALL_INDEX(IC,-1))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP
      IF (DR%X>XF .AND. BOUNDARY_TYPE(WALL_INDEX(IC, 1))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP
      IF (DR%Y<YS .AND. BOUNDARY_TYPE(WALL_INDEX(IC,-2))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP
      IF (DR%Y>YF .AND. BOUNDARY_TYPE(WALL_INDEX(IC, 2))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP
      IF (DR%Z<ZS .AND. BOUNDARY_TYPE(WALL_INDEX(IC,-3))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP
      IF (DR%Z>ZF .AND. BOUNDARY_TYPE(WALL_INDEX(IC, 3))/=SOLID_BOUNDARY) CYCLE DROPLET_LOOP

      ! If droplet hits an obstacle, change its properties

      AIR_TO_SOLID: IF (II/=IIN .OR. JJ/=JJN .OR. KK/=KKN) THEN

         IOR_OLD   = DR%IOR
         HIT_SOLID = .FALSE.

         ! Check if any solid boundaries of original grid cell have been crossed

         IWP1 = WALL_INDEX(IC, 1) 
         IWM1 = WALL_INDEX(IC,-1)
         IWP2 = WALL_INDEX(IC, 2)
         IWM2 = WALL_INDEX(IC,-2)
         IWP3 = WALL_INDEX(IC, 3)
         IWM3 = WALL_INDEX(IC,-3)
         STEP_FRACTION = 1._EB

         IF (KKN>KK .AND. BOUNDARY_TYPE(IWP3)==SOLID_BOUNDARY) THEN
            DR%IOR=-3
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KK)-Z_OLD-0.05_EB*DZ(KK))/(DR%Z-Z_OLD))
         ENDIF
         IF (KKN<KK .AND. BOUNDARY_TYPE(IWM3)==SOLID_BOUNDARY) THEN
            DR%IOR= 3
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KKN)-Z_OLD+0.05_EB*DZ(KKN))/(DR%Z-Z_OLD))
         ENDIF
         IF (IIN>II .AND. BOUNDARY_TYPE(IWP1)==SOLID_BOUNDARY) THEN
            DR%IOR=-1
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(II)-X_OLD-0.05_EB*DX(II))/(DR%X-X_OLD))
         ENDIF
         IF (IIN<II .AND. BOUNDARY_TYPE(IWM1)==SOLID_BOUNDARY) THEN
            DR%IOR= 1
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(IIN)-X_OLD+0.05_EB*DX(IIN))/(DR%X-X_OLD))
         ENDIF
         IF (JJN>JJ .AND. BOUNDARY_TYPE(IWP2)==SOLID_BOUNDARY) THEN
            DR%IOR=-2
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJ)-Y_OLD-0.05_EB*DY(JJ))/(DR%Y-Y_OLD))
         ENDIF
         IF (JJN<JJ .AND. BOUNDARY_TYPE(IWM2)==SOLID_BOUNDARY) THEN
            DR%IOR= 2
            HIT_SOLID = .TRUE.
            STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJN)-Y_OLD+0.05_EB*DY(JJN))/(DR%Y-Y_OLD))
         ENDIF


         IF (DR%IOR/=0 .AND. .NOT.ALLOW_SURFACE_DROPLETS) THEN
            DR%R = 0.9_EB*PC%KILL_RADIUS
            CYCLE DROPLET_LOOP
         ENDIF

         IML = MINLOC(STEP_FRACTION,DIM=1)
         IOR_FIRST = 0
         SELECT CASE(IML)
            CASE(1)
               IOR_FIRST = -3
            CASE(2)
               IOR_FIRST = -2
            CASE(3)
               IOR_FIRST = -1
            CASE(5)
               IOR_FIRST =  1
            CASE(6)
               IOR_FIRST =  2
            CASE(7)
               IOR_FIRST =  3
         END SELECT
         DR%WALL_INDEX = WALL_INDEX(IC,-IOR_FIRST)

         ! If no solids boundaries of original cell have been crossed, check boundaries of new grid cell
    
         IF (DR%WALL_INDEX==0) THEN
            IWP1 = WALL_INDEX(ICN, 1)
            IWM1 = WALL_INDEX(ICN,-1)
            IWP2 = WALL_INDEX(ICN, 2)
            IWM2 = WALL_INDEX(ICN,-2)
            IWP3 = WALL_INDEX(ICN, 3)
            IWM3 = WALL_INDEX(ICN,-3)
            HIT_SOLID = .FALSE.
            STEP_FRACTION = 1._EB
            IF (KKN>KK .AND. BOUNDARY_TYPE(IWM3)==SOLID_BOUNDARY) THEN
               DR%IOR=-3
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KK)-Z_OLD-0.05_EB*DZ(KK))/(DR%Z-Z_OLD))
            ENDIF
            IF (KKN<KK .AND. BOUNDARY_TYPE(IWP3)==SOLID_BOUNDARY) THEN
               DR%IOR= 3
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Z(KKN)-Z_OLD+0.05_EB*DZ(KKN))/(DR%Z-Z_OLD))
            ENDIF
            IF (IIN>II .AND. BOUNDARY_TYPE(IWM1)==SOLID_BOUNDARY) THEN
               DR%IOR=-1
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(II)-X_OLD-0.05_EB*DX(II))/(DR%X-X_OLD))
            ENDIF
            IF (IIN<II .AND. BOUNDARY_TYPE(IWP1)==SOLID_BOUNDARY) THEN
               DR%IOR= 1
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(X(IIN)-X_OLD+0.05_EB*DX(IIN))/(DR%X-X_OLD))
            ENDIF
            IF (JJN>JJ .AND. BOUNDARY_TYPE(IWM2)==SOLID_BOUNDARY) THEN
               DR%IOR=-2
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJ)-Y_OLD-0.05_EB*DY(JJ))/(DR%Y-Y_OLD))
            ENDIF
            IF (JJN<JJ .AND. BOUNDARY_TYPE(IWP2)==SOLID_BOUNDARY) THEN
               DR%IOR= 2
               HIT_SOLID = .TRUE.
               STEP_FRACTION(DR%IOR) = MAX(0._EB,(Y(JJN)-Y_OLD+0.05_EB*DY(JJN))/(DR%Y-Y_OLD))
            ENDIF

            IML = MINLOC(STEP_FRACTION,DIM=1)
            IOR_FIRST = 0
            SELECT CASE(IML)
               CASE(1)
                  IOR_FIRST = -3
               CASE(2)
                  IOR_FIRST = -2
               CASE(3)
                  IOR_FIRST = -1
               CASE(5)
                  IOR_FIRST =  1
               CASE(6)
                  IOR_FIRST =  2
               CASE(7)
                  IOR_FIRST =  3
            END SELECT
            DR%WALL_INDEX = WALL_INDEX(ICN,IOR_FIRST)
         ENDIF
   
         ! Check if droplet has crossed no solid planes or too many
  
         IF_HIT_SOLID: IF (HIT_SOLID) THEN
 
            IF (DR%WALL_INDEX==0) CYCLE SUB_TIME_STEP_ITERATIONS

            ! Add droplet mass to accumulated liquid array

            IF (ACCUMULATE_WATER .AND. HIT_SOLID .AND. .NOT.DR%SPLAT) THEN
               AWMPUA(DR%WALL_INDEX,PC%EVAP_INDEX) = AWMPUA(DR%WALL_INDEX,PC%EVAP_INDEX) + DR%PWT*PC%FTPR*DR%R**3*RAW(DR%WALL_INDEX)
               DR%SPLAT = .TRUE.
            ENDIF

            ! Adjust the size of the droplet and weighting factor 

            DR%R   = MIN(0.5_EB*SURFACE_DROPLET_DIAMETER,(DR%PWT*RDC)**ONTH)
            DR%PWT = DR%PWT*RDC/DR%R**3

            ! Move particle to where it almost hits solid

            DR%X = X_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%U
            DR%Y = Y_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%V
            DR%Z = Z_OLD + MINVAL(STEP_FRACTION)*DTSP*DR%W
            XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
            YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
            ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
            IIN = FLOOR(XI+1._EB)
            JJN = FLOOR(YJ+1._EB)
            KKN = FLOOR(ZK+1._EB)
            ICN = CELL_INDEX(IIN,JJN,KKN)
            IF (IOR_OLD==DR%IOR) CYCLE SUB_TIME_STEP_ITERATIONS

         ! Check if droplet has not found surface. Simply remove for now. Todo: search algorith

            IW = WALL_INDEX(ICN, -DR%IOR)
            IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) THEN
               DR%R = 0.9_EB*PC%KILL_RADIUS
               CYCLE DROPLET_LOOP
            ENDIF

            ! Choose a direction for the droplets to move

            DIRECTION: SELECT CASE(DR%IOR)
               CASE (-2:-1,1:2) DIRECTION  
                  DR%U = 0._EB
                  DR%V = 0._EB
                  DR%W = -PC%VERTICAL_VELOCITY 
               CASE (-3) DIRECTION 
                  IF (.NOT.ALLOW_UNDERSIDE_DROPLETS) THEN 
                  DR%U = 0._EB
                  DR%V = 0._EB
                  DR%W = -PC%VERTICAL_VELOCITY 
                  DR%IOR = 0
                  ELSE 
                     CALL RANDOM_NUMBER(RN)
                     THETA_RN = TWOPI*RN
                     DR%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
                     DR%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
                     DR%W = 0._EB
                  ENDIF
               CASE (3) DIRECTION
                  CALL RANDOM_NUMBER(RN)
                  THETA_RN = TWOPI*RN
                  DR%U = PC%HORIZONTAL_VELOCITY*COS(THETA_RN)
                  DR%V = PC%HORIZONTAL_VELOCITY*SIN(THETA_RN)
                  DR%W = 0._EB
            END SELECT DIRECTION

         ENDIF IF_HIT_SOLID
      ENDIF AIR_TO_SOLID 

      ! Check if droplets that were attached to a solid are still attached after the time update

      IW = WALL_INDEX(ICN, -DR%IOR)
      
      SELECT CASE(DR%IOR)
         CASE( 1)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%X = DR%X - 0.2_EB*DX(II)
               DR%W = -DR%W
            ENDIF
         CASE(-1)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%X = DR%X + 0.2_EB*DX(II)
               DR%W = -DR%W
            ENDIF
         CASE( 2)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%Y = DR%Y - 0.2_EB*DY(JJ)
               DR%W = -DR%W
            ENDIF
         CASE(-2)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
               DR%Y = DR%Y + 0.2_EB*DY(JJ)
               DR%W = -DR%W
            ENDIF
         CASE( 3)
            IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN  ! Particle has reached the edge of a horizontal surface
               DR%U = -DR%U
               DR%V = -DR%V
               DR%Z =  DR%Z - 0.2_EB*DZ(KK)
            ENDIF
         CASE(-3)

      END SELECT

      IF (DR%IOR/=0 .AND. BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY) THEN
         DR%IOR = 0
         DR%WALL_INDEX = 0
         IF (RECOUNT_DRIP) DR%SPLAT=.FALSE.
      ELSE
         DR%WALL_INDEX = WALL_INDEX(ICN,-DR%IOR)
      ENDIF

   ENDDO SUB_TIME_STEP_ITERATIONS


ENDDO DROPLET_LOOP
!$OMP END DO NOWAIT
!$OMP CRITICAL
UVWMAX = MAX(UVWMAX,P_UVWMAX)
!$OMP END CRITICAL
!$OMP END PARALLEL 

PART_CFL = DT*UVWMAX

! Remove out-of-bounds particles
 
CALL REMOVE_DROPLETS(T,NM)

CONTAINS

REAL(EB) FUNCTION  WAKE_REDUCTION(DROP_VOL_FRAC,RE,C_DRAG)
REAL(EB)DROP_VOL_FRAC, RE, C_DRAG
REAL(EB) WAKE_VEL, LODM, RELOD
! Compute C_DRAG reduction due to the wake effect
! Ramírez-Muñoz et al. 2007.
IF (DROP_VOL_FRAC <= 0._EB) THEN
   WAKE_REDUCTION = C_DRAG
ELSE
   LODM = (PI/(6._EB*DROP_VOL_FRAC))**(1./3.) - 0.5_EB
   RELOD = RE/(16._EB * LODM)
   WAKE_VEL = 1._EB - 0.5_EB*C_DRAG*(1._EB - EXP(-RELOD))
   WAKE_VEL = MAX(WAKE_VEL,0.15_EB)
   WAKE_REDUCTION = C_DRAG * WAKE_VEL * (1._EB + (RELOD/LODM)*EXP(-RELOD))
ENDIF
RETURN
END FUNCTION WAKE_REDUCTION

END SUBROUTINE MOVE_PARTICLES



SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,NM)
    
! Mass and energy transfer between gas and droplets

USE PHYSICAL_FUNCTIONS, ONLY : GET_MASS_FRACTION,GET_AVERAGE_SPECIFIC_HEAT,GET_MOLECULAR_WEIGHT,GET_SPECIFIC_GAS_CONSTANT,&
GET_SPECIFIC_ENTHALPY
USE COMP_FUNCTIONS, ONLY: SHUTDOWN

REAL(EB), POINTER, DIMENSION(:,:,:) :: DROP_DEN,DROP_RAD,DROP_TMP,MVAP_TOT
REAL(EB), POINTER, DIMENSION(:) :: FILM_THICKNESS
REAL(EB) :: R_DROP,NUSSELT,K_AIR,H_V,H_V_REF, H_L,&
            RVC,WGT,Q_CON_GAS,Q_CON_WALL,Q_RAD,H_HEAT,H_MASS,SH_FAC_GAS,SH_FAC_WALL,NU_FAC_GAS,NU_FAC_WALL, &
            T,PR_AIR,M_VAP,M_VAP_MAX,XI,YJ,ZK,RDT,MU_AIR,H_SOLID,Q_DOT_RAD,DEN_ADD, &
            Y_DROP,Y_GAS,LENGTH,U2,V2,W2,VEL,DENOM,DY_DTMP_DROP,TMP_DROP_NEW,TMP_WALL,H_WALL, &
            SC_AIR,D_AIR,DHOR,SHERWOOD,X_DROP,M_DROP,RHO_G,MW_RATIO,MW_DROP,FTPR,&
            C_DROP,M_GAS,A_DROP,TMP_G,TMP_DROP,TMP_MELT,TMP_BOIL,MINIMUM_FILM_THICKNESS,RE_L,OMRAF,Q_FRAC,Q_TOT,DT_SUBSTEP, &
            CP,H_NEW,YY_GET(1:N_SPECIES),M_GAS_NEW,MW_GAS,CP2,VEL_REL,DELTA_H_G,TMP_G_I,H_G_OLD,H_L_REF,&
            TMP_G_NEW,DT_SUM
INTEGER :: I,II,JJ,KK,IW,IGAS,N_PC,EVAP_INDEX,N_SUBSTEPS,ITMP
INTEGER, INTENT(IN) :: NM
REAL(EB), PARAMETER :: RUN_AVG_FAC=0.5_EB
LOGICAL :: TEMPITER
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC
TYPE (SURFACE_TYPE), POINTER :: SF

! Initializations
RDT    = 1._EB/DT
OMRAF  = 1._EB - RUN_AVG_FAC
FUEL_DROPLET_MLR(NM) = 0._EB

! Rough estimates

MINIMUM_FILM_THICKNESS = 1.E-5_EB   ! Minimum thickness of liquid film on the surface (m)
H_SOLID                = 300._EB    ! Heat transfer coefficient from solid surface to drop (W/m2/K)

! Empirical coefficients

D_AIR                  = 2.6E-5_EB  ! Water Vapor - Air binary diffusion (m2/s at 25 C, Incropera & DeWitt, Table A.8) 
SC_AIR                 = 0.6_EB     ! NU_AIR/D_AIR (Incropera & DeWitt, Chap 7, External Flow)
PR_AIR                 = 0.7_EB     
SH_FAC_GAS             = 0.6_EB*SC_AIR**ONTH
NU_FAC_GAS             = 0.6_EB*PR_AIR**ONTH        
SH_FAC_WALL            = 0.037_EB*SC_AIR**ONTH
NU_FAC_WALL            = 0.037_EB*PR_AIR**ONTH        

! Working arrays

IF (N_EVAP_INDICES>0) THEN
   WCPUA  = RUN_AVG_FAC*WCPUA
   WMPUA  = RUN_AVG_FAC*WMPUA
   DROP_DEN => WORK4
   DROP_RAD => WORK5
   DROP_TMP => WORK6
   MVAP_TOT => WORK7
   MVAP_TOT = 0._EB
ENDIF

! Loop over all types of evaporative species

EVAP_INDEX_LOOP: DO EVAP_INDEX = 1,N_EVAP_INDICES

   FILM_THICKNESS => WALL_WORK2
   FILM_THICKNESS =  0._EB

   ! Loop over all particle/droplet classes that have the given evaporative index
   ! First loop is for evaporation and energy transfer
   
   PART_CLASS_LOOP: DO N_PC = 1,N_PART

      PC => PARTICLE_CLASS(N_PC)

      ! If the particle class is associated with a particular SURF, just get its temperature and cycle
      
      IF (PC%SURF_INDEX>0) THEN
         SF => SURFACE(PC%SURF_INDEX)
         DO I=1,NLP
            DR=>DROPLET(I)
            IF (DR%CLASS==N_PC) THEN
               IW = DR%WALL_INDEX
               DR%TMP = TMP_F(IW)
               IF (SF%SHRINK) THEN
                  DR%R = SUM(WALL(IW)%LAYER_THICKNESS)
               ELSE
                  IF (SF%THERMALLY_THICK) THEN
                     DR%R = MAXVAL(SF%X_S)
                  ELSE
                     DR%R = SF%RADIUS
                  ENDIF
               ENDIF
            ENDIF
         ENDDO
         CYCLE PART_CLASS_LOOP
      ENDIF

      ! Check to see of the particles/droplets evaporate
      
      IF (.NOT.PC%EVAPORATE) CYCLE PART_CLASS_LOOP
      IF (PC%EVAP_INDEX/=EVAP_INDEX) CYCLE PART_CLASS_LOOP
      IF (PC%TREE) CYCLE PART_CLASS_LOOP

      ! Initialize quantities common to the PARTICLE_CLASS 

      FTPR     = PC%FTPR
      TMP_MELT = PC%TMP_MELT
      TMP_BOIL = PC%TMP_V
      IGAS     = PC%SPEC_INDEX
      MW_DROP  = SPECIES(IGAS)%MW
      H_V_REF  = PC%H_V(NINT(PC%H_V_REFERENCE_TEMPERATURE))
      H_L_REF  = PC%C_P_BAR(NINT(TMP_MELT))*TMP_MELT
      ! Loop through all droplets in the class and determine the depth of the liquid film on each surface cell

      FILM_SUMMING_LOOP: DO I=1,NLP
         DR => DROPLET(I)
         IF (DR%IOR==0)        CYCLE FILM_SUMMING_LOOP
         IF (DR%WALL_INDEX==0) CYCLE FILM_SUMMING_LOOP
         IF (DR%CLASS /= N_PC) CYCLE FILM_SUMMING_LOOP
         IF (DR%R<=0._EB)      CYCLE FILM_SUMMING_LOOP
         IW = DR%WALL_INDEX
         FILM_THICKNESS(IW) = FILM_THICKNESS(IW) + DR%PWT*DR%R**3/AW(IW)
      ENDDO FILM_SUMMING_LOOP
      FILM_THICKNESS = FILM_THICKNESS*FTPR/PC%DENSITY

      FILM_THICKNESS = MAX(MINIMUM_FILM_THICKNESS,FILM_THICKNESS) 

      ! Loop through all droplets within the class and determine mass/energy transfer

      DROPLET_LOOP: DO I=1,NLP

         DR => DROPLET(I)
         IF (DR%CLASS /= N_PC) CYCLE DROPLET_LOOP
         IF (DR%R<=0._EB)      CYCLE DROPLET_LOOP

         ! Determine the current coordinates of the particle

         XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
         YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
         ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
         II  = FLOOR(XI+1._EB)
         JJ  = FLOOR(YJ+1._EB)
         KK  = FLOOR(ZK+1._EB)
         RVC = RDX(II)*RDY(JJ)*RDZ(KK)
         
         ! Determine how many sub-time step iterations are needed and then iterate over the time step.
         ! This is not fully functional. Keep as a placeholder for now.

    !!!  N_SUBSTEPS = 0  
    !!!  GET_N: DO ITER=1,100
    !!!     N_SUBSTEPS = N_SUBSTEPS + 1
    !!!     IF (DT/REAL(N_SUBSTEPS,EB)<=10000000.*DR%R) EXIT GET_N
    !!!  ENDDO GET_N

         N_SUBSTEPS = 1
         DT_SUBSTEP = DT/REAL(N_SUBSTEPS,EB) 
         DT_SUM = 0._EB
         TIME_ITERATION_LOOP: DO WHILE (DT_SUM < DT)

            IF (N_SPECIES==0) THEN
               MW_GAS = SPECIES(0)%MW
            ELSE
               YY_GET = YY(II,JJ,KK,:)
               CALL GET_MOLECULAR_WEIGHT(YY_GET,MW_GAS)
            ENDIF
            MW_RATIO = MW_GAS/MW_DROP


            ! Initialize droplet thermophysical data

            R_DROP   = DR%R
            M_DROP   = FTPR*R_DROP**3
            TMP_DROP = DR%TMP
            ITMP     = NINT(TMP_DROP)
            H_V      = PC%H_V(ITMP)
            C_DROP   = PC%C_P(ITMP)
            H_L      = PC%C_P_BAR(ITMP)*TMP_DROP-H_L_REF
            WGT      = DR%PWT
            DHOR     = H_V*MW_DROP/R0   
            ! Gas conditions

            TMP_G  = TMP(II,JJ,KK)
            RHO_G  = RHO(II,JJ,KK)
            MU_AIR = Y2MU_C(MIN(5000,NINT(TMP_G)))*SPECIES(0)%MW
            M_GAS  = RHO_G/RVC        
            M_VAP_MAX = (0.33_EB * M_GAS - MVAP_TOT(II,JJ,KK)) / WGT!limit to avoid diveregence errors
            K_AIR  = CPOPR*MU_AIR
            IF (IGAS>0) THEN
               Y_GAS = YY(II,JJ,KK,IGAS)
            ELSE
               Y_GAS = 0._EB
            ENDIF
            U2 = 0.5_EB*(U(II,JJ,KK)+U(II-1,JJ,KK))
            V2 = 0.5_EB*(V(II,JJ,KK)+V(II,JJ-1,KK))
            W2 = 0.5_EB*(W(II,JJ,KK)+W(II,JJ,KK-1))
            VEL_REL = SQRT((U2-DR%U)**2+(V2-DR%V)**2+(W2-DR%W)**2)
            ! Set variables for heat transfer on solid

            SOLID_OR_GAS_PHASE: IF (DR%IOR/=0 .AND. DR%WALL_INDEX>0) THEN

               IW   = DR%WALL_INDEX
               A_DROP = M_DROP/(FILM_THICKNESS(IW)*PC%DENSITY)
               TMP_WALL = TMP_F(IW) 
               SELECT CASE(ABS(DR%IOR))
                  CASE(1)
                     VEL = SQRT(V2**2+W2**2)
                  CASE(2)
                     VEL = SQRT(U2**2+W2**2)
                  CASE(3)
                     VEL = SQRT(U2**2+V2**2)
               END SELECT
               LENGTH   = 1._EB
               RE_L     = MAX(5.E5_EB,RHO_G*VEL*LENGTH/MU_AIR)
               NUSSELT  = NU_FAC_WALL*RE_L**0.8_EB
               SHERWOOD = SH_FAC_WALL*RE_L**0.8_EB
               H_HEAT   = NUSSELT*K_AIR/LENGTH
               IF (PC%EVAPORATE) THEN
                  H_MASS = SHERWOOD*D_AIR/LENGTH
               ELSE
                  H_MASS = 0._EB
               ENDIF
               H_WALL    = H_SOLID
               Q_DOT_RAD = A_DROP*QRADIN(IW)

            ELSE SOLID_OR_GAS_PHASE

               A_DROP   = 4._EB*PI*R_DROP**2
               NUSSELT  = 2._EB + NU_FAC_GAS*SQRT(DR%RE)
               SHERWOOD = 2._EB + SH_FAC_GAS*SQRT(DR%RE)
               H_HEAT   = NUSSELT *K_AIR/(2._EB*R_DROP)
               IF (PC%EVAPORATE) THEN
                  H_MASS = SHERWOOD*D_AIR/(2._EB*R_DROP)
               ELSE
                  H_MASS = 0._EB
               ENDIF
               H_WALL   = 0._EB
               TMP_WALL = TMPA
               IF (AVG_DROP_DEN(II,JJ,KK,EVAP_INDEX )>0._EB) THEN
                  Q_DOT_RAD = QR_W(II,JJ,KK)/SUM(AVG_DROP_DEN(II,JJ,KK,:))*M_DROP
               ELSE
                  Q_DOT_RAD = 0._EB
               ENDIF

            ENDIF SOLID_OR_GAS_PHASE
         
            ! Compute equilibrium droplet vapor mass fraction, Y_DROP, and its derivative w.r.t. droplet temperature
    
            IF (PC%EVAPORATE) THEN
               X_DROP  = MIN(1._EB,EXP(DHOR*(1._EB/TMP_BOIL-1._EB/TMP_DROP)))
               Y_DROP  = X_DROP/(MW_RATIO + (1._EB-MW_RATIO)*X_DROP)
               IF (TMP_DROP < TMP_BOIL) THEN
                  DY_DTMP_DROP = (MW_RATIO/(X_DROP*(1._EB-MW_RATIO)+MW_RATIO)**2)*DHOR*X_DROP/TMP_DROP**2
               ELSE
                  DY_DTMP_DROP = 0._EB
               ENDIF
               IF (Y_DROP<=Y_GAS) H_MASS = 0._EB
            ELSE
               DY_DTMP_DROP = 0._EB
               Y_DROP       = 0._EB
            ENDIF
            ! Update the droplet temperature semi_implicitly
    
            DENOM = 1._EB + (H_HEAT + H_WALL + H_MASS*RHO_G*H_V*DY_DTMP_DROP)*DT_SUBSTEP*A_DROP/(2._EB*M_DROP*C_DROP) 
            TMP_DROP_NEW = ( TMP_DROP + DT_SUBSTEP*( Q_DOT_RAD + &
                             A_DROP*(H_HEAT*(TMP_G   -0.5_EB*TMP_DROP) + H_WALL*(TMP_WALL-0.5_EB*TMP_DROP) -  &
                             H_MASS*RHO_G*H_V*(Y_DROP-0.5_EB*DY_DTMP_DROP*TMP_DROP-Y_GAS))/(M_DROP*C_DROP)) ) / DENOM
            ! Compute the total amount of heat extracted from the gas, wall and radiative fields
            Q_RAD      = DT_SUBSTEP*Q_DOT_RAD
            Q_CON_GAS  = DT_SUBSTEP*A_DROP*H_HEAT*(TMP_G   -0.5_EB*(TMP_DROP+TMP_DROP_NEW))
            Q_CON_WALL = DT_SUBSTEP*A_DROP*H_WALL*(TMP_WALL-0.5_EB*(TMP_DROP+TMP_DROP_NEW))
            Q_TOT      = Q_RAD+Q_CON_GAS+Q_CON_WALL
            ! Compute the total amount of liquid evaporated
            M_VAP = DT_SUBSTEP*A_DROP*H_MASS*RHO_G*(Y_DROP+0.5_EB*DY_DTMP_DROP*(TMP_DROP_NEW-TMP_DROP)-Y_GAS) 
            M_VAP = MAX(0._EB,MIN(M_VAP,M_DROP,M_VAP_MAX))
            
            ! Evaporate completely small droplets
            IF (PC%EVAPORATE .AND. R_DROP<0.5_EB*PC%MINIMUM_DIAMETER) THEN  
               M_VAP  = M_DROP
               IF (Q_TOT>0._EB) THEN
                  Q_FRAC = M_VAP*H_V/Q_TOT 
                  Q_CON_GAS  = Q_CON_GAS*Q_FRAC
                  Q_CON_WALL = Q_CON_WALL*Q_FRAC
                  Q_RAD      = Q_RAD*Q_FRAC
                  Q_TOT  = Q_RAD+Q_CON_GAS+Q_CON_WALL
               ENDIF
            ENDIF
            IF (M_VAP < M_DROP) TMP_DROP_NEW = TMP_DROP + (Q_TOT - M_VAP * H_V)/(C_DROP * (M_DROP - M_VAP))

            ! If the droplet temperature drops below its freezing point, just reset it

            IF (PC%EVAPORATE .AND. TMP_DROP_NEW<TMP_MELT) TMP_DROP_NEW = TMP_MELT

            ! If the droplet temperature reaches boiling, use only enough energy from gas to vaporize liquid

            IF (PC%EVAPORATE .AND. TMP_DROP_NEW>=TMP_BOIL) THEN  
               M_VAP  = MIN(M_VAP_MAX,M_DROP,M_VAP + (TMP_DROP_NEW - TMP_BOIL)*C_DROP*M_DROP/H_V)
               TMP_DROP_NEW = TMP_BOIL
               IF (Q_TOT>0._EB) THEN
                  Q_FRAC = M_VAP*H_V/Q_TOT
                  Q_CON_GAS  = Q_CON_GAS*Q_FRAC
                  Q_CON_WALL = Q_CON_WALL*Q_FRAC
                  Q_RAD      = Q_RAD*Q_FRAC
                  Q_TOT  = Q_RAD+Q_CON_GAS+Q_CON_WALL
               ENDIF
            ENDIF
            M_DROP = M_DROP - M_VAP
        
            ! Compute surface cooling and density
            IF (DR%IOR/=0 .AND. DR%WALL_INDEX>0) THEN
               WCPUA(IW,EVAP_INDEX) = WCPUA(IW,EVAP_INDEX) + OMRAF*WGT*(Q_RAD+Q_CON_WALL)*RAW(IW)/DT_SUBSTEP
               WMPUA(IW,EVAP_INDEX) = WMPUA(IW,EVAP_INDEX) + OMRAF*WGT*M_DROP*RAW(IW)/REAL(N_SUBSTEPS,EB) 
            ENDIF

            !Update gas temperature and determine new subtimestep
            DELTA_H_G = (H_L + H_V - PC%H_V_CORRECTOR)
            YY_GET(:) = YY(II,JJ,KK,:)
            ITMP = MAX(1,MIN(5000,NINT(TMP_G)))
            CALL GET_AVERAGE_SPECIFIC_HEAT(YY_GET,CP,ITMP)            
            H_G_OLD = M_GAS*CP*TMP_G
            M_GAS_NEW = M_GAS + WGT*M_VAP
            TMP_G_NEW = TMP_G
            H_NEW = H_G_OLD + (DELTA_H_G * M_VAP - Q_CON_GAS) * WGT 
            IF (H_NEW > 0._EB) THEN
               YY_GET = YY_GET * M_GAS/M_GAS_NEW
               YY_GET(IGAS) = YY_GET(IGAS) + WGT*M_VAP/M_GAS_NEW
               TMP_G_I = TMP_G
               TEMPITER = .TRUE.
               ITERATE_TEMP: DO WHILE (TEMPITER)
                  TEMPITER=.FALSE.
                  ITMP = MAX(1,MIN(5000,NINT(TMP_G)))
                  IF (N_SPECIES == 0 ) THEN
                     CP2 = Y2CPBAR_C(ITMP)
                  ELSE
                     CALL GET_AVERAGE_SPECIFIC_HEAT(YY_GET,CP2,ITMP)
                  ENDIF
                  TMP_G_I = TMP_G_I+(H_NEW-CP2*TMP_G_I*M_GAS_NEW)/(CP2*M_GAS_NEW)
                  IF ((TMP_G_NEW-TMP_G_I) > 0.5_EB) TEMPITER = .TRUE.
                  TMP_G_NEW = TMP_G_I
               ENDDO ITERATE_TEMP
            ELSE
               DT_SUBSTEP = DT_SUBSTEP * 0.5_EB
               CYCLE TIME_ITERATION_LOOP
            ENDIF
            IF (ABS(TMP_G_NEW/TMP_G - 1._EB) > 0.2_EB) THEN
               DT_SUBSTEP = DT_SUBSTEP * 0.5_EB            
               IF (DT_SUBSTEP <= 0.00001_EB*DT) CALL SHUTDOWN('Numerical instability in particle energy transport')
               CYCLE TIME_ITERATION_LOOP
            ENDIF
            RHO(II,JJ,KK) = M_GAS_NEW*RVC
            YY(II,JJ,KK,:) = YY_GET(:)
            CALL GET_SPECIFIC_GAS_CONSTANT(YY_GET,RSUM(II,JJ,KK))
            TMP(II,JJ,KK) = TMP_G_NEW
            TMP(II,JJ,KK) = MIN(TMPMAX,MAX(TMPMIN,TMP(II,JJ,KK)))
            YY_GET = 0._EB
            YY_GET(IGAS) = 1._EB
            ITMP = MIN(5000,NINT(TMP_G))
            CALL GET_AVERAGE_SPECIFIC_HEAT(YY_GET,CP2,ITMP)
            DELTA_H_G = (DELTA_H_G - CP2 * TMP_G) 
            
            ! Compute contribution to the divergence
            D_LAGRANGIAN(II,JJ,KK) =  D_LAGRANGIAN(II,JJ,KK) + (MW_RATIO *M_VAP /M_GAS + &
                                      (M_VAP*DELTA_H_G- Q_CON_GAS)/H_G_OLD) * WGT / DT_SUBSTEP 

            ! Add fuel evaporation rate to running counter before adjusting its value
            IF (IGAS>0 .AND. IGAS==I_FUEL) FUEL_DROPLET_MLR(NM) = FUEL_DROPLET_MLR(NM) + WGT*M_VAP/DT_SUBSTEP

            ! Adjust mass of evaporated liquid to account for different Heat of Combustion between fuel droplet and gas
            M_VAP = PC%ADJUST_EVAPORATION*M_VAP


            !Track total mass evaporate in cell
            MVAP_TOT(II,JJ,KK) = MVAP_TOT(II,JJ,KK) + WGT*M_VAP
            ! Update droplet quantities
            DR%R   = (M_DROP/FTPR)**ONTH
            DR%TMP = TMP_DROP_NEW
            ! Get out of the loop if the droplet has evaporated completely
            IF (DR%R<=0._EB) CYCLE DROPLET_LOOP
            DT_SUM = DT_SUM + DT_SUBSTEP
            DT_SUBSTEP = MIN(DT-DT_SUM,DT_SUBSTEP * 1.5_EB)

         ENDDO TIME_ITERATION_LOOP
         
      ENDDO DROPLET_LOOP
      
   ENDDO PART_CLASS_LOOP

   ! Second loop is for summing the part quantities

   PART_CLASS_SUM_LOOP: DO N_PC = 1,N_PART

      PC => PARTICLE_CLASS(N_PC)
      IF (PC%EVAP_INDEX/=EVAP_INDEX) CYCLE PART_CLASS_SUM_LOOP

      DROP_DEN = 0._EB
      DROP_TMP = 0._EB
      DROP_RAD = 0._EB

      DROPLET_LOOP_2: DO I=1,NLP

         DR => DROPLET(I)
         IF (DR%CLASS /= N_PC) CYCLE DROPLET_LOOP_2
         IF (DR%R<=0._EB)      CYCLE DROPLET_LOOP_2

         ! Determine the current coordinates of the particle

         XI = CELLSI(FLOOR((DR%X-XS)*RDXINT))
         YJ = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
         ZK = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
         II  = FLOOR(XI+1._EB)
         JJ  = FLOOR(YJ+1._EB)
         KK  = FLOOR(ZK+1._EB)
         RVC = RDX(II)*RDY(JJ)*RDZ(KK)

         M_DROP   = PC%FTPR*DR%R**3

         ! Assign liquid mass to the cell for airborne drops

         IF (DR%IOR==0) THEN
            DEN_ADD = DR%PWT*M_DROP*RVC
            DROP_DEN(II,JJ,KK) = DROP_DEN(II,JJ,KK) + DEN_ADD
            DROP_TMP(II,JJ,KK) = DROP_TMP(II,JJ,KK) + DEN_ADD*DR%TMP
            DROP_RAD(II,JJ,KK) = DROP_RAD(II,JJ,KK) + DEN_ADD*DR%R
         ENDIF

         ! Compute surface density

         IF (DR%IOR/=0 .AND. DR%WALL_INDEX>0) THEN
            IW = DR%WALL_INDEX
            WMPUA(IW,EVAP_INDEX) = WMPUA(IW,EVAP_INDEX) + OMRAF*DR%PWT*M_DROP*RAW(IW)
         ENDIF

      ENDDO DROPLET_LOOP_2

      ! Compute cumulative quantities for droplet "clouds"
    
      DROP_RAD = DROP_RAD/(DROP_DEN+1.E-10_EB)
      DROP_TMP = DROP_TMP/(DROP_DEN+1.E-10_EB)
    
      AVG_DROP_RAD(:,:,:,EVAP_INDEX ) = (AVG_DROP_DEN(:,:,:,EVAP_INDEX )*AVG_DROP_RAD(:,:,:,EVAP_INDEX )+DROP_DEN*DROP_RAD) &
                                          /(AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + DROP_DEN + 1.0E-10_EB)
      AVG_DROP_TMP(:,:,:,EVAP_INDEX ) = (AVG_DROP_DEN(:,:,:,EVAP_INDEX )*AVG_DROP_TMP(:,:,:,EVAP_INDEX )+DROP_DEN*DROP_TMP) &
                                          /(AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + DROP_DEN + 1.0E-10_EB)
      AVG_DROP_TMP(:,:,:,EVAP_INDEX ) = MAX(TMPM,AVG_DROP_TMP(:,:,:,EVAP_INDEX ))
      AVG_DROP_DEN(:,:,:,EVAP_INDEX ) = RUN_AVG_FAC*AVG_DROP_DEN(:,:,:,EVAP_INDEX ) + OMRAF*DROP_DEN
      WHERE (AVG_DROP_DEN(:,:,:,EVAP_INDEX )<0.0001_EB .AND. DROP_DEN==0._EB) AVG_DROP_DEN(:,:,:,EVAP_INDEX ) = 0.0_EB

   ENDDO PART_CLASS_SUM_LOOP

ENDDO EVAP_INDEX_LOOP

! Remove droplets that have completely evaporated
 
CALL REMOVE_DROPLETS(T,NM)
 
END SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER



SUBROUTINE PARTICLE_MOMENTUM_TRANSFER

! Add droplet momentum as a force term in momentum equation

USE COMP_FUNCTIONS, ONLY : SECOND
REAL(EB), POINTER, DIMENSION(:,:,:) :: FVXS,FVYS,FVZS
REAL(EB) :: FLUXMIN,XI,YJ,ZK
INTEGER :: II,JJ,KK,IIX,JJY,KKZ,I,J,K,IC,IW
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC

FVXS  => WORK1
FVYS  => WORK2
FVZS  => WORK3
!$OMP PARALLEL
!$OMP WORKSHARE
FVXS  = 0._EB
FVYS  = 0._EB
FVZS  = 0._EB
!$OMP END WORKSHARE

!$OMP DO PRIVATE(I,DR,PC,XI,YJ,ZK,II,JJ,KK,IIX,JJY,KKZ,IC,IW)
SUM_MOMENTUM_LOOP: DO I=1,NLP
   DR=>DROPLET(I)
   PC=>PARTICLE_CLASS(DR%CLASS)
   IF (DR%IOR/=0)   CYCLE SUM_MOMENTUM_LOOP
   IF (PC%MASSLESS) CYCLE SUM_MOMENTUM_LOOP
   XI  = CELLSI(FLOOR((DR%X-XS)*RDXINT))
   YJ  = CELLSJ(FLOOR((DR%Y-YS)*RDYINT))
   ZK  = CELLSK(FLOOR((DR%Z-ZS)*RDZINT))
   II  = FLOOR(XI+1._EB)
   JJ  = FLOOR(YJ+1._EB)
   KK  = FLOOR(ZK+1._EB)
   IF (SOLID(CELL_INDEX(II,JJ,KK))) CYCLE SUM_MOMENTUM_LOOP
   IIX = FLOOR(XI+.5_EB)
   JJY = FLOOR(YJ+.5_EB)
   KKZ = FLOOR(ZK+.5_EB)
   IC = CELL_INDEX(IIX,JJ,KK)
   IW = WALL_INDEX(IC,1)
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) THEN
      !$OMP ATOMIC
      FVXS(IIX,JJ,KK) = FVXS(IIX,JJ,KK) - DR%A_X
   ENDIF
   IC = CELL_INDEX(II,JJY,KK)
   IW = WALL_INDEX(IC,2) 
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) THEN
      !$OMP ATOMIC
      FVYS(II,JJY,KK) = FVYS(II,JJY,KK) - DR%A_Y
   ENDIF
   IC = CELL_INDEX(II,JJ,KKZ)
   IW = WALL_INDEX(IC,3) 
   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY) THEN
      !$OMP ATOMIC
      FVZS(II,JJ,KKZ) = FVZS(II,JJ,KKZ) - DR%A_Z
   ENDIF
ENDDO SUM_MOMENTUM_LOOP
!$OMP END DO NOWAIT
 
!$OMP SINGLE
FLUXMIN = -FLUXMAX
!$OMP END SINGLE

!$OMP DO COLLAPSE(3) PRIVATE(K,J,I) 
DO K=0,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         IF (FVXS(I,J,K)>FLUXMAX) FVXS(I,J,K) = FLUXMAX
         IF (FVXS(I,J,K)<FLUXMIN) FVXS(I,J,K) = FLUXMIN
         FVX(I,J,K) = FVX(I,J,K) + FVXS(I,J,K)
         IF (FVYS(I,J,K)>FLUXMAX) FVYS(I,J,K) = FLUXMAX
         IF (FVYS(I,J,K)<FLUXMIN) FVYS(I,J,K) = FLUXMIN
         FVY(I,J,K) = FVY(I,J,K) + FVYS(I,J,K)
         IF (FVZS(I,J,K)>FLUXMAX) FVZS(I,J,K) = FLUXMAX
         IF (FVZS(I,J,K)<FLUXMIN) FVZS(I,J,K) = FLUXMIN
         FVZ(I,J,K) = FVZ(I,J,K) + FVZS(I,J,K)
      ENDDO
   ENDDO
ENDDO
!$OMP END DO NOWAIT
!$OMP END PARALLEL
 
END SUBROUTINE PARTICLE_MOMENTUM_TRANSFER
 
 
 
SUBROUTINE REMOVE_DROPLETS(T,NM)

! Remove droplets that do not lie in any mesh
 
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_DROPLETS 
INTEGER :: IKILL,I,NM,IPC
REAL(EB) :: T,MAX_DROP_LIFE
TYPE (DROPLET_TYPE), POINTER :: DR
TYPE (PARTICLE_CLASS_TYPE), POINTER :: PC

IF (NLP + INSERT_COUNT > MAXIMUM_DROPLETS) THEN
   MAX_DROP_LIFE = REAL(MAXIMUM_DROPLETS,EB) / INSERT_RATE
ELSE
   MAX_DROP_LIFE = 1.E6
ENDIF
IKILL = 0
DROP_LOOP: DO I=1,NLP
   WEED_LOOP: DO
      DR=>DROPLET(I) 
      IPC=DR%CLASS 
      PC=>PARTICLE_CLASS(IPC)
      IF (I>NLP-IKILL) EXIT DROP_LOOP
      IF (DR%R<=PC%KILL_RADIUS .AND. .NOT.PC%MASSLESS) THEN 
         CALL REPLACE
         CYCLE WEED_LOOP
      ENDIF
      IF ((T-DR%T>PC%LIFETIME) .OR. (T-DR%T>MAX_DROP_LIFE .AND. IKILL <= INSERT_COUNT)) THEN
         CALL REPLACE
         CYCLE WEED_LOOP
      ENDIF
      IF (DR%X>XS .AND. DR%X<XF .AND.DR%Y>YS .AND. DR%Y<YF .AND. DR%Z>ZS .AND. DR%Z<ZF)  CYCLE DROP_LOOP
      CALL REPLACE
   ENDDO WEED_LOOP
ENDDO DROP_LOOP
 
NLP = NLP - IKILL
 
CONTAINS
 
SUBROUTINE REPLACE
 
INTEGER :: OM,NOM
TYPE (MESH_TYPE), POINTER :: M
TYPE (OMESH_TYPE), POINTER :: M2


NOM = 0
SEARCH_LOOP: DO OM=1,NMESHES
   IF (NIC(NM,OM)==0) CYCLE SEARCH_LOOP
   IF (EVACUATION_ONLY(OM)) CYCLE SEARCH_LOOP
   M=>MESHES(OM)
   IF (DR%X>M%XS .AND. DR%X<M%XF .AND.  DR%Y>M%YS .AND. DR%Y<M%YF .AND.  DR%Z>M%ZS .AND. DR%Z<M%ZF) THEN
      NOM = OM
      EXIT SEARCH_LOOP
   ENDIF
ENDDO SEARCH_LOOP
 
IF (NOM/=0) THEN
   M2=>MESHES(NM)%OMESH(NOM)
   M2%N_DROP_ORPHANS = M2%N_DROP_ORPHANS + 1
   IF (M2%N_DROP_ORPHANS>M2%N_DROP_ORPHANS_DIM) CALL RE_ALLOCATE_DROPLETS(2,NM,NOM,1000)
   M2%DROPLET(M2%N_DROP_ORPHANS) = DROPLET(I)
ENDIF
 
DROPLET(I) = DROPLET(NLP-IKILL)
IKILL = IKILL + 1
 
END SUBROUTINE REPLACE
 
END SUBROUTINE REMOVE_DROPLETS


 
SUBROUTINE GET_REV_part(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') partrev(INDEX(partrev,':')+1:LEN_TRIM(partrev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') partdate

END SUBROUTINE GET_REV_part
END MODULE PART

Show details Hide details

Change log

r5401 by mcgratta on Dec 28, 2009 Diff

FDS Source: Issue 939. Add iteration loop
in pressure solver to better enforce no-
flux BC.

Go to:

Project members, sign in to write a code review

Older revisions

r5379 by drjfloyd on Dec 24, 2009 Diff

FDS Source: Modify ITMP in part.f90 to
prevent ITMP=0.   Issue #910

r5365 by shostikk on Dec 22, 2009 Diff

FDS Source: Implemented PARTICLE_CFL
and particle drag reduction models.

r5324 by mcgratta on Dec 17, 2009 Diff

FDS Source: Add convective heating for
lagrangian particles tied to a SURF
line, and change D_VAP to
D_LAGRANGIAN.

All revisions of this file

File info

Size: 67817 bytes, 1916 lines

View raw file

File properties

svn:eol-style: LF
svn:keywords: Date Revision Author URL Id