C    Numerical calculation of Cahn-Hilliard equation
C
C       2-dimensional(!-Dimentional for nx=1)
C
C                         Y.Saito, Waseda University
C                          
      parameter(nt=2,ix=201,iy=201)
      common/ran0/idummy
      dimension c(nt,ix,iy),g(ix,iy),cout(0:100,ix), itim(10)
      dimension st(10,10000),sti(0:100,0:100),str(0:100,0:100)
      dimension ifq(2500)
      data itim/900,1800,3600,7200,18000,36000,72000, 
     * 180000,360000,720000/
      open(unit=4,file='outch.dat')
      idummy=139475347
      nx=1
      ny=101
      a0=2.57e-8
      nsys=nx*ny
      asys=nsys
      pi=3.14159
        fac=2.0*pi/50
           dt=10.
        dx=2.5e-8
  	  r=8.3144
C	  tc=1103
 1001    write(6,*) "Temperature, K(750<T< 850)?"
	  read(5,*) t
          if(t.gt.850) go to 1001
	  write(6,*) "Cr, wt.frac(c<0.49)?"
	  read(5,*) cwt
          write(6,*) "Holding time(=<720000)"
          read(5,*) nzeit
      ak=0.16667*a0**2*(25104.-11.72*t)*1.1
	  afe=55.84
	  acr=51.996
 	  cav=cwt*afe/(acr*(1-cwt)+afe*cwt)
	  cav0=cav*asys
       dcr=0.19*exp(-246000/(r*t))*0.01
C              w=2.0*tc
       am=dcr/(2.0*(25104.-11.72*t)-4.0*r*t)
      dx1=ak/(dx*dx)
      dx2=am/(dx*dx)
      cav1=0.0	
          write(*,*)
          write(*,*)
          write(*,*) "Initial distribution of solute"
	  do 120 i=1,nx
          do 130 j=1,ny
	  call randu(xran)
	  c(1,i,j)=cav+(xran-0.5)*0.05
	  cav1=cav1+c(1,i,j)
  130 continue
  120 continue
      do 140 i=1,nx
      do 150 j=1,ny
	  c(1,i,j)=c(1,i,j)*cav0/cav1
  150 continue
  140 continue
      write(6,*) (c(1,1,j),j=1,ny)
	  do 160 i=1,ny
	  cout(0,i)=c(1,1,i)
  160 continue
      iout=1
      do 110 izeit=1,nzeit
      do 10 i=1,nx
      do 20 j=1,ny
      i1=i-1
      i2=i+1
      j1=j-1
      j2=j+1
	  if(c(1,i,j).le.0.0) c(1,i,j)=0.00001
 	  if(c(1,i,j).ge.1.0) c(1,i,j)=0.99999
	  c0=c(1,i,j)
	  c1=1.-c0
	  dfc=r*t*alog(c0/c1)-2.*(c0-0.5)*(25104.-11.72*t)
      if(i1.lt.1)  i1=nx
      if(i2.gt.nx) i2=1
      if(j1.lt.1)  j1=ny
      if(j2.gt.ny) j2=1
      g1=c(1,i1,j)+c(1,i2,j)+c(1,i,j1)+c(1,i,j2)-4.0*c(1,i,j)
      g(i,j)=(-dx1*g1+dfc)*c0*c1
   20 continue
   10 continue
      do 30 i=1,nx
      do 40 j=1,ny
      i1=i-1
      i2=i+1
      j1=j-1
      j2=j+1
      if(i1.lt.1)  i1=nx
      if(i2.gt.nx) i2=1
      if(j1.lt.1)  j1=ny
      if(j2.gt.ny) j2=1
      f1=g(i1,j)+g(i2,j)+g(i,j1)+g(i,j2)-4.0*g(i,j)
      c(2,i,j)=c(1,i,j)+dt*dx2*f1
   40 continue
   30 continue
      do 50 i=1,nx
      do 60 j=1,ny
	  c(1,i,j)=c(2,i,j)
   60 continue
   50 continue
      if(izeit.eq.itim(iout)) then
	  write(6,1000) izeit
 1000 format("  Time:", i8)
      write(6,*) (c(1,1,i),i=1,ny)
	  do 70 j=1,ny
	  cout(iout,j)=c(1,1,j)
   70 continue
      do 405 kx=0,(nx-1)/2
      do 406 ky=0,(ny-1)/2
      str(kx,ky)=0.0
      sti(kx,ky)=0.0
  406 continue
  405 continue
      do 407 kx=1,625
      st(iout,kx)=0.0
      ifq(kx)=0
  407 continue
      do 310 kx=0,(nx-1)/2
      do 315 ky=0,(ny-1)/2
      ksq=kx*kx+ky*ky
      if(ksq.eq.0) go to 315
      if(ksq.gt.625) go to 315
      do 320 i=1,nx
      do 325 j=1,ny
      xx=kx*i+ky*j
      akr=xx*fac
       str(kx,ky)=str(kx,ky)+cos(akr)*c(1,i,j)
       sti(kx,ky)=sti(kx,ky)+sin(akr)*c(1,i,j)
  325 continue
  320 continue
  315 continue
  310 continue
C
      do 425 kx=0,(nx-1)/2
      do 435 ky=0,(ny-1)/2
      ksq=kx*kx+ky*ky
      if(ksq.eq.0) go to 435
      if(ksq.gt.625) go to 435
      ifq(ksq)=ifq(ksq)+1
       st(iout,ksq)=st(iout,ksq)+(str(kx,ky)**2+sti(kx,ky)**2)/asys
  435 continue
  425 continue
      	  iout=iout+1
	  endif 
  110 continue
       delta=dx*1.0e7

      write(4,*)"d 0 .25hr .5hr 1hr 2hr 5hr 10hr 20hr 50hr 100hr 200hr"

      do 210 i=1,ny

	  disp=delta*i-delta

	  write(4,2000) disp,(cout(j,i),j=0,iout-1)

 2000 format(f10.3,11f8.4)

  210 continue

      write(3,*) "k 0 .25hr .5hr 1hr 2hr 5hr 10hr 20hr 50hr 100hr 200hr"

      do 220 kx=1,625

      if(ifq(kx).ge.1) then

       fsq=ifq(kx)
       st(iout,kx)=st(iout,kx)/fsq
        ak=kx
        akk=fac*sqrt(ak)
	  write(3,3000) akk,(st(j,kx),j=1,iout-1)
 3000 format(f8.2,10f9.2)
      endif
  220 continue
      stop
      end
	  subroutine randu(xran)
	  common/ran0/idummy
	  anum=0.4656612e-9
	  iy=idummy
	  iy=iy*16807
	  if(iy) 15,16,16
   15 iy=iy+2147483647+1
   16 y=iy
      xran=y*anum
	  idummy=iy
	  return
	  end