program CN2

c   Kreider 5/27/03  CN on 2d square
c   Block SOR version

      implicit none
      double precision AL(100),AM(100),AR(100),u(100,100),b(100)
      double precision delta(100),x(100),y(100)
      double precision dx,dy,dt,alpha,rx,ry,d,omega
      integer i,j,k,N,M,row,ksor

c     this is just ONE time step
c     assume initial condition is u(x,y,0) = 0 so righthand side is
c     source term only

c     assume homogeneous Dirichlet boundary conditions as well

      N = 61
      M = 41
      dx = 1.d0/dfloat(N-1)
      dy = 1.d0/dfloat(M-1)

      do i=1,N
         x(i) = dfloat(i-1)*dx
      end do
      do j=1,M
         y(j) = dfloat(j-1)*dy
      end do

      dt = 1.d-3
      alpha = 1.d0
      rx = alpha*dt/dx**2
      ry = alpha*dt/dy**2
      d = 1.d0 + rx + ry

c     build main diagonal block
      do i=2,N-1
         AL(i) = -rx/2.d0
         AM(i) = d
         AR(i) = -rx/2.d0
      end do
      AM(1) = 1.d0
      AM(N) = 1.d0
      AL(1) = 0.d0
      AL(N) = 0.d0
      AR(1) = 0.d0
      AR(N) = 0.d0

c     initialization
      do i=1,N
         do j=1,M
            u(i,j) = 0.d0
         end do
      end do

      open(1,file='CN2.input')
      read(1,*) omega
      close(1)

c     SOR loop
      do ksor=1,50
c        block j=1 -- diagonal, zero for this case
         j = 1
         do i=1,N
            delta(i) = 0.d0
         end do
         do i=1,N
            u(i,j) = u(i,j) + delta(i)
         end do
c        blocks j=2,M-1 -- tridiagonal main block
         do j=2,M-1
c           build RHS
            i = 1
            b(i) = -( AM(i)*u(i,j)+AR(i)*u(i+1,j) ) + 0.d0
            b(i) = omega*b(i)
            do i=2,N-1
               b(i) = -(AL(i)*u(i-1,j)+AM(i)*u(i,j)+AR(i)*u(i+1,j))
     .                + 2.d0*exp(x(i)*y(j))
     .                -( -ry/2.d0*u(i,j-1)-ry/2.d0*u(i,j+1) )
            end do
            i = N
            b(i) = -( AL(i)*u(i-1,j)+AM(i)*u(i,j) ) + 0.d0
            b(i) = omega*b(i)   
            call tridiag(AL,AM,AR,b,delta,N)
c           if (j.eq.3) then
c           write(99,*) 'ksor,j',ksor,j
c           write(99,*) (delta(i),i=1,N)
c           write(99,*)
c           end if
            do i=1,N
               u(i,j) = u(i,j) + delta(i)
            end do
         end do
c        block j=M -- diagonal, zero for this case
         j = M
         do i=1,N
            delta(i) = 0.d0
         end do
         do i=1,N
            u(i,j) = u(i,j) + delta(i)
         end do
      end do

      open(1,file='CN2.dat')
      do i=1,N
         do j=1,M
            write(1,*) x(i),y(j),u(i,j)
         end do
         write(1,*)
      end do

      stop
      end

      SUBROUTINE tridiag(a,b,c,r,u,n)
      implicit none
      INTEGER n,NMAX
      double precision a(n),b(n),c(n),r(n),u(n)
      PARAMETER (NMAX=5000)
      INTEGER j
      double precision bet,gam(NMAX)
      if(b(1).eq.0.)pause 'tridag: rewrite equations'
      bet=b(1)
      u(1)=r(1)/bet
      do 11 j=2,n
        gam(j)=c(j-1)/bet
        bet=b(j)-a(j)*gam(j)
        if(bet.eq.0.)pause 'tridag failed'
        u(j)=(r(j)-a(j)*u(j-1))/bet
11    continue
      do 12 j=n-1,1,-1
        u(j)=u(j)-gam(j+1)*u(j+1)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software +%6V+j)D2.