c*********************************************************************
c
c      Program Thermal 2A (implicit differencing on T)
c      Solves the phase-field equations for a pure material
c 
c      FTCS is used on phi
c      Last Modified (11/12/96) - Uses new parameter definitions
c     
c*********************************************************************
c
c
	  implicit double precision (a-h,o-z)

      parameter(nx=1200, ny=600, iskp=100, jskp=50)
c     parameter(nx=80, ny=40, iskp=10, jskp=5)
      parameter(nxs=2, nys=2)
      common/soln/ p(0:nx,0:ny),t(0:nx,0:ny),time
      common/solv/ ax(0:nx),bx(0:nx),cx(0:nx),dx(0:nx),
     *             ay(0:ny),by(0:ny),cy(0:ny),dy(0:ny)
      common/peqn/ pcx,pcy,pcxx,pcyy,pcxy,pepc,pefc
      common/teqn/ rx,ry,sx,sy
      common/prms/ ramp,afac,stn,adel,tol,delt,ndelt
      character*20 fname
      data fname/'tem800_400.bin'/
      data tbulk/-1.0/, tol/1.0e-10/
      tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
c
c read the data...
c
      open(10,file='phase1200x600_input.dat',status='unknown')
      read(10,*)
      read(10,*) nctn,nstep,iprt,iout,ndelt
      read(10,*)
      read(10,*) delt,xl,yl,adel,ramp
      read(10,*)
      read(10,*) eps,parm,alpha,stn,pdel
      read(10,*)
      read(10,*) nics,r0,y0,amp,ts
      close(10)
c
c write out the input data
c
	write(*,*) 'Phase-Field Solution - Code 2 FTCS implicit/explicit'
	write(*,10) eps,parm,alpha,stn,adel,ramp,delt,ts,xl,yl,
     *          r0,y0,amp,pdel,nstep,ndelt,nics,iprt,iout,nctn
 10     format(/1x,
     *  '  eps = ',1pe13.5,5x,' parm = ',1pe13.5/1x,
     *  'alpha = ',1pe13.5,5x,'  stn = ',1pe13.5/1x,
     *  ' adel = ',1pe13.5,5x,' ramp = ',1pe13.5/1x,
     *  ' delt = ',1pe13.5,5x,'   ts = ',1pe13.5/1x,
     *  '   xl = ',1pe13.5,5x,'   yl = ',1pe13.5/1x,
     *  '   r0 = ',1pe13.5,5x,'   y0 = ',1pe13.5/1x,
     *  '  amp = ',1pe13.5,5x,' pdel = ',1pe13.5/1x,
     *  'nstep = ',i6,7x,'ndelt = ',i6,7x,' nics = ',i6/1x,
     *  ' iprt = ',i6,7x,' iout = ',i6,7x,' nctn = ',i2/)
c
c dump the data .. first open dump channel (10)
c
c        open(11,file=fname,form='unformatted',
c     *  status='unknown')
c 	write(11) eps,parm,alpha,stn,adel,ramp,delt,ts,
c     *            xl,yl,r0,y0,amp,pdel,nx,ny,nstep,iout
c
c define required constants
c
      pi = 4.0*atan(1.0)
c
      eps2=eps**2
      pmob = parm/eps2
      afac = eps*alpha*stn
c
      delx=xl/float(nx)
      dely=yl/float(ny)
c
      delt1 = delt/real(ndelt)
c**** modify stefan number for different time steps:
      stn = real(ndelt)*stn
c     
      pcx  =  0.5/delx
      pcy  =  0.5/dely
      pcxx =  1.0/delx**2
      pcyy =  1.0/dely**2
      pcxy =  pcx*pcy
      pepc =  delt1*pmob*eps2
      pefc =  delt1*pmob
c
c     define ADI constants and arrays
c
      rx =  delt/(delx**2)
      ry =  delt/(dely**2)
      sx = 2.0*(1.0 - rx)
      sy = 2.0*(1.0 - ry)
c
      do 20 i=0,nx 
        ax(i) = -rx 
        bx(i) = 2.*(1.+rx) 
        cx(i) = -rx
 20   continue
c
      cx(0)  = cx(0) + ax(0) 
      ax(nx) = ax(nx) + cx(nx)
c      
      do 22 i=1,nx
        dx(i) = ax(i)/bx(i-1)
        bx(i) = bx(i) - dx(i)*cx(i-1)
 22   continue
c
      do 25 j=0,ny 
        ay(j) = -ry 
        by(j) = 2.*(1+ry) 
        cy(j) = -ry
 25   continue
c
      cy(0)  = cy(0) + ay(0)
      ay(ny) = ay(ny) + cy(ny)
c      
      do 28 j=1,ny
        dy(j) = ay(j)/by(j-1)
        by(j) = by(j) - dy(j)*cy(j-1)
 28   continue
c
c initialize the fields
c
      if(nctn.eq.0) then
      time = 0.0
c
      nyic = ny
c     if(nics.eq.3) nyic = ny/2
      do 40 j = 0,nyic
      y = j*dely
c
      if(nics.eq.1) then
        do 30 i = 0,nx
        x = i*delx
          if(x.eq.0.0) then
          theta = pi/2.0
          else
          theta = atan(amp*y/x)
          endif
        r = sqrt((amp*x)**2 + (y)**2)
        rp = r0*(1.0 + pdel*cos(y0*theta))
        p(i,j) = 0.5*(1.0 + tanh((r-rp)/(2.*eps)))
        t(i,j) = -p(i,j)
c     if(p(i,j).lt.0.95) then
c       t(i,j) = 0.0
c       rt = r
c     else
c       t(i,j) = -(1.0 - exp(-y0*(r-rt)))
c     end if
 30     continue
c
      else
c
        do 35 i = 0,nx
	  x = i*delx
          if(x.eq.0.0) then
          theta = pi/2.0
          else
          theta = atan(amp*y/x)
          endif
          r = sqrt((amp*x)**2 + (y)**2)
          rp = r0*(1.0 + pdel*cos(y0*theta))
        if(r.le.rp) then
   	  p(i,j) = 0.0
	  t(i,j) = ts
        else
   	  p(i,j) = 1.0
	  t(i,j) = tbulk	
        end if
 35     continue
      endif
      if(nics.eq.3) then
      jr = j
      nx2 = nx/2   
      do 38 i=0,nx2
      ir = nx-i 
      p(ir,jr) = p(i,j)
      t(ir,jr) = t(i,j)
 38   continue
      endif
 40   continue
      else
      open(12,file='temi2a.in',form='unformatted',
     $  status='unknown')
        read(12) time
        read(12) ((t(ii,jj), ii = 0,nx),jj=0,ny)
        read(12) ((p(ii,jj), ii = 0,nx),jj=0,ny)
      endif

c
c output the initial conditions
c
c     call print(iskp,jskp)
c
c dump initial profile
c
c        write(11) time
c        write(11) ((t(ii,jj), ii = 0,nx,nxs),jj=0,ny,nys)
c        write(11) ((p(ii,jj), ii = 0,nx,nxs),jj=0,ny,nys)
c
c now time step
c
      do 100 indx = 1,nstep

      time = time + delt

      call update
c
c output the results
c
c      if(mod(indx,iprt).eq.0) call print(iskp,jskp)
c
c output tip data
c
      call tip(p,t,plast,nx,ny,indx,iprt,time,delx,dely)
c
c dump the results as well on channel 11
c
c      if(mod(indx,iout).eq.0) then
c        write(11) time
c        write(11) ((t(ii,jj), ii = 0,nx,nxs),jj=0,ny,nys)
c        write(11) ((p(ii,jj), ii = 0,nx,nxs),jj=0,ny,nys)
c      end if
100   continue
c
c dump solution for restart
c 
c       open(14,file='temi2a.out',form='unformatted',
c     $ status='unknown')
c       write(14) time
c       write(14) ((t(ii,jj), ii = 0,nx),jj=0,ny)
c       write(14) ((p(ii,jj), ii = 0,nx),jj=0,ny)
c
      stop
      end
c***************************************************************************
c
c this subroutine advances the solution one time step using simple
c      explicit euler
c
c***************************************************************************
      subroutine update
       parameter(nx=1200, ny=600, nx1=nx-1, ny1=ny-1, 
     *          nxp1 = nx+1, nyp1 = ny+1)

c      parameter(nx=80, ny=40, nx1=nx-1, ny1=ny-1, nxp1=nx+1, 
c    *   nyp1=ny+1)

c
	  implicit double precision (a-h,o-z)
      real ps(0:nx,0:ny),ts(0:nx,0:ny)
      real uu(0:nx,0:ny),ff(0:nx,0:ny)
c
      common/soln/ p(0:nx,0:ny),t(0:nx,0:ny),time
      common/solv/ ax(0:nx),bx(0:nx),cx(0:nx),dx(0:nx),
     *             ay(0:ny),by(0:ny),cy(0:ny),dy(0:ny)
      common/peqn/ pcx,pcy,pcxx,pcyy,pcxy,pepc,pefc
      common/teqn/ rx,ry,sx,sy
      common/prms/ ramp,afac,stn,adel,tol,delt,ndelt
c
c     phase field:
c
c**** calculate the nonlinear term:
      do 25 j=0,ny
      do 25 i=0,nx
      g = p(i,j)*(1.0 - p(i,j))
c     xn = ramp*2.0*(ranf() - 0.5) - 0.5
      xn = -0.5
      ff(i,j) = g*(p(i,j) + xn + 30.*afac*g*t(i,j))
 25   continue
c
c interior points
c
      do 125 j = 1,ny1
      do 125 i = 1,nx1

      dpx = pcx*(p(i+1,j) - p(i-1,j))
      dpy = pcy*(p(i,j+1) - p(i,j-1))
      gradm2 = dpx*dpx + dpy*dpy

      dpxx = pcxx*(p(i+1,j) - 2.0*p(i,j) + p(i-1,j))
      dpyy = pcyy*(p(i,j+1) - 2.0*p(i,j) + p(i,j-1))
      del2p = dpxx + dpyy

      ps(i,j) = p(i,j) + pepc*del2p + pefc*ff(i,j)

      if(gradm2.gt.tol) then
      gradm = sqrt(gradm2)
      dpxy  = pcxy*(p(i+1,j+1) - p(i+1,j-1) - p(i-1,j+1) + p(i-1,j-1))
      dthex = (dpx*dpxy - dpy*dpxx)/gradm2
      dthey = (dpx*dpyy - dpy*dpxy)/gradm2
      sthe = dpy/gradm
      cthe = dpx/gradm
      cthe2 = cthe*cthe
      cnthe = 8.0*cthe2*(cthe2 - 1.0) + 1.0
      snthe = 4.0*cthe*sthe*(2.0*cthe2 - 1.0)

      ps(i,j) = ps(i,j)
     *        + pepc*adel*(cnthe*(2.0 + adel*cnthe)*del2p
     *        -  8.0*snthe*(1.0 + adel*cnthe)*(dthex*dpx + dthey*dpy)
     *        - 16.0*(cnthe + adel*(cnthe**2 - snthe**2))
     *        *(dthey*dpx - dthex*dpy))
      endif
 125    continue
c
c edges
c
	do 150 j = 0,ny
        jm1 = j-1 
        jp1 = j+1 
        if(j.eq.0)  jm1 = 1
        if(j.eq.ny) jp1 = ny1
c
c left edge
c
        dpy  = pcy*(p(0,jp1) - p(0,jm1))
        dpxx = pcxx*2.0*(p(1,j) - p(0,j))
        dpyy = pcyy*(p(0,jp1) - 2.0*p(0,j) + p(0,jm1))
        del2p = dpxx + dpyy
        gradm2 = dpy*dpy
 
      ps(0,j) = p(0,j) + pepc*del2p + pefc*ff(0,j)

      if(gradm2.gt.tol) then
      dthex =  -dpy*dpxx/gradm2
c
      ps(0,j) = ps(0,j) + pepc*adel*((2.0 + adel)*del2p 
     *          + 16.0*(1.0 + adel)*dthex*dpy)
      endif
c
c right edge
c
        dpy  = pcy*(p(nx,jp1) - p(nx,jm1))
        dpxx = pcxx*2.0*(p(nx1,j) - p(nx,j))
        dpyy = pcyy*(p(nx,jp1) - 2.0*p(nx,j) + p(nx,jm1))
        del2p = dpxx + dpyy
        gradm2 = dpy*dpy
 
      ps(nx,j) = p(nx,j) + pepc*del2p + pefc*ff(nx,j)

      if(gradm2.gt.tol) then
      dthex =  -dpy*dpxx/gradm2
c
      ps(nx,j) = ps(nx,j) + pepc*adel*((2.0 + adel)*del2p 
     *           + 16.0*(1.0 + adel)*dthex*dpy)
      endif


 150     continue

	do 175 i = 0,nx
        im1 = i-1
        ip1 = i+1
        if(i.eq.0)  im1 = 1
        if(i.eq.nx) ip1 = nx1
      dpx  = pcx*(p(ip1,0) - p(im1,0))
      dpxx = pcxx*(p(ip1,0) - 2.0*p(i,0) + p(im1,0))
      dpyy = pcyy*2.0*(p(i,1) - p(i,0))
      del2p = dpxx + dpyy
      gradm2 = dpx*dpx
c
c bottom edge
c
      ps(i,0) = p(i,0) + pepc*del2p + pefc*ff(i,0)

      if(gradm2.gt.tol) then
      dthey =  dpx*dpyy/gradm2
c
      ps(i,0) = ps(i,0) + pepc*adel*((2.0 + adel)*del2p 
     *          - 16.0*(1.0 + adel)*dthey*dpx)
      endif
c
c top edge
c
      dpx  = pcx*(p(ip1,ny) - p(im1,ny))
      dpxx = pcxx*(p(ip1,ny) - 2.0*p(i,ny) + p(im1,ny))
      dpyy = pcyy*2.0*(p(i,ny1) - p(i,ny))
      del2p = dpxx + dpyy
      gradm2 = dpx*dpx

      ps(i,ny) = p(i,ny) + pepc*del2p + pefc*ff(i,ny)

      if(gradm2.gt.tol) then
      dthey =  dpx*dpyy/gradm2
c
      ps(i,ny) = ps(i,ny) + pepc*adel*((2.0 + adel)*del2p 
     *           - 16.0*(1.0 + adel)*dthey*dpx)
      endif

 175     continue
c
        do 180 j = 0,ny
        do 180 i = 0,nx
c        g = p(i,j)*(1.0 - p(i,j))
c        uu(i,j) = 30.*stn*g*g*(ps(i,j) - p(i,j))
        go = 1.0 - p(i,j)
        gn = 1.0 - ps(i,j)
        uold = -go*go*go*(6.0*go*go - 15.0*go + 10.0)
        unew = -gn*gn*gn*(6.0*gn*gn - 15.0*gn + 10.0)
        uu(i,j) = -(unew - uold)/stn
 180    continue 
c   
       do 200 j = 0,ny
       do 200 i = 0,nx
       p(i,j) = ps(i,j)
200    continue
c
c     temp field:
c
c**** gradient bc's in both directions. 
c     x sweep 
      do 220 j=0,ny 
      jm1 = j-1 
      jp1 = j+1 
      if(j.eq.0)  jm1 = 1
      if(j.eq.ny) jp1 = ny1
      do 220 i=0,nx 
      ff(i,j) = uu(i,j) + ry*t(i,jp1) + sy*t(i,j) + ry*t(i,jm1) 
220   continue
      do 230 i=1,nx 
      do 230 j=0,ny
      ff(i,j) = ff(i,j) - dx(i)*ff(i-1,j) 
230   continue
      do 240 j=0,ny
      ts(nx,j) = ff(nx,j)/bx(nx)
240   continue
      do 250 i=nx1,0,-1
      do 250 j=0,ny
      ts(i,j) = (ff(i,j) - cx(i)*ts(i+1,j))/bx(i)
250   continue
c
c**** y sweep 
c 
      do 270 i=0,nx 
      im1 = i-1 
      ip1 = i+1 
      if(i.eq.0)  im1 = 1
      if(i.eq.nx) ip1 = nx1
      do 270 j=0,ny 
      ff(i,j) = uu(i,j) + rx*ts(ip1,j) + sx*ts(i,j) + rx*ts(im1,j)
270   continue
      do 280 j=1,ny
      do 280 i=0,nx 
      ff(i,j) = ff(i,j) - dy(j)*ff(i,j-1) 
280   continue
      do 290 i=0,nx
      t(i,ny) = ff(i,ny)/by(ny)
290   continue
      do 300 j=ny1,0,-1
      do 300 i=0,nx
      t(i,j) = (ff(i,j) - cy(j)*t(i,j+1))/by(j)
300   continue
c
      return
      end
c***************************************************************************
c
c     subroutine to print the data at specific time levels:
c
c***************************************************************************
      subroutine print(iskp,jskp) 
       parameter(nx=1200, ny=600)
c      parameter(nx=80, ny=40)
c
	  implicit double precision (a-h,o-z)
      common/soln/ p(0:nx,0:ny),t(0:nx,0:ny),time
      data iw/6/
c
      write(iw,101) 
      write(iw,103) time,(i,i=0,nx,iskp)
      do 25 j=0,ny,jskp
        write(iw,105) j,(p(i,j),i=0,nx,iskp)
 25   continue
c
      write(iw,107) 
      write(iw,103) time,(i,i=0,nx,iskp)
      do 50 j=0,ny,jskp
        write(iw,105) j,(t(i,j),i=0,nx,iskp)
 50   continue
c
101   format(/1x,'phase-field values:'/)
103   format(1x,'time = ',1pe11.4/,1x,'j',11(8x,i3)) 
105   format(1x,i3,2x,11(1x,e10.3))
107   format(/1x,'temperature values:'/)
      return
      end
c****************************************************************************
c
c subroutine to determine tip velocity and radius of curvature
c
c****************************************************************************
      subroutine tip(p,t,plast,nx,ny,indx,iprt,time,delx,dely)
	  implicit double precision (a-h,o-z)
      dimension p(0:nx,0:ny),t(0:nx,0:ny)
      data ientry/0/
      save ientry

      ientry = ientry + 1
      if(ientry.eq.1)  open(50,file='tip.dat',status='unknown',
     $                      form='unformatted')
c
c first search for interval containing the interface at phi = 0.5
c
      i = -1
 10   i = i + 1
       test = (p(i,0)-0.5)*(p(i+1,0)-0.5)
        if(test.le.0.0) then
          ileft = i
	  iright = ileft + 1
        goto 20
        end if

      if(i.le.nx-1) goto 10
      if(i.eq.nx) return

 20   continue
c
c find difference from each end
c
	dphil = p(ileft,0) - 0.5
	dphir = 0.5 - p(iright,0)
        dphi = p(ileft,0) - p(iright,0)
c
c compute position of phi = 1/2
c
        tpos = delx*ileft + dphil*delx/dphi 
c
c compute tip temp
c
        tt = t(ileft,0) + dphil*(t(iright,0)-t(ileft,0))/dphi 
c        write(*,*) ileft,iright,dphil,dphi,dphir
c
c compute interface tip curvature at left and right points
c
        phixl = (p(ileft+1,0) - p(ileft-1,0))/(2.0*delx)
        phixr = (p(iright+1,0) - p(iright-1,0))/(2.0*delx)

       phiyyl = 2.0*(p(ileft,1) - p(ileft,0))/(dely**2)
       cl  = phiyyl/phixl

       phiyyr = 2.0*(p(iright,1) - p(iright,0))/(dely**2)
       cr  = phiyyr/phixr
        
c
c compute interpolant of curvature
c
       c = cl + dphil*(cr-cl)/dphi
c
c
        if(mod(indx,iprt).eq.0) then

        write(*,'(//1x,'' Time = '',1pe13.5)') time
        write(*,'(/1x,'' Tip Position = '',1pe13.5/1x,
     $  '' Tip Temperature = ''1pe13.5/1x,
     $  '' Tip Radius = '',3(1pe13.5,1x))//') 
     $  tpos,tt,1.0/cl,1.0/c,1.0/cr

c
c write out phase field and temp trhu tip
c
        write(*,'(10(1x,f8.3))') (i*delx,i=ileft-4,iright+4)  
        write(*,'(10(1x,f8.3))') (p(i,0),i=ileft-4,iright+4)  
        write(*,'(10(1x,f8.3))') (t(i,0),i=ileft-4,iright+4)
c
        end if
        write(50) time,tpos,1/c,tt
c
       return
       end

