Listing 1:  FORTRAN Program fdtemp.for, implicit method finite difference method to solve [12]

c     fdtemp

      implicit double precision (a-h, k, o-z)
      parameter (iz=100)
      dimension aa(0:iz), bb(0:iz), cc(0:iz), dd(0:iz)
      dimension tmj(0:iz), tm(0:iz)

      open (1, file='outtem.dat', status='unknown')
      zl = 2.d0
      tl = 3600.d0*24.d0
      np = 80
      dt = 1.0d0
      dift = 6.d-7
      t0 = 15.d0
      tnp = 25.d0
      tini = 16.d0
      t = 0.d0
      dz = zl/dble(np)
      rz = dift*dt/(dz*dz)
      rz2 = 1.d0 + 2.d0*rz
      do i = 1, np-1
         tmj(i) = tini
         tm(i)  = tini
      end do
      tm(0) = t0
      tm(np) = tnp
      do while (t .le. tl) 
         t = t + dt
         do i = 1, np-1
            aa(i) = -rz
            bb(i) = rz2
            cc(i) = -rz
            dd(i) = tmj(i)
         end do
         dd(1) = dd(1) + rz*t0
         dd(np-1) = dd(np-1) + rz*tnp
         
         call tridia (np-1, aa, bb, cc, dd)

         do i = 1, np-1
            tm(i) = dd(i)
            tmj(i) = tm(i)
         end do

      end do

      do i = 0, np
         z = dble(i)*dz - zl
         write (1, '(f6.3,f7.3)') z, tm(i)
      end do
      close (1)

      end


      subroutine tridia (nn, a, b, c, d)

      implicit double precision (a-h, o-z)
      parameter (iz=100)
      dimension a(0:iz), b(0:iz), c(0:iz), d(0:iz)

      do n = 2, nn
         r = a(n)/b(n-1)
         b(n) = b(n) - r*c(n-1)
         d(n) = d(n) - r*d(n-1)
      end do
      d(nn) = d(nn)/b(nn)
      do n = nn-1, 1, -1
         d(n) = (d(n) - c(n)*d(n+1))/b(n)
      end do

      return
      end