ccccccccccccccccccccccccc     Program 6.4     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM N_T_FIELD
C
C This program solves the time dependent temperature field
C around a nuclear waste rod in a two-dimensional model.
C
      PARAMETER (M1=1001,M=2001,N=19)
      DIMENSION A(N),B(N),C(N),S(M,N),X(M,N+1),Y(N),
     *          G(N),W(N),V(N),U(N)
C
      HT  = 1.0/(M1-1)
      TC  = 1.0
      XK  = 3153.6
      RA  = 25.0
      RB  = 100.0
      H   = RB/(N+1)
      H2  = H*H
      T0  = 10.0
      S0  = HT*XK*T0/RA**2
      CS  = HT*XK/H2
C
      DO    90  I = 1, M
        DO  80  J = 1, N+1
          X(I,J) = 0.0
   80   CONTINUE
   90 CONTINUE

      DO   100  I = 1, N
        R = H*I
        A(I) =  2.0*(1.0+CS)*R
        B(I) = -(1.0+0.5/I)*CS*R
        C(I) = -(1.0-0.5/I)*CS*R
        IF (R.LT.RA) THEN
          S(1,I) = S0*R
        ELSE
          S(1,I) = 0.0
        END IF
  100 CONTINUE
C
C Assign the source of the radiation heat
C
      DO   200  I = 2, M
        T   = HT*(I-1)
        DO   150  J = 1, N
          R = H*J
          IF (R.LT.RA) THEN
            S(I,J) = R*S0/EXP(T/TC)
          ELSE
            S(I,J) = 0.0
          END IF
  150   CONTINUE
  200 CONTINUE
C
      DO    1000  I = 2, M
C
C Find the values from the last time step
C
        A0 = 2.0*(1.0-CS)*H
        B0 = (1.0+0.5)*CS*H
        C0 = (1.0-0.5)*CS*H
        G(1) = A0*X(I-1,1)
     *        +C0*(4.0*X(I-1,1)-X(I-1,2))/3.0
     *        +B0*X(I-1,2)+S(I-1,1)+S(I,1)
        DO   300  J = 2, N
          R  = J*H
          A0 = 2.0*(1.0-CS)*R
          B0 = (1.0+0.5/J)*CS*R
          C0 = (1.0-0.5/J)*CS*R
          G(J) = A0*X(I-1,J)+B0*X(I-1,J+1)
     *          +C0*X(I-1,J-1)+(S(I-1,J)+S(I,J))
  300   CONTINUE
C
C Find the elements in L and U
C
        W(1) = A(1)
        V(1) = C(1)
        U(1) = B(1)/W(1)
        DO    400 J = 2, N
          V(J) = C(J)
          W(J) = A(J)-V(J-1)*U(J-1)
          U(J) = B(J)/W(J)
  400   CONTINUE
C
C Find the solution of the temperature
C
        Y(1) = G(1)/W(1)
        DO   500 J = 2, N
          Y(J) = (G(J)-V(J-1)*Y(J-1))/W(J)
  500   CONTINUE
C
        X(I,N) = Y(N)
        DO   600 J = N-1,1,-1
          X(I,J) = Y(J)-U(J)*X(I,J+1)
  600   CONTINUE
 1000 CONTINUE
      I = M
      DO    800 J = 1, N+1
        R = J*H
        WRITE (6,9999) R,X(I,J)
  800 CONTINUE
      STOP
 9999 FORMAT (2F16.8)
      END