
      SUBROUTINE timestep(tphys,t,tprint,tmovie)
      include 'undefined.h'
      INTEGER tprint,tvelo,tmovie,t,tfront
      INTEGER printdiv,i,j,nfilter,nres
      DOUBLE PRECISION  velmax, sumfield,h,tau1,momc,amplitud,tphys,
     * sumfield2,velmaxliq
      real*8 findninf,findmax,findmin,findninf1
      character*80 c,b
      DOUBLE PRECISION mineps,maxeps,norm_old
      INTEGER cycles,ncycle,npre,npost,printres,relax_times,
     * relax_times_old,ncycle_old
c      
      DOUBLE PRECISION mass_before,mass_after,delta_mass,usang
      DOUBLE PRECISION m_v_x,m_v_y,tauSI,fmaxcfl
      DOUBLE PRECISION tgraph,delta_0 !,m_x_0,m_y_0,y_c_0,
c      DOUBLE PRECISION x_c_before,y_c_before,x_c_after,y_c_after
      DOUBLE PRECISION accellength,acceltime,g0 
      DOUBLE PRECISION gp_DFETCH
      integer gp_IFETCH, ifetch, dfetch
c      
      INCLUDE 'size.h'
      INCLUDE 'global.f'
      INCLUDE 'pointers.f'

      DOUBLE PRECISION mag2(NMX,NMY),mag3(NMX,NMY)

c
      COMMON /cycling/ mineps,maxeps,cycles,ncycle,npre,npost,nres,
     *relax_times,norm_old,printres,relax_times_old,ncycle_old
c      
c
      if(dfetch("maxcfl",fmaxcfl).ne.1) fmaxcfl=MAXCFL
      printdiv=1
      i=ifetch("printdiv",printdiv)
c      
      if(dfetch("tau-si",tauSI).ne.1)   stop 'no tau-si in time step'
c      
c
c  Send important output to a file (unit 66)
c  that will always be opened and closed, so that it can be read while
c  the program runs. 
c
c  First delete the file at t=0

      if(t.eq.t0) then 
         open(unit=66,file='msg',status="unknown")
         close(unit=66,status="delete")

         open(unit=11,file='graph',status="unknown")
         close(unit=11,status="delete") 

         open(unit=55,file='velmaxliq',status="unknown")
         close(unit=55,status="delete") 
      endif

c Read to the end of the file

 100  format(a80)
      open(unit=66,file='msg',status="unknown")
 7771 read(66,100,end=7772,err=7773) c
      go to 7771
 7773 write(6,*) "error reading file msg"
      stop 'msg'
 7772 continue

 900  format(a80)
      open(unit=55,file='velmaxliq',status="unknown")
 971  read(55,900,end=972,err=973) b
      go to 971
 973  write(6,*) "error reading file velmaxliq"
      stop 'velmaxliq'
 972  continue

 700  format(2E18.11,I6,I6) 
      open(unit=11,file='graph',status="unknown")
 771  read(11,700,end=772,err=773)tgraph,delta_0,i,j 
      go to 771
 773  write(6,*) "error reading file graph"
      stop 'graph'
 772  continue



      if(t.eq.t0)  write(6,*) "printdiv=",printdiv
c
c Control CFL condition
c
      velmax = findninf(u,v,nx,ny)

      if( velmax.gt. fmaxcfl) then
         write(6,*) "surfer:warning> t velmax",t, velmax
         call datbfile(u,v,z(icc(ng)),z(iu(ng)),nx,ny,nz,t
     *        ,sigma,tau)
         stop 'CFL number too large. Try smaller tau'
      endif
c
c Some output
c     

c     IRA: 21.10.97 TESTS OF MASS CONSERVATION AND POSITION OF CENTER
c     
      mass_before=sumfield2(z(icc(ng)),nx,ny) 
c      call center_bubble(z(icc(ng)),nx,ny,x_c_before,y_c_before)
      
c     
c     IRA: 21.10.97 TESTS OF MASS CONSERVATION 	
      if(mod(t,tprint).eq.0) then
         write(6,*) ":surfer,CFL> ",t,velmax
         write(66,*)  ":surfer,CFL> ",t,velmax
        velmaxliq = findninf1(u,v,z(icc(ng)),nx,ny)
        write(55,*)  ":surfer,velmaxliq> ",t,velmaxliq
        write(6,*) "Usang =",usang(tau,v,z(icc(ng)),nx,ny)
         write(6,*) "mass=",sumfield2(z(icc(ng)),nx,ny)         
         write(6,*) "zmomentum=",sumfield(v,nx,ny)
         write(6,*) "max C", findmax(z(icc(ng)),nx,ny,i,j), 
     $              "min C", findmin(z(icc(ng)),nx,ny,i,j)
c         call hofx(z(icc(ng)),z(ires(ng)),nx,ny)
c	 amplitud=2.0d0*(findmax2(z(ires(ng)),nx,ny,i)-(ny-2)/2.0d0)
c         amplitud=findmax2(z(ires(ng)),nx,ny,i)-findmin2(z(ires(ng))
c     $    ,nx,ny,j)
c         write(6,*) "amplitude = ",amplitud
         call flush(6)
      endif
c
c Output some fields
c
c Writes a datb file format 
c (http://www.lmm.jussieu.fr/MEMBRES/POPINET/popiplot/index.html)
c
      i=ifetch("tvelo",tvelo)
      if(mod(t,tvelo).eq.0) then
	call datbfile(u,v,z(icc(ng)),z(iu(ng)),nx,ny,nz,t,sigma,tau)
      endif

      tfront = GP_ifetch("tfront",-1,"DEFAULT")
      if(tfront.ne.-1) then
         if(mod(t,tfront).eq.0) then
            call front(z(icc(ng)),nx,ny,tphys)
         endif
      endif
c
c  Define the space step value, deduced from nx as everywhere else
c
      h=1.d0/(nx-2)
c
c  Calls CIAM
c     
      if(mod(t,2).eq.0) then
         call x_swp(u,z(icc(ng)),z(iei(ng)),
     $        z(ires(ng)),z(irhs(ng)),nx,ny)
         call z_swp(v,z(icc(ng)),z(iei(ng)),
     $        z(ires(ng)),z(irhs(ng)),nx,ny)
      else
         call z_swp(v,z(icc(ng)),z(iei(ng)),
     $        z(ires(ng)),z(irhs(ng)),nx,ny)
         call x_swp(u,z(icc(ng)),z(iei(ng)),
     $        z(ires(ng)),z(irhs(ng)),nx,ny)
      end if
c        CIAM_Raus deprecated
c     
c: 21.10.97 TESTS OF MASS CONSERVATION 
	mass_after=sumfield2(z(icc(ng)),nx,ny) 
   	delta_mass=DABS((mass_before-mass_after)/mass_before)
        if(delta_mass.eq.0.0d0) then 
        delta_mass=-16.0d0
        else
        delta_mass=DLOG10(delta_mass)
        endif
        
        if ( ( relax_times_old .ne. relax_times) .or.
     * (cycles .ne. ncycle_old) ) then
        write(11,700)(t-1)*tauSI,delta_mass,ncycle_old,relax_times_old
        write(11,700)t*tauSI,delta_mass,cycles,relax_times
        endif

      if(sigma.eq.0.d0) then
         call fill0(S11,nx,ny)
         call fill0(S22,nx,ny)
         call fill0(S12,nx,ny)
      else
c     
c     gets a parameter to determine whether we filter C before 
c     computing surface tension. Then compute surface tension
c     in two steps. 
c
c     1) The filtered C goes into z(irhs(ng))which  will now be a pointer
c        to the filtered volume fraction field
c     
         nfilter=0
         i=ifetch("nfilter",nfilter) 
         if(nfilter.eq.0) then
            call stresses(S11,S22,S12,z(icc(ng)),z(iei(ng)),
     $                    z(ires(ng)),sigma,tau,nx,ny)    
         else
c
c  z(irhs(ng)) will now be a pointer to the filtered volume fraction field
c
            call filter(z(icc(ng)),z(irhs(ng)),nx,ny)
            call makebccf(z(irhs(ng)),nx,ny)
            do i=1,nfilter-1
               call copy(z(iei(ng)),z(irhs(ng)),nx,ny)
               call filter(z(iei(ng)),z(irhs(ng)),nx,ny)
               call makebccf(z(irhs(ng)),nx,ny)
            enddo
            call stresses(S11,S22,S12,z(irhs(ng)),z(iei(ng)),
     $                    z(ires(ng)),sigma,tau,nx,ny)    
         endif
      endif
c
c     adapt gravity for an acceleration phase for an impacting droplet
c
      if(GP_ifetch("acceleration",0,"DEFAULT").eq.1) then
c  by default make the accel length larger than L_y 
         accellength=GP_dfetch("acceleration_length",0.d0,"SEVERE")
         acceltime= accellength/du
         g0 = (- 2.)*du**2/accellength
c  x = (1/2) g0 t**2 , x = accellength, t=acceltime
         if(tphys.lt.acceltime) then 
            g=g0
         else
            g=0.d0
         endif
      endif

      call momentum(mag2,mag3,u,v,z(icc(ng)),z(iu(ng)),z(ia(ng)),
     *     z(ic(ng)),z(iei(ng)),S11,S22,S12,
     *     z(ires(ng)),z(irhs(ng)),tau,mugas,muliq,
     *     rhogas,g,sigma,nx,ny,printdiv,t,t0) 
c      call adaptative_tau(tau,u,v,z(icc(ng)),z(iu(ng)),z(irhs(ng))
c     * ,nx,ny)
      tphys=tphys+tau
      close(unit=66)
      close(unit=55)
      close(unit=11)
        
      return
      end
