c***********************************************************************
      SUBROUTINE timestep(t,tswap)
c***  
      INCLUDE 'undefined.h'
      INCLUDE 'dimension.h'
      INCLUDE 'physics.h'
      INCLUDE 'flags.h'
      INCLUDE 'global.h'
      INCLUDE 'pointers.h'
      INTEGER t, tswap, d, ncut, GP_IFETCH
      DOUBLE PRECISION velmax,SUMFIELD,FINDNINF,CFL
      INTRINSIC MOD
      EXTERNAL SUMFIELD,FINDNINF,CFL
      double precision cleandebris, meanmass, remMass

      ncut = GP_ifetch("zapout",0,"DEFAULT")
c     SUBROUTINES: vof, bc_scalar, stresses, momentum, fill0 
c                  (phydebug, printps)

C
C.... Real physical boundaries.
C
      CALL bc_vector(u,v,w,nx,ny,nz)
      CALL bc_c(zcc(ia(ng)),nx,ny,nz)

c Output some fields
c
      if(mod(t,tstore_c).eq.0) then
        call datbfile(u,v,w,zcc(ia(ng)),nx,ny,nz,t)
c        call v5dconv(u,v,w,zcc(ia(ng)),nx,ny,nz,t)
      endif
c***  
c     Checkpoint every tcheck steps
c***
      if (MOD(t, tcheck) .eq. 0) then
         velmax = FINDNINF(u,v,w,nx,ny,nz)
         write (*,*) 'Surfer: check point at time: ',t
         write (*,*) 'max |v|               = ',velmax
         write (*,*) 'total volume fraction = ',
     %        SUMFIELD(zcc(ia(ng)),nx,ny,nz)
         write (*,*) 'total x momentum      = ',
     %        SUMFIELD(u,nx,ny,nz)
         call phy_debug(u,v,w,zu(ia(ng)),zcc(ia(ng)),
     %        nx,ny,nz,t,'at checkpoint')
      endif
c***  
c     alternate first direction of interface motion
c***  
      tswap = tswap + 1
      if (tswap .gt. 3) tswap = 1
      if (MOD(tswap,3) .eq. 0) then
         call swp(w,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),3)
         call swp(u,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),1)
         call swp(v,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),2)
      elseif (MOD(tswap,2) .eq. 0) then
         call swp(v,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),2)
         call swp(w,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),3)
         call swp(u,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),1)
      else 
         call swp(u,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),1)
         call swp(v,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),2)
         call swp(w,zcc(ia(ng)),zei(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),3)
      endif 
      call bc_c(zcc(ia(ng)),nx,ny,nz)

      remMass = cleandebris(zcc(ia(ng)),meanmass)
      call cutout(zcc(ia(ng)),ncut)
      if (MOD(t,tcheck).eq.0) write(*,*) "time = ", t ,
     $     " removed = ",remMass," mean mass = ", mean mass
c***
c     get surface tension contribution to the stress tensor
c***
      if (sigma .eq. 0.d0) then
         call fill0(S11,nx,ny,nz)
         call fill0(S22,nx,ny,nz)
         call fill0(S33,nx,ny,nz)
         call fill0(S12,nx,ny,nz)
         call fill0(S13,nx,ny,nz)
         call fill0(S23,nx,ny,nz)
      else
         call stresses(zcc(ia(ng)),S11,S22,S33,S12,S13,S23,
     %        zd(ia(ng)),zres(ia(ng)),zrhs(ia(ng)),za(ia(ng)),
     %        zb(ia(ng)),zc(ia(ng)),zei(ia(ng)))
      endif
c***
c     N-S + multigrid
c***

      call momentum(u,v,w,zcc(ia(ng)),zu(ia(ng)),S11,S22,S33,S12,
     %        S13,S23,zres(ia(ng)),zrhs(ia(ng)),za(ia(ng)),zb(ia(ng)),
     %        zc(ia(ng)),zei(ia(ng)),zd(ia(ng)),t)
c***
c     check on CFL condition before advancing the interface
c***
      velmax = CFL(u,v,w,nx,ny,nz)
      if (velmax .ge. MAXCFL) then
         write(6,*) "Surfer:warning> t velmax",t, velmax
         stop 'CFL number too large. Try smaller tau'
      endif
c***  
      return
      end