program circle
*     **************************************************************
*     main program to solve the Hele-Shaw boussinesq combustion
*     without the heat loss effect.
*     We use a single layer representation for the pressure. The
*     evolution of the front is done in the theta-L variables.
*
*     input:
*         n:    initial number of points (must be 2**k, for some k);
*         nt:   number of time steps;
*         dt:   time step;
*         sl:   laminar burning speed;
*         pei:  1/Pe where Pe is the Peclet number;
*         nprint: number of printing output.
*
*     **************************************************************
      integer i,j,n,n2,nt,k,nprint,kprint,kunit,ktime
      double precision x(0:8192),y(0:8192),th(0:8192),phi(0:8192),
     +                 curv(0:8192),vf(0:8192),utan(0:8192),ta(0:8192)
      double precision fphio(0:8192),fphi(0:8192),fphin(0:8192)
      double precision pi,h,hp,p2,dcut,t,dt,tol,sl,pei,so,
     +                 xst,yst,xsto,ysto,xstn,ystn,flt,flto,sa,san,
     +                 coef1,v0,a0,b0,ak,bk,omak,opak
      common/darcyc/coef1
      common/salpha/sa
      common/gridsz/h,hp
      common/valupi/pi,p2
*
      data kunit/20/,t/0.D0/

      print *,'enter n,nt,dt,sl,pei,nprint'
      read *,n,nt,dt,sl,pei,nprint
*
*     obtain the constants:
      tol = 1.D-13
      dcut = 10.D0**(-6)
      pi = 4.D0*datan(1.D0)
      p2 = 2.D0*pi
      h = 1.D0/dble(n)
      hp = h*pi
      coef1 = 1.1D-1
      kprint = nt/nprint
*
*     initialize the front and check initial equal-distance condition:
      call initia(n,phi,dcut,tol,x,y)
      so = sa
      call chleng(x,y,n)
      write(kunit,103) (x(j),y(j),j=0,n)
*
*     ***************************************************
*     first step, using the Euler's method:
      
      t = t+dt
      print 102, t,dt,sa
      write(12,101) t,sa

      do 10 j=0,n
         fphio(j) = phi(j)
 10   continue
      call fast(fphio,n)

      call fd2(n,phi,ta,h,dcut)

      do 12 j=0,n
         ta(j) = 1.D0 + ta(j)
         curv(j) = ta(j)/sa
 12   continue

      do 14 j=0,n
         th(j) = phi(j) + p2*dfloat(j)*h
 14   continue

      call potent(x,y,th,curv,vf,n)
      call tangen(n,ta,vf,curv,utan,sl,pei,flto)

      v0 = vf(0) - sl*(1.D0 - pei*curv(0))

      call fast(vf,n)

      xsto = x(0)
      ysto = y(0)
      xst = xsto - v0*dsin(th(0))*dt
      yst = ysto + v0*dcos(th(0))*dt

      san = sa - dt*flto/p2

      n2 = n/2

      a0 = dt*sl*pei/sa/sa
      b0 = dt/sa
      fphi(0) = fphio(0) + b0*utan(0)
      fphi(1) = 0.D0

      do 15 j=1,n2
         k=2*j
         ak = 1.D0 - a0*dfloat(j*j)
         bk = b0*dfloat(j)
         fphi(k) = ak*fphio(k) - bk*vf(k+1) + b0*utan(k)
         fphi(k+1) = ak*fphio(k+1) + bk*vf(k) + b0*utan(k+1)
 15   continue

      do 16 j=0,n
         phi(j) = fphi(j)
 16   continue
      call fsst(phi,n)
      phi(n) = phi(0)

      call recons(n,phi,x,y,xst,yst,dcut)

      sa = san

c     **********************************************************
c     time stepping for the rest:

      do 1000 ktime=2,nt

      t = t+dt
      print 102, t,dt,sa
      write(12,101) t,sa

*     obtain the derivative of phi, ta, through FFT:
      call fd2(n,phi,ta,h,dcut)
      
*     obtain theta_alpha and the curvature:
      do 22 j=0,n
         ta(j) = ta(j) + 1.D0
         curv(j) = ta(j)/sa
 22   continue

*     obtain theta:
      do 23 i=0,n
         th(i) = phi(i) + p2*dble(i)*h
 23   continue

*     calculate the normal velocity (fluid convection + propagation)
*     explicitly, and the tangential motion part theta_alpha U:
      call potent(x,y,th,curv,vf,n)
      call tangen(n,ta,vf,curv,utan,sl,pei,flt)

*     the normal velocity for the reference point, then move the point:
      v0 = vf(0) - sl*(1.D0 - pei*curv(0))

      xstn = xsto - 2.D0*v0*dsin(th(0))*dt
      ystn = ysto + 2.D0*v0*dcos(th(0))*dt

      call fast(vf,n)

      san = sa - 0.5D0*dt*(3.D0*flt-flto)/p2

*     evolve phi in one time step in Fourier space, implicit C-N scheme
*     is used:

      a0 = dt*sl*pei/sa/sa
      b0 = 2.D0*dt/sa
      fphin(0) = fphio(0) + b0*utan(0)
      fphin(1) = 0.D0

      do 25 j=1,n2
         k=2*j
         ak = a0*dfloat(j*j)
         bk = b0*dfloat(j)
         opak = 1.D0 + ak
         omak = 1.D0 - ak
         fphin(k) = (omak*fphio(k) - bk*vf(k+1) + 
     +               b0*utan(k))/opak
         fphin(k+1) = (omak*fphio(k+1) + bk*vf(k) + 
     +               b0*utan(k+1))/opak
 25   continue

c     add points if necessary:

      if( san .gt. 2.D0*so ) then
         print *,'sa and so ',sa,so
         print *,'double the points.'
         do 30 j=0,n
            fphin(j) = 2.D0*fphin(j)
            fphi(j) = 2.D0*fphi(j)
 30      continue
         do 35 j=n+1,2*n
            fphin(j) = 0.D0
            fphi(j) = 0.D0
 35      continue
         n = 2*n
         n2 = 2*n2
         h = 0.5D0*h
         hp = 0.5D0*hp
         so = 2.D0*so
c            dt = 0.5D0*dt
      endif

*     recover phi
      do 40 j=0,n
         phi(j) = fphin(j)
 40   continue
      call fsst(phi,n)
      phi(n) = phi(0)

*     reconstruct the front:
      call recons(n,phi,x,y,xstn,ystn,dcut)

c     renew the information about the front:
      do 45 j=0,n
         fphio(j) = fphi(j)
         fphi(j) = fphin(j)
 45   continue
      sa = san
      flto = flt
      xsto = xst
      ysto = yst
      xst = xstn
      yst = ystn

      if ( mod(ktime,kprint) .eq. 0 ) then
         call chleng(x,y,n)
         kunit = kunit + 1
         print *,'t= ',t,' write to unit ',kunit
         print *,'sa= ',sa
         write(kunit,103) (x(i),y(i),i=0,n)
      endif

 1000 continue

      write(10,101) (i,vf(i),i=1,n)
      write(11,101) (i,curv(i),i=1,n)
        
 101  format('t= ',f12.8,' sa= ',f12.8)
 102  format('t= ',f12.8,' dt= ',f12.8,' sa= ',f12.8)
 103  format(2(f16.12,2x))
      end