c
c
c
c
      PROGRAM Chute
      include 'undefined.h'
      INTEGER tmax,tmovie,tprint,tsave
      INTEGER t,i, ifetch
      DOUBLE PRECISION tphys, sumfield
      INCLUDE 'size.h'
      INCLUDE 'global.f'
      INCLUDE 'pointers.f'
      call initlog()
c      CALL machineid()
      open(unit=11)
      call read_stdin()
      call initialize(tphys,tmax,tprint)
      call closelog()

      i = ifetch("tsave",tsave)
      call flush(6)
      do t=t0,tmax
         if(mod(t,tsave).eq.0) then
            call save(u,v,z(icc(ng)),z(iu(ng)),nx,ny,t)
	    write (*,*) ' --- Writing dump.dat at t = ',t,' ---'
            write(6,*) "zmomentum=",sumfield(v,nx,ny)
	    write (*,*) ' --- Wrote dump.dat at t = ',t,' ---'
            call flush(6)
         endif
         call timestep(tphys,t,tprint,tmovie)
      enddo

      STOP
      end

      SUBROUTINE initialize(tphys,tmax,tprint)
      include 'undefined.h'
      INTEGER tmax,tprint,restart,i, ifetch
      DOUBLE PRECISION tphys, sumfield
      INCLUDE 'size.h'
      INCLUDE 'global.f'
      INCLUDE 'pointers.f'
      call readglob(tmax,tprint)
      restart=0
      i=ifetch("restart",restart)
      write(6,*) "restart=",restart
      call initpvel(nx,ny)
c
      call printhist(tprint)
      call initpressure(nx,ny)
      write(6,*) "pressure initialized"

      if(restart.eq.0) then
	t0=0;
c
c Include here relevant initialisation module
c
c It should be called ini"something". 
c
         tphys=0.0d0
c         call ini1disk(u,v,z(icc(ng)),z(iccol),z(iu(ng)),tau,nx,ny)
c        call ininothing(u,v,z(icc(ng)),z(iu(ng)),tau,nx,ny)
c         call inisquare(u,v,z(icc(ng)),z(iu(ng)),tau,nx,ny)
           write(6,*) "tau=",tau
           write(6,*) "initializing disks"
         call ini1disk(u,v,z(icc(ng)),z(iccol),z(iu(ng)),
     $          xvel,yvel,tau,nx,ny)
c$$$        write(6,*) "x symmetry check for disk initialisation c"
c$$$        call  chkx(cc,nx,ny,0,1)
c$$$        write(6,*) "y symmetry check for disk initialisation c"
c$$$        call  chky(cc,nx,ny,0,1)
         write(6,*) v(12,70)
      else if(restart.eq.1) then
         write(6,*) "reading saved data"
         call iniget(u,v,z(icc(ng)),z(iu(ng)),nx,ny,t0)
         write(6,*) "zmomentum=",sumfield(v,nx,ny)
         write(6,*) "read saved data from t = ",t0
         call datbfile(u,v,z(icc(ng)),z(iu(ng)),nx,ny,nz,t0
     *        ,sigma,tau)
         write(6,*) "wrote saved data from t = ",t0
      else if(restart.eq.2) then
         write(6,*) "reading saved data"
         call iniget(u,v,z(icc(ng)),z(iu(ng)),nx,ny,t0)
         write(6,*) "zmomentum=",sumfield(v,nx,ny)
         write(6,*) "read saved data from t = ",t0
         call datbfile(u,v,z(icc(ng)),z(iu(ng)),nx,ny,nz,t0
     *        ,sigma,tau)
         write(6,*) "wrote saved data from t = ",t0
         stop

c fill pressure with 0 (work around corrupted dump of p filed)

      else if(restart.eq.3) then
         write(6,*) "reading saved data"
         call iniget(u,v,z(icc(ng)),z(iu(ng)),nx,ny,t0)
         call fill0(z(iu(NG)),nx,ny)
         write(6,*) "read saved data from t = ",t0
         write(6,*) " and filled p with 0"
         call datbfile(u,v,z(icc(ng)),z(iu(ng)),nx,ny,nz,t0
     *        ,sigma,tau)
         write(6,*) "wrote saved data from t = ",t0
      endif

      call makebccf(z(icc(ng)),nx,ny)
      call makebcf(u,v,nx,ny)

c
c Used in rhouv version
c
c      call acfromcc(z(icc(ng)),z(ia(ng)),z(ic(ng)),
c     %              z(iei(ng)),z(ires(ng)),nx,ny)
c      call rhouv(rhou,rhov,u,v,z(ia(ng)),z(ic(ng)),rhogas,nx,ny)
c      call makebcf(rhou,rhov,nx,ny)

      return
      end