integrate_time_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: CFNEWVALUES                                            !
!                                                                      !
!  Purpose: DES - Calculate the new values of particle velocity,       !
!           position, angular velocity etc                             !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE integrate_time_pic

      USE discretelement
      USE mfix_pic

      IMPLICIT NONE

! call the coloring function like approach
      IF(mppic_solid_stress_snider) THEN
         CALL integrate_time_pic_snider
      ELSE
         CALL integrate_time_pic_garg
      ENDIF

      RETURN
      END SUBROUTINE integrate_time_pic

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: INTEGRATE_TIME_PIC_SNIDER                              !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE integrate_time_pic_snider

      USE param
      USE param1
      USE parallel
      USE physprop
      USE constant
      USE fldvar
      USE discretelement
      USE des_bc
      USE mpi_utility
      USE fldvar
      USE cutcell
      USE mfix_pic
      USE randomno
      USE geometry, only: do_k, no_k
      USE fun_avg
      USE functions

      IMPLICIT NONE

! Local Variables
!---------------------------------------------------------------------//
! Loop counters
      INTEGER :: NP, LC, PC
! Drag coefficient divided by parcel mass
      DOUBLE PRECISION :: DP_BAR
! Restitution coefficient plus one
      DOUBLE PRECISION :: ENp1
! Integrated particle velocity without particle normal stress
      DOUBLE PRECISION :: VEL(3)
! Estimated parcel velocity from normal particle stress
      DOUBLE PRECISION :: DELUP(3)
! Actual parcel velocity from particle normal stress.
      DOUBLE PRECISION :: UPRIMETAU(3)
! Slip velocity between parcel and bulk solids
      DOUBLE PRECISION :: SLIPVEL(3)
! Distance traveled by the parcel (and magnitude)
      DOUBLE PRECISION :: DIST(3), DIST_MAG
! Max velocity magnitude of all parcels
      DOUBLE PRECISION :: MAX_VEL
! Mininum fluid cell dimenstion
      DOUBLE PRECISION :: DXYZ_MIN
! Three times the square root of two.
      DOUBLE PRECISION :: l3SQT2
! Volume fraction used to limit parcel movement in close pack regions.
      DOUBLE PRECISION :: EPg_MFP
! Mean free path of parcels based on volume fraction.
      DOUBLE PRECISION :: MFP
! Solids volume fraction at the parcel
      DOUBLE PRECISION :: EPs
! Loop bound
      INTEGER :: LC_BND
! Gravity normal and magnitude
      DOUBLE PRECISION :: GRAV_NORM(3)
!......................................................................!

! Volume fraction used to limit parcel movement in close pack regions.
      epg_mfp = 0.95d0*ep_star
! Initialize max parcel velocity
      max_vel = -undefined
! Restitution coefficient plus one.
      enp1 = mppic_coeff_en1 + 1.0d0
! Three times the square root of two.
      l3sqt2 = merge(3.0d0, 2.0d0, do_k)*sqrt(2.0d0)
! Dimenion loop bound
      lc_bnd = merge(2,3,no_k)

! Calculate the normal of gravity
      IF(grav_mag > small_number) grav_norm = grav/grav_mag

      pc = 1
      DO np = 1, max_pip

         IF(pc.GT.pip) EXIT
         IF(is_nonexistent(np)) cycle

         pc = pc+1

! If a particle is classified as new, then forces are ignored.
! Classification from new to existing is performed in routine
! des_check_new_particle.f
         IF(is_entering(np))THEN
            fc(np,:) = zero
         ELSE
            fc(np,:) = fc(np,:)/pmass(np) + grav(:)
         ENDIF

! DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
! By comparing MFIX equations and equations in those papers,
! D_p = Beta/(EP_S*RHOP)
! F_gp in drag_fgs.f  = Beta*PVOL/EP_S
! Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
         IF(des_continuum_coupled .AND. mppic_pdrag_implicit)&
            dp_bar = f_gp(np)/pmass(np)

! Solids volume fraction
!         EPs = max(ONE-EP_G(PIJK(NP,4)), 1.0d-4)
         eps = max(one-epg_p(np), 1.0d-4)

! Numerically integrated particle velocity without particle normal
! stress force :: Snider (Eq 38)
         vel(:) = (des_vel_new(np,:) + &
            fc(np,:)*dtsolid)/(1.d0+dp_bar*dtsolid)

! Estimated discrete particle velocity from continuum particle
! normal stress graident :: Snider (Eq 39)
         delup(:) = -(dtsolid*ps_grad(:,np)) / &
            (eps*ro_s0(1)*(1.d0+dp_bar*dtsolid))

! Slip velocity between parcel and average solids velocity adjusted
! for gravity driven flows.
            slipvel = avg_vel_wgrav(np,vel) - vel

! Calculate particle velocity from particle normal stress.
         DO lc = 1, lc_bnd
! Snider (40)
            IF(ps_grad(lc,np) <= zero) THEN
               uprimetau(lc) = min(delup(lc), enp1*slipvel(lc))
               uprimetau(lc) = max(uprimetau(lc), zero)
! Snider (41)
            ELSE
               uprimetau(lc) = max(delup(lc), enp1*slipvel(lc))
               uprimetau(lc) = min(uprimetau(lc), zero)
            ENDIF
         ENDDO

! Total particle velocity is the sum of the particle velocity from
! normal stress and velocity from all other foces :: Sinder (Eq 37)
         des_vel_new(np,:) = vel + uprimetau
! Total distance traveled over the current time step
         dist(:) = des_vel_new(np,:)*dtsolid

! Limit parcel movement in close pack regions to the mean free path. 
         IF(epg_p(np) < epg_mfp) THEN
            IF(dot_product(des_vel_new(np,:),ps_grad(:,np)) > zero) THEN
! Total distance traveled given current velocity.
               dist_mag = dot_product(dist, dist)
! Mean free path based on volume fraction :: Snider (43)
               mfp = 3.7d0*des_radius(np)/(l3sqt2*eps)
! Impose the distance limit and calculate the corresponding velocity.
               IF(mfp**2 < dist_mag) THEN
                  dist = mfp*dist/sqrt(dist_mag)
                  des_vel_new(np,:) = dist/dtsolid
               ENDIF
            ENDIF
         ENDIF


! Update the parcel position.
         des_pos_new(np,:) = des_pos_new(np,:) + dist
! Clear out the force array.
         fc(np,:) = zero

! Determine the maximum distance traveled by any single parcel.
         dist_mag = dot_product(des_vel_new(np,:),des_vel_new(np,:)) 
         max_vel = max(max_vel, dist_mag)

      ENDDO

! Get the minimum fluid cell dimension
      dxyz_min = min(minval(dx(imin1:imax1)),minval(dy(jmin1:jmax1)))
      IF(do_k) dxyz_min = min(dxyz_min,minval(dz(kmin1:kmax1)))

! Calculate the maximum time step to prevent a parcel from moving more
! than a fluid cell in one time step.  The global_all_max could get
! expensive in large MPI applications.
      dtpic_cfl = cfl_pic*dxyz_min/(max_vel+small_number)
      CALL global_all_max(dtpic_cfl)

      RETURN
      contains


!......................................................................!
!                                                                      !
! Function: AVG_SOLIDS_vGRAV                                           !
! Author: J.Musser                                                     !
!                                                                      !
! Purpose: Return the average solids velocity around the parcel with   !
! a correction for gravity dominated flows.                            !
!                                                                      !
! Notes: The average solids velocity is reduced by 1/2 when it is      !
! moving within 5 degrees of parall of the direction of gravity. This  !
! is helps prevent overpacking while solids are settling. The value is !
! blended from AVG to 0.5AVG from 0.1 to 0.0 radians.                  !
!                                                                      !
!''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''!
      FUNCTION avg_vel_wgrav(lNP,lVEL) RESULT(aVEL)

      IMPLICIT NONE

! Particle index
      INTEGER, INTENT(IN) :: lNP
! Intermediate parcel velocity
      DOUBLE PRECISION, INTENT(IN) :: lVEL(3)
! Average solids velocity modified for gravity
      DOUBLE PRECISION :: aVEL(3)
! Temp double precision variables.
      DOUBLE PRECISION :: tDP1, tDP2

      avel = avgsolvel_p(:,lnp)

! Reduce the avreage solids velocity for gravity driven flows.
      IF(grav_mag > small_number) THEN
         tdp1 = dot_product(grav_norm,avel)
         IF(tdp1 > zero) THEN
            IF(dot_product(avel,lvel) > zero) THEN
               tdp2 = dot_product(avel, avel)
               IF(tdp1**2 > 0.99*tdp2) avel = &
                  (50.75d0- 50.0d0*tdp1/sqrt(tdp2))*avel   ! 3/4

!               IF(tDP1**2 > 0.99*tDP2) aVEL = aVEL * &
!                  (100.50d0-100.0d0*tDP1/sqrt(tDP2))*aVEL  ! 1/2

!               IF(tDP1**2 > 0.99*tDP2) aVEL = 0.75d0* aVEL
            ENDIF
         ENDIF
      ENDIF

      RETURN
      END FUNCTION avg_vel_wgrav

      END SUBROUTINE integrate_time_pic_snider






























      SUBROUTINE integrate_time_pic_garg

      USE param
      USE param1
      USE parallel
      USE physprop
      USE constant
      USE fldvar
      USE discretelement
      USE des_bc
      USE mpi_utility
      USE mfix_pic
      USE randomno
      USE cutcell
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
      INTEGER NP, M, LC
      INTEGER I, J, K, IJK, IJK_OLD

      DOUBLE PRECISION DD(3), DIST, &
                       DP_BAR, MEANVEL(3), D_GRIDUNITS(3)

      DOUBLE PRECISION MEAN_FREE_PATH
! index to track accounted for particles
      INTEGER PC

! Logical for local debug warnings
      LOGICAL DES_LOC_DEBUG

! maximum distance particles can move in MPPIC
      DOUBLE PRECISION UPRIMEMOD

! dt's in each direction  based on cfl_pic for the mppic case

      DOUBLE PRECISION :: DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ
      DOUBLE PRECISION :: THREEINTOSQRT2, RAD_EFF, POS_Z
      DOUBLE PRECISION :: VELP_INT(3)

      LOGICAL :: DELETE_PART, INSIDE_DOMAIN
      INTEGER :: PIP_DEL_COUNT

      INTEGER :: LPIP_DEL_COUNT_ALL(0:numpes-1), PIJK_OLD(5)

      double precision  sig_u, mean_u
      DOUBLE PRECISION :: DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
      double precision, allocatable, dimension(:,:) ::  rand_vel

      double precision :: norm1, norm2, norm3
      Logical :: OUTER_STABILITY_COND, DES_FIXED_BED
!-----------------------------------------------

      outer_stability_cond = .false.
      des_fixed_bed = .false.
      pc = 1
      dtpic_cfl = large_number
      dtpic_min_x =  large_number
      dtpic_min_y =  large_number
      dtpic_min_z =  large_number

      if(no_k) threeintosqrt2 = 2.d0*sqrt(2.d0)
      if(do_k) threeintosqrt2 = 3.d0*sqrt(2.d0)
      threeintosqrt2 = 3.d0*sqrt(2.d0)
      pip_del_count = 0

      !EPG_MIN2 = MINVAL(EP_G(:))
      !epg_min_loc = MINLOC(EP_G)
      !IJK = epg_min_loc(1)
      allocate(rand_vel(3, max_pip))
      do lc = 1, merge(2,3,no_k)
         mean_u = zero
         sig_u = 1.d0
         CALL nor_rno(rand_vel(lc, 1:max_pip), mean_u, sig_u)
      enddo
      !WRITE(*, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(3,:)), MAXVAL(FC(3,:))
      !WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(3,:)), MAXVAL(FC(3,:))

      DO np = 1, max_pip
         delete_part = .false.
! pradeep skip ghost particles
         if(pc.gt.pip) exit
         if(is_nonexistent(np)) cycle
         pc = pc+1
         if(is_ghost(np) .or. is_entering_ghost(np) .or. is_exiting_ghost(np)) cycle

         des_loc_debug = .false.

         IF(.NOT.is_entering(np))THEN
            fc(np,:) = fc(np,:)/pmass(np) + grav(:)
         ELSE
            fc(np,:) = zero
         ENDIF

         IF(des_fixed_bed) fc(np,:) = zero

         !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
         !By comparing the MFIX and equations in those papers,
         !D_p = Beta/(EP_S*RHOP)
         !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
         !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
         dp_bar = f_gp(np)/(pmass(np))
         IF(.NOT.des_continuum_coupled) dp_bar = zero

         if(.not.mppic_pdrag_implicit) dp_bar = zero
         m = pijk(np,5)
         ijk = pijk(np,4)
         ijk_old = ijk
         pijk_old(:) = pijk(np,:)

         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)

         des_vel_new(np,:) = (des_vel_new(np,:) + &
         & fc(np,:)*dtsolid)/(1.d0+dp_bar*dtsolid)

         IF(des_fixed_bed) des_vel_new(np,:) = zero

         !MPPIC_VPTAU(NP,:) = DES_VEL_NEW(NP,:)

         velp_int(:) = des_vel_new(np,:)

         meanvel(1) = u_s(ijk_old,m)
         meanvel(2) = v_s(ijk_old,m)
         IF(do_k) meanvel(3) = w_s(ijk_old,m)

         rad_eff = des_radius(np)
         !RAD_EFF = (DES_STAT_WT(NP)**(1.d0/3.d0))*DES_RADIUS(NP)
         IF(.not.des_oneway_coupled) then
            mean_free_path  = max(1.d0/(1.d0-ep_star), 1.d0/(1.d0-ep_g(ijk_old)))
            mean_free_path  = mean_free_path*rad_eff/threeintosqrt2
         endif

         DO lc = 1, merge(2,3,no_k)
            !SIG_U = 0.05D0*MEANVEL(LC)
            !DES_VEL_NEW(NP,LC) = DES_VEL_NEW(NP,LC) + SIG_U*RAND_VEL(LC, L )
            !PART_TAUP = RO_Sol(NP)*((2.d0*DES_RADIUS(NP))**2.d0)/(18.d0* MU_G(IJK))
            sig_u = 0.005d0      !*MEANVEL(LC)
            rand_vel(lc, np)  = sig_u*rand_vel(lc, np)*des_vel_new(np,lc)
            IF(des_fixed_bed) rand_vel(lc,np) = zero
            !rand_vel(LC, NP)  = sig_u*rand_vel(LC, NP)/part_taup
            !rand_vel(LC, NP)  = sig_u* mean_free_path*rand_vel(LC, NP)/part_taup
            des_vel_new(np,lc) = des_vel_new(np,lc) + rand_vel(lc, np)
         enddo

         IF(.not.des_oneway_coupled.and.(.not.des_fixed_bed)) CALL mppic_apply_ps_grad_part(np)

         uprimemod = sqrt(dot_product(des_vel_new(np,1:), des_vel_new(np,1:)))

         !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
         !   DES_VEL_NEW(NP,:) = (DES_VEL_NEW(NP,:)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
         !ENDIF

         IF(des_fixed_bed) THEN
            dd(:) = zero
         ELSE
            dd(:) = des_vel_new(np,:)*dtsolid !+ rand_vel(:, NP)*dtsolid
         ENDIF

         CALL pic_find_new_cell(np)

         ijk = pijk(np,4)

         !IF((EP_G(IJK).LT.0.35.and.fluid_at(ijk)).or.(ep_g(ijk_old).lt.0.35)) then !.and.(ijk.ne.ijk_old)) then
         !IF((EP_G(IJK).LT.EP_STAR.and.fluid_at(ijk)).and.(ijk.ne.ijk_old)) then
         inside_domain = .true.
         inside_domain = fluid_at(ijk)

         IF(cut_cell_at(ijk)) THEN
            pos_z = zero
            IF(do_k) pos_z = des_pos_new(np,3)
            CALL get_del_h_des(ijk,'SCALAR', &
            & des_pos_new(np,1),  des_pos_new(np,2), &
            & pos_z, &
            & dist, norm1, norm2, norm3, .true.)

            IF(dist.LE.zero) inside_domain = .false.
         ENDIF

         !IF((EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN)) then
         !IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN) then
         IF(1.d0 - ep_g(ijk).GT. 1.3d0*(1.d0 - ep_star).and.inside_domain.and.ijk.NE.ijk_old.and.outer_stability_cond) then
            dd(:) = zero
            des_vel_new(np,:) = 0.8d0*des_vel_new(np,:)
         ENDIF

         pijk(np,:) = pijk_old(:)

         des_pos_new(np,:) = des_pos_new(np,:) + dd(:)

         dist = sqrt(dot_product(dd,dd))

         d_gridunits(1) = abs(dd(1))/dx(pijk(np,1))
         d_gridunits(2) = abs(dd(2))/dy(pijk(np,2))
         d_gridunits(3) = 1.d0
         IF(do_k) d_gridunits(3) = abs(dd(3))/dz(pijk(np,3))

         dist = sqrt(dot_product(dd,dd))
         dtpic_tmpx = (cfl_pic*dx(pijk(np,1)))/(abs(des_vel_new(np,1))+small_number)
         dtpic_tmpy = (cfl_pic*dy(pijk(np,2)))/(abs(des_vel_new(np,2))+small_number)
         dtpic_tmpz = large_number
         IF(do_k) dtpic_tmpz = (cfl_pic*dz(pijk(np,3)))/(abs(des_vel_new(np,3))+small_number)


! Check if the particle has moved a distance greater than or equal to grid spacing
! if so, then exit

         DO lc = 1, merge(2,3,no_k)
            IF(d_gridunits(lc).GT.one) THEN
               IF(dmp_log.OR.mype.eq.pe_io) THEN

                  WRITE(unit_log, 2001) np, d_gridunits(:), des_vel_new(np,:)
                  WRITE(unit_log, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,np)
                  IF (do_old) WRITE(unit_log, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(np,:)
                  WRITE(unit_log, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(np,:)
                  WRITE(unit_log, '(A,2x,3(g17.8))') 'FC          = ', fc(np,:)

                  WRITE(*, 2001) np, d_gridunits(:), des_vel_new(np,:)

                  WRITE(*, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,np)
                  IF (do_old) WRITE(*, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(np,:)
                  WRITE(*, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(np,:)
                  WRITE(*, '(A,2x,3(g17.8))') 'FC          = ', fc(np,:)
!                  read(*,*)
                  delete_part = .true.

               ENDIF
               !CALL mfix_exit(myPE)
            END IF

         END DO

         IF(.not.delete_part) then
            dtpic_cfl = min(dtpic_cfl, dtpic_tmpx, dtpic_tmpy, dtpic_tmpz)
            dtpic_min_x = min(dtpic_min_x, dtpic_tmpx)
            dtpic_min_y = min(dtpic_min_y, dtpic_tmpy)
            dtpic_min_z = min(dtpic_min_z, dtpic_tmpz)
         ENDIF
         fc(np,:) = zero

         IF(delete_part) THEN
            CALL set_nonexistent(np)
            pip_del_count = pip_del_count + 1
         ENDIF
         IF (des_loc_debug) WRITE(*,1001)
      ENDDO


      DEALLOCATE(rand_vel)

      CALL global_all_max(dtpic_cfl)
      pip = pip - pip_del_count

      lpip_del_count_all(:) = 0
      lpip_del_count_all(mype) = pip_del_count
      CALL global_all_sum(lpip_del_count_all)
      IF((dmp_log).AND.sum(lpip_del_count_all(:)).GT.0) THEN
         IF(print_des_screen) THEN
            WRITE(*,'(/,2x,A,2x,i10,/,A)') &
            'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ',&
             sum(lpip_del_count_all(:)), &
            'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
         ENDIF

         WRITE(unit_log,'(/,2x,A,2x,i10,/,A)') &
            'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC=',&
            sum(lpip_del_count_all(:)), &
            'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
         !DO IPROC = 0, NUMPES-1
         !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
         !     'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
         !ENDDO

      ENDIF

      dtpic_max = min(dtpic_cfl, dtpic_taup)

      RETURN
      IF(dtsolid.GT.dtpic_max) THEN
         !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
         IF(mype.eq.pe_io) WRITE(*, 2004) dtsolid
      ELSEIF(dtsolid.LT.dtpic_max) THEN

         !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
         IF(mype.eq.pe_io) WRITE(*, 2002) dtpic_max
      ELSE
        !WRITE(*,'(A40,2x,g17.8)') 'DT
        !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
      END IF

      WRITE(unit_log, '(10x,A,2x,3(g17.8))')&
         'DTPIC MINS IN EACH DIRECTION = ', dtpic_min_x, &
         dtpic_min_y, dtpic_min_z
      WRITE(*, '(10x,A,2x,3(g17.8))')'DTPIC MINS IN EACH DIRECTION = ',&
        dtpic_min_x, dtpic_min_y, dtpic_min_z

 1001 FORMAT(3x,'<---------- END INTEGRATE_TIME_PIC ----------')

2001  FORMAT(/1x,70('*'),//,10x,  &
           & 'MOVEMENT UNDESIRED IN INTEGRATE_TIME_PIC: PARTICLE', i5, /,10x, &
           & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10x, &
           & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10x,  &
           & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
           & /1x,70('*'), /10x, &
           & 'DES_VEL_NEW = ',  3(2x, g17.8))

 2004 FORMAT(/10x, &
      & 'DTSOLID SHUD BE REDUCED TO', g17.8)

 2002 FORMAT(/10x, &
      & 'DTSOLID CAN BE INCREASED TO', g17.8)

      RETURN
      END SUBROUTINE integrate_time_pic_garg


!------------------------------------------------------------------------
! subroutine       : des_dbgpic
! Author           : Pradeep G.
! Purpose          : For debugging the pic values
! Parameters       : pstart - start indices of the particle
!                    pend - end indices of the particle
!                    pfreq - optional frequency (when the local count matches the
!                    frequency the filw will be written)
!                    if not send then it prints the file
!------------------------------------------------------------------------

      subroutine des_dbgpic (pstart,pend,pfreq)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      use discretelement
      USE fldvar
      use physprop, only: mmax
      use functions
      implicit none
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
      INTEGER, INTENT(IN) :: pstart,pend
      INTEGER, INTENT(IN), OPTIONAL :: pfreq
!-----------------------------------------------
! Local variables
!-----------------------------------------------
      integer lp,lijk
      integer, save :: lfcount = 0 ,lfreq =0
      character(255) :: filename
!-----------------------------------------------
      if (present(pfreq)) then
         lfreq = lfreq+1
         if (lfreq .ne. pfreq) return
         lfreq =0
      end if
      lfcount = lfcount + 1
      write(filename,'("debug",I3.3)') lfcount
      open (unit=100,file=filename)
      do lp = pstart,pend
         if (is_normal(lp) .or. is_entering(lp) .or. is_exiting(lp)) then
            lijk = pijk(lp,4)
            write(100,*)"positon =",lijk,pijk(lp,1),pijk(lp,2), &
               pijk(lp,3),ep_g(lijk),u_s(lijk,mmax+1)
            write(100,*)"forces =", fc(lp,2),tow(lp,1)
         end if
      end do
      close (100)

      RETURN
      END SUBROUTINE des_dbgpic


!------------------------------------------------------------------------
! subroutine       : des_dbgtecplot
! Author           : Pradeep G.
! Purpose          : prints the tecplot file for particle location
! Parameters       : pstart - start indices of the particle
!                    pend - end indices of the particle
!                    pfreq - optional frequency (when the local count matches the
!                    frequency the filw will be written)
!                    if not send then it prints the file
!------------------------------------------------------------------------

      subroutine des_dbgtecplot (pstart,pend,pfreq)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      use discretelement
      USE fldvar
      use functions
      implicit none
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
      INTEGER, INTENT(IN) :: pstart,pend
      INTEGER, INTENT(IN), OPTIONAL :: pfreq
!-----------------------------------------------
! Local variables
!-----------------------------------------------
      integer lp,lijk
      integer, save :: lfcount = 0 ,lfreq =0
      character(255) :: filename
!-----------------------------------------------

      if (present(pfreq)) then
         lfreq = lfreq+1
         if (lfreq .ne. pfreq) return
         lfreq =0
      end if
      lfcount = lfcount + 1
      write(filename,'("new_tec",I3.3,".dat")') lfcount
      open (unit=100,file=filename)
      write(100,'(9(A,3X),A)') &
         'VARIABLES = ', '"ijk"', '"x"', '"y"', '"vx"', '"vy"', &
         '"ep_g"', '"FCX"' ,'"FCY"', '"TOW"'
      write(100,'(A,F14.7,A)') 'zone T = "' , s_time , '"'
      do lp = pstart,pend
         if (.not.is_nonexistent(lp)) then
            lijk = pijk(lp,4)
            write(100,*)lijk,des_pos_new(lp,1),des_pos_new(lp,2), &
               des_vel_new(lp,1),des_vel_new(lp,2),ep_g(lijk),&
               fc(lp,1),fc(lp,2),tow(1,lp)
         endif
      enddo
      close (100)
      RETURN
      END SUBROUTINE des_dbgtecplot