ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c		E C 2 4 D F D 1
c
c          FINITE DIFFERENCE SCHEME for 1-D TRANSIENT
c          HEAT CONDUCTION with FILM BOUNDARY CONDITION
c
c          OUTPUT FILE CREATED IS: history.dat (distance vs. temp)
c          AT THE END OF THE TRANSIENT ANALYSIS
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      parameter (npoints = 100)
      dimension temp(npoints+1)
c
c                give the material constants
c
      conductk = 1.40
      spheat   = 260.0
      density  = 7800.0
c
c                linear variation of film coefficient
c                h = h0 + h1 * temp
c
      h0       = 10.0
      h1       = 0.02
c
      alpha    = conductk / (spheat * density)
      onealpha = 1.0 / alpha
c
c                compute stable time increment
c
      tmax      = 3600.0
      nmax      = 360000
      totlength = 0.1 
      ds        = totlength / float(npoints)
      dtime     = tmax / float(nmax)
      onestable = alpha * dtime / (ds * ds)
c
c                initialize temp
c
      do i = 1,101
        temp(i) = 0.0
      end do
c
c                open data files
c
      open(unit=8,file='history.dat')
c
      time = 0.0
c
c                "it" do loop for time-space
c
      do it = 1,nmax+1
c
c                film boundary condition
c
        tinf    = 100.0 * (1.0 + time / tmax)
        filmcon = h0 + h1 * temp(1)
        hdsk    = filmcon * ds / conductk
        temp(1) = 2.0 * onestable * (temp(2) + hdsk * tinf)
     $       + (1.0 - 2.0 * onestable * (hdsk + 1.0)) * temp(1)
c
c                "i" do loop for spatial integration
c
        do i = 2,100
          im1     = i - 1
          ip1     = i + 1
          temp(i) = onestable * (temp(im1) + temp(ip1))
     $         + (1.0 - 2.0 * onestable) * temp(i)
        end do

        temp(101) = onestable * 2.0 * temp(100)
     $       + (1.0 - 2.0 * onestable) * temp(101)

        if (it.eq.nmax) then
          do iw = 1,101
            dist = 0.1 - float(iw - 1) * ds
            write(8,222) dist, temp(iw)
          end do
        end if

        time = time + dtime
      end do
c
 222  format(2f10.4)
c
      stop
      end