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