program circle
*     **************************************************************
*
*     main program to solve the problem of front propagation in a
*     Hele-Shaw cell with the Boussinesq approximation.
*     Heat loss effect is NOT considered here.
*
*     We use a single layer representation for the pressure. The
*     evolution of the front is described by the theta-L variables.
*
*     The burning speed has a curvature dependence:
*         s = sl*exp(-kappa/Pe)
*     so the cusps are smoothed out by the curvature.
*
*     A potential flow can be added to correct the normal velocity
*     along the boundary. The elliptic problem is solved by a
*     multigrid code.
*
*     *time discretization is first order in this code.*
*
*     input:
*         nbd:  flag to check if we add a potential flow to correct the
*               boundary velocity. (0 for no correction and 1 for
*               correction).
*         n:    initial number of points (must be 2**k);
*         nt:   number of time steps;
*         dt:   time step;
*         sl:   laminar burning speed for planar front;
*         pei:  1/Pe where Pe is the Peclet number;
*         nprint: number of times to output front positions.
*
*     key variables:
*         (x,y): front positions;
*         phi:   theta - 2 pi alpha (the periodic part of theta);
*         fphi:  Fourier transform of phi;
*         vf:    normal velocity of the front due to the flow field;
*         ut:    tangential velocity of the front;
*         tauf:  Fourier transform of (theta_alpha ut).
*
*     **************************************************************
      character*10 frame
      integer i,j,n,n2,nt,k,nprint,kprint,kunit,ktime,nx,ny,mxy,nxy2,
     +        nbd,iframe
      parameter(nx=128,ny=128)
      parameter(nxy2=(nx+2)*(ny+2),mxy=4*nxy2)
      double precision x(0:8192),y(0:8192),th(0:8192),phi(0:8192),
     +                 curv(0:8192),vf(0:8192),ut(0:8192),
     +                 tauf(0:8192),ta(0:8192)
      double precision fphi(0:8192)
      double precision padd(mxy),hsrc(mxy),res(mxy)
      double precision bl(ny),br(ny),bb(nx),bt(nx)
      double precision uadd(0:nx,ny),vadd(nx,0:ny)
      double precision pi,h,hp,p2,dcut,t,dtin,dt,tol,sl,pei,
     +                 xst,yst,flt,sa,so,
     +                 coef1,v0,a0,b0,ak,bk,hm,hinv,
     +                 rad,pt
      parameter(tol=1.D-13,dcut=1.D-10)
      common/darcyc/coef1
      common/salpha/sa
      common/gridsz/h,hp
      common/meshsz/hm,hinv
      common/valupi/pi,p2
      common/inicur/rad,pt
*
      data kunit,iframe/20,0/,t/0.D0/

      print *,'enter nbd (0 no and 1 yes)'
      read *,nbd
      print *,'enter n,nt,dtin,sl,pei,nprint'
      read *, n,nt,dtin,sl,pei,nprint
      print *,'enter initial radius (rad) and perturbation amount (pt):'
      read *, rad,pt
*
*     write input to an output file:
      open(12,file='run.info',status='unknown')
      open(13,file='area_t',status='unknown')
      open(15,file='cengra',status='unknown')
      if(nbd .eq. 0) then
        write(12,*) 'no boundary correction applied.'
      else
        write(12,*) 'boundary correction applied.'
      endif
      write(12,*) 'initial number of points = ',n
      write(12,*) 'number of time steps = ',nt
      write(12,*) 'dt_max = ',dtin
      write(12,*) 'laminar burning speed = ',sl
      write(12,*) '1/Pe = ',pei
      write(12,*) 'number of output frames = ',nprint
      write(12,*) 'initial radius = ',rad
      write(12,*) 'perturbation amount = ',pt
      close(12)
*
*     obtain the constants:
      pi = 4.D0*datan(1.D0)
      p2 = 2.D0*pi
      h = 1.D0/dfloat(n)
      hp = h*pi
      hinv = dfloat(ny)
      hm = 1.D0/hinv
      coef1 = 1.1D-1
      kprint = nt/nprint
*
*     initialize the front and check initial equal-distance condition:
      call initia(n,phi,dcut,tol,x,y)
c      call contin(n,phi,dcut,tol,x,y)
      so = sa
      call chleng(x,y,n)
      write(frame,'(''front.'',i4.4)') iframe
      open(20,file=frame,status='unknown')
      write(20,103) (x(i),y(i),i=0,n)
      close(20)
      xst = x(0)
      yst = y(0)
*     obtain the initial phi^:
      do 10 j=0,n
         fphi(j) = phi(j)
 10   continue
      call fast(fphi,n)

*
*     ***************************************************
*     first-order forward Euler time stepping:
      
      do 1000 ktime=1,nt

*     d (phi) / d alpha:
      call fd2(n,phi,ta,h,dcut)

*     d theta / d alpha and the curvature:
      do 12 j=0,n
         ta(j) = 1.D0 + ta(j)
         curv(j) = ta(j)/sa
 12   continue

*     theta as a function of alpha
      do 14 j=0,n
         th(j) = phi(j) + p2*dfloat(j)*h
 14   continue

*     evaluate the single layer potential:
      call potent(x,y,th,curv,vf,n,dcut)
      if(nbd .eq. 1) then
         call boundf(x,y,th,bl,br,bb,bt,n,nx,ny)
         call setzro(hsrc,mxy)
         call poisson(hsrc,padd,res,bl,br,bb,bt,nx,ny,mxy)
         call evapot(x,y,th,vf,bl,br,bb,bt,padd,uadd,vadd,nx,ny,n)
         call kfilter(n,dcut,vf)
      endif
*     evaluate the tangential motion:
      call tangen(n,ta,vf,curv,ut,tauf,sl,pei,flt,dcut)
*     correct the curvature dependence of the burning velocity:
      call curcor(n,curv,vf,sl,pei,dcut)

      v0 = vf(0) - sl*dexp(-pei*curv(0))
*
*     check the time step:
      call timest(vf,ut,curv,n,dtin,dt)

*     calculate the total area via the div theorem:
      call clarea(t,n,x,y,th)

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

      call fast(vf,n)

      xst = xst - v0*dsin(th(0))*dt
      yst = yst + v0*dcos(th(0))*dt

      sa = sa - dt*flt/p2

      n2 = n/2

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

*     time stepping of phi^:
      do 15 j=1,n2
         k=2*j
         ak = 1.D0 + a0*dfloat(j*j)
         bk = b0*dfloat(j)
         fphi(k) = (fphi(k) - bk*vf(k+1) + b0*tauf(k))/ak
         fphi(k+1) = (fphi(k+1) + bk*vf(k) + b0*tauf(k+1))/ak
 15   continue

c     add points if necessary:

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

*     recover phi from phi^:
      do 16 j=0,n
         phi(j) = fphi(j)
 16   continue
      call fsst(phi,n)
      phi(n) = phi(0)

*     reconstruct the curve from phi:
      call recons(n,phi,x,y,xst,yst,dcut)


      if ( mod(ktime,kprint) .eq. 0 .or.
     +     dt .le. 0.5D0*dtin ) then
         call chleng(x,y,n)
         iframe = iframe + 1
c
c        output front positions:
         write(frame,'(''front.'',i4.4)') iframe
         print 105, frame
         print *,'fphi(n) ',fphi(n)
         open(20,file=frame,status='unknown')
         write(20,103) (x(i),y(i),i=0,n)
         close(20)
c
c        output Fourier transform coefficients:
         write(frame,'(''fcoef.'',i4.4)') iframe
         print 105, frame
         open(20,file=frame,status='unknown')
         write(20,103) (dfloat(i),2.D0*h*
     +                  dsqrt(fphi(2*i)**2+fphi(2*i+1)**2),i=0,n2)
         close(20)
c
c        output curvature:
         write(frame,'(''cvt.'',i4.4)') iframe
         print 105, frame
         open(20,file=frame,status='unknown')
         write(20,103) (dfloat(2*i)*hp,curv(i),i=0,n)
         close(20)

      endif

 1000 continue

 101  format('t= ',f12.8,' sa= ',f12.8)
 102  format('t= ',f12.8,' dt= ',f12.8,' sa= ',f12.8)
 103  format(2(f16.10,2x))
 105  format('write to file ',a10)
      end