      PROGRAM flame

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'param2.txt'
      INCLUDE 'common.txt'

      INTEGER ib,ie,jb,je,kb,ke,ini

      REAL time1
      
      LOGICAL mkflam

      REAL upipe(0:pi1,0:pj1,0:pk1)
      REAL vpipe(0:pi1,0:pj1,0:pk1)
      REAL wpipe(0:pi1,0:pj1,0:pk1)
      REAL drpipe
      INTEGER aup(0:i1)

      REAL rhs(imax,jmax,kmax)

      REAL alpha,beta,preexp,incrho,rhoi
      INTEGER iout

      COMMON /indata/ upipe,vpipe,wpipe,drpipe,aup
      COMMON /heatrl/ alpha,beta,preexp,incrho,rhoi

      DATA ib,ie,jb,je,kb,ke /1,imax,1,jmax,1,kmax/

c-------------------------------------------------------------------

c.....read in the stuff from the file "invoer".....

      READ(5,*) istart
      READ(5,*) nstep
      READ(5,*) nsteps
      READ(5,*) ioutp

      READ(5,*) time1
      READ(5,*) time
      READ(5,*) dtmax

      READ(5,*) Re
      READ(5,*) Pr
      READ(5,*) Le
      READ(5,*) Cs
      READ(5,*) Cl

      rdfile='start'
      flname='start'

c.....start main program.....

      CALL mkgrid
      
      CALL pipe(ijet)
      CALL interp

      incrho = 1.0
      rhoi   = 1.0
      alpha  = 0.4
      beta   = 7.0
      preexp = 0! 2000.0 
c      preexp = 0.0 

      IF (istart.eq.3) CALL fkdat
      IF (istart.eq.2) THEN
         CALL loadd(0)
      END IF

      CALL bound(unew,vnew,wnew,tnew,ynew,rhonew,time,ijet,dt,zp,zw,1,ru,rp)

      istep  = nsteps
      dt     = dtmax
      mkflam = .false.

c      CALL dynamic(ib,ie,jb,je,kb,ke)
c      CALL smagor(ib,ie,jb,je,kb,ke)
      CALL viscos(ib,ie,jb,je,kb,ke)

      CALL chkdt

c...................START MAIN LOOP..........................
c............................................................

      ini  = 0
      iout = 0
      DO istep=1+nsteps,nstep+nsteps

         IF(MOD(istep,1) .EQ. 0) CALL chkdt

         time=time+dt

            incrho=incrho+dt
            if(incrho .gt. 1.0) incrho=1.0

c         CALL dynamic(ib,ie,jb,je,kb,ke)
c         CALL smagor(ib,ie,jb,je,kb,ke)
         CALL viscos(ib,ie,jb,je,kb,ke)

c...predictor from sandia scheme...

         CALL advec(ib,ie,jb,je,kb,ke)
         CALL diffu(ib,ie,jb,je,kb,ke)

         CALL predic

         CALL fillps
         CALL psolver(ini,p,ru,rp,dphi,zp,zw)
         ini = 1
         CALL correc

         CALL bound(unew,vnew,wnew,tnew,ynew,rhonew,time,ijet,dt,zp,zw,1,ru,rp)

c...corrector from sandia scheme...

         CALL viscos(ib,ie,jb,je,kb,ke)

         CALL advec2(ib,ie,jb,je,kb,ke)
         CALL diffu2(ib,ie,jb,je,kb,ke)

         CALL predc2

         CALL fillps
         CALL psolver(ini,p,ru,rp,dphi,zp,zw)
         CALL correc

         CALL bound(unew,vnew,wnew,tnew,ynew,rhonew,time,ijet,dt,zp,zw,1,ru,rp)

c...dump fields...

         IF( MOD(istep,1).EQ. 0)THEN

            PRINT*,'dt=',dt,' time=',time
            PRINT*,'alpha=',alpha,'incrho=',incrho,'rhoi=',rhoi

            CALL chkdiv

c            IF(time .GT. 4.) STOP

         END IF

c...filter the scalar fields...

         IF(MOD(istep,1) .EQ. 0 .AND. 1.EQ.0)THEN
            DO k=1,kmax
               DO j=1,jmax
                  DO i=1,imax

                     IF(ynew(i,j,k).gt.1.0) ynew(i,j,k)=1.0
                     IF(ynew(i,j,k).lt.0.0) ynew(i,j,k)=0.0

                     IF(tnew(i,j,k).gt.1.0)THEN
                        tnew(i,j,k)=1.0
                        rhonew(i,j,k)=1.0/((incrho*alpha/(1.0-incrho
     $                       *alpha))*tnew(i,j,k)+rhoi)
                     END IF

                     IF(tnew(i,j,k).lt.0.0)THEN
                        tnew(i,j,k)=0.0
                        rhonew(i,j,k)=1.0/((incrho*alpha/(1.0-incrho
     $                       *alpha))*tnew(i,j,k)+rhoi)
                     END IF

                  END DO
               END DO
            END DO

            CALL scbndy(tnew,ynew,rhonew,ijet)
            
         END IF

c...update fields...

         CALL update

         IF(MOD(istep,100) .EQ. 0)THEN
            PRINT*,'output'

            iout=iout+1

            CALL output(iout)
c            CALL loadd(1)
c            CALL tec3d
         END IF

c.....................DONE ONE TIME ADVANCEMENT.............
c...........................................................

      END DO

      END
