c-----------------------------------------------------------------------
c Find steady-state temperature profile across the pipe wall set up by heat
c conduction.
c-----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Declare variables (Divide into 10 subintervals) ----------------------
      parameter (n=10)
      real r(n+1)
      real*8 T(n+1)

c Output a program header ----------------------------------------------
      print *, '-------------------------------------------------------'
      print *, 'Find steady-state temperature profile across the pipe  '
      print *, 'wall set up by heat conduction.                        '
      print *, '-------------------------------------------------------'
      print *, ' '

c Specify the boundary condition ---------------------------------------
      r_in  = 0.05
      r_out = 0.10
      T_in  = 200.
      T_out =  80.

c Call a routine to find the temperature profile -----------------------
      call temp(r_in, r_out, T_in, T_out, n, r, T)

c Output the answer ----------------------------------------------------
      print *, 'The temperature profile is '
      print *, '  -----------------------------'
      print *, '  Point    Radius   Temperature'
      print *, '    #       (m)         (C)    '
      print *, '  -----------------------------'
      print 650, (i, r(i), T(i), i=1, n+1)
      print *, '  -----------------------------'

c Some formats ---------------------------------------------------------
650   format(i6, f11.3, f11.1)
      end


c-----------------------------------------------------------------------
      subroutine temp(r_in, r_out, T_in, T_out, n, r, T)
c-----------------------------------------------------------------------
c The following differential equation describes the steady-state temperature
c profile across the pipe wall set up by heat conduction.
c
c   dT/dr + r (d^2T/dr^2) = 0       BC: T (r_in)  = T_in
c                                       T (r_out) = T_out
c
c The finite difference representation is:
c
c   (1/2+r/dr)*T(i+1) - 2*r/dr*T(i) + (-1/2+r/dr)*T(i-1) = 0   i=2,3,...,10
c
c where r = r_in+dr*(i-1)       i=1, 2,..., n+1
c
c An example with n=10
c i=1:  T(1) = T_in        ... specified
c i=2:  (-0.5+0.055/0.005)*T(1) -2*0.055/0.005*T(2) + (0.5+0.055/0.005)*T(3) = 0
c i=3:  (-0.5+0.060/0.005)*T(2) -2*0.060/0.005*T(3) + (0.5+0.060/0.005)*T(4) = 0
c i=4:  (-0.5+0.065/0.005)*T(3) -2*0.065/0.005*T(4) + (0.5+0.065/0.005)*T(5) = 0
c  :                                                                         :
c i=10: (-0.5+0.095/0.005)*T(9) -2*0.095/0.005*T(10)+ (0.5+0.095/0.005)*T(11)= 0
c i=11: T(11) = T_out       ... specified
c
c Variables used ...
c   r_in  = radius of the inner wall in meter (input)
c   r_out = radius of the outer wall in meter (input)
c   T_in  = temperature at the inner wall in Celsius (input)
c   T_out = temperature at the outer wall in Celsius (input)
c   n     = number of divisions or intervals (input)
c   r(i)  = radius at point i (output)
c   T(i)  = temperature at point i (output)
c-----------------------------------------------------------------------
c Programming Note: call TM of Riggs to solve A*x=b
c   Vectors Al, Ad, & Ar are the coefficients to the left, center, and right of
c   the diagonal, respectively.   Beta and gam are working vectors.
c   All the input and output variables to TM must be double precision!
c MS Linker Bug: "link/e" creates an invalid .exe file because real*4 and
c                real*8 are all mixed up.  Use plain "link" without packing
c                with the "/e" attribute.
c-----------------------------------------------------------------------

c Declare variables ----------------------------------------------------
      parameter (nmax=101)
      real r(1)
      real*8 T(1)
      real*8 Al(nmax), Ad(nmax), Ar(nmax), b(nmax), beta(nmax),gam(nmax)

c Assign constants in the linear equation ------------------------------
      dr = (r_out-r_in)/float(n)
      do i=2, n
        r(i)  = r_in + dr*(i-1)
        Al(i) = -1./2. + r(i)/dr
        Ad(i) = -2.*r(i)/dr
        Ar(i) =  1./2. + r(i)/dr
      enddo

c B.C.: at r=r_in, T=T_in ----------------------------------------------
      r(1)   = r_in
      Ad(1)  = 1.
      b(1)   = T_in
c B.C.: at r=r_out, T=T_out --------------------------------------------
      r(n+1) = r_out
      Ad(n+1)= 1.
      b(n+1) = T_out

c call a subroutine to solve for a banded matrix with the Thomas method
      call TM(n+1, Al, Ad, Ar, b, T, beta, gam)

      end