c
c fusion.f - simulation of the temperature field in an optical fiber
c            that is being locally heated by a time-dependent external
c            source.  The method is finite element solution to the
c            diffusion equation.
c
c                   2                    dT
c             K grad  T(x,y,z,t) - c rho --  =  S(x,y,z,t)
c                                        dt
c
c            where T(x,y,z) is the temperature field inside the fiber,
c            c is the heat capacity of the fiber material, rho is its
c            density, and K is its thermal conductivity.  The field S,
c            which has dimensions power/volume, is the heat source.
c            The boundary conditions stipulate that the temperature
c            inside the fiber is clamped to room temperature at the
c            extreme ends of the fiber, and that convective cooling
c            and radiative cooling occurs at the boundaries in x,y.
c
c            Convective heat transport at the fiber surface is modeled
c            using Newton's Law of Cooling.
c
c              flux = h (T  - T )
c                         s    b
c
c            where T_s is the temperature of the solid at the surface
c            and T_b is the temperature of the bulk gas far from the
c            surface.  The convective heat transfer coefficient h is
c            sensitive to airflow conditions (see convection.pdf in
c            this directory) but for free convection at sea level the
c            value of 10 W/m^2/K is an often quoted value.

c            Radiative heat transport at the fiber surface is modeled
c            by the Stefan-Bolzmann equation.
c
c                                 4    4
c              flux = ems sigma (T  - T )
c                                 s    b
c
c            where sigma is the Stefan-Boltzmann constant, ems is the
c            emissivity of polystyrene, and T_s,T_b are the same as
c            defined above for convection.
c
c author: Richard Jones at uconn.edu
c version: Oct. 1, 2009

      subroutine fusion(T0)
*	implicit none
      real T0
      include 'fusion.inc'
      write(6,*) 'fusion - simulation of local heating in an',
     +           ' optical fiber'
      write(6,*) '  initializing fiber at uniform temperature ',
     +           T0,' K'
      call reset(T0)
      end

      subroutine stepper(step)
*	implicit none
      real step
      include 'fusion.inc'
      real Source
      external Source
      integer ix,iy,iz
      real x,y,z
      real x0,y0,z0
      real dx,dy,dz
      real dx2,dy2,dz2
      real LaplacianT,dT
      real Fconv,Frad
      call increment_time(step)
      dx = fiberdim(1)/ndivx
      dy = fiberdim(2)/ndivy
      dz = fiberdim(3)/ndivz
      x0 = -(fiberdim(1)+dx)/2
      y0 = -(fiberdim(2)+dy)/2
      z0 = -(fiberdim(3)+dz)/2
      dx2 = dx**2
      dy2 = dy**2
      dz2 = dz**2
       
c  do the internal cells first, they are easier
      do ix=2,ndivx-1
        x = ix*dx+x0
        do iy=2,ndivy-1
          y = iy*dy+y0
          do iz=2,ndivz-1
            z = iz*dz+z0
            LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                    + Tfield(ix-1,iy,iz,nbuf)
     +                    - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +                   (Tfield(ix,iy+1,iz,nbuf) 
     +                    + Tfield(ix,iy-1,iz,nbuf)
     +                    - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +                   (Tfield(ix,iy,iz+1,nbuf) 
     +                    + Tfield(ix,iy,iz-1,nbuf)
     +                    - 2*Tfield(ix,iy,iz,nbuf))/dz2
            dT = (Source(x,y,z) + kcond*LaplacianT)*tstep/(cheat*rho)
            Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          enddo
        enddo
      enddo

c  next do the thermal clamps at the ends of the fiber
      do ix=2,ndivx-1
        x = ix*dx+x0
        do iy=2,ndivy-1
          y = iy*dy+y0
          iz = 1
          z = iz*dz+z0
          LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                  + Tfield(ix-1,iy,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy+1,iz,nbuf) 
     +                  + Tfield(ix,iy-1,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tfield(ix,iy,iz+1,nbuf) 
     +                  + Tamb
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          dT = (Source(x,y,z) + kcond*LaplacianT)*tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qcond = Qcond+kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                       *tstep*dx*dy
          iz = ndivz
          z = iz*dz+z0
          LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                  + Tfield(ix-1,iy,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy+1,iz,nbuf) 
     +                  + Tfield(ix,iy-1,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tamb
     +                  + Tfield(ix,iy,iz-1,nbuf) 
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          dT = (Source(x,y,z) + kcond*LaplacianT)*tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qcond = Qcond+kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                       *tstep*dx*dy
        enddo
      enddo

c  apply the boundary conditions on the fiber surfaces (-x,+x)
      do iy=2,ndivy-1
        y = iy*dy+y0
        do iz=2,ndivz-1
          z = iz*dz+z0
          ix = 1
          x = ix*dx+x0
          LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                  - Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy+1,iz,nbuf) 
     +                  + Tfield(ix,iy-1,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tfield(ix,iy,iz+1,nbuf) 
     +                  + Tfield(ix,iy,iz-1,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
          Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
          dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +         *tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qconv = Qconv + Fconv*tstep*dy*dz
          Qrad = Qrad + Frad*tstep*dy*dz
          ix = ndivx
          x = ix*dx+x0
          LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +                  - Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy+1,iz,nbuf) 
     +                  + Tfield(ix,iy-1,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tfield(ix,iy,iz+1,nbuf) 
     +                  + Tfield(ix,iy,iz-1,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
          Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
          dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +         *tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qconv = Qconv + Fconv*tstep*dy*dz
          Qrad = Qrad + Frad*tstep*dy*dz
        enddo
      enddo

c  apply the boundary conditions on the fiber surfaces (-y,+y)
      do ix=2,ndivx-1
        x = ix*dx+x0
        do iz=2,ndivz-1
          z = iz*dz+z0
          iy = 1
          y = iy*dy+y0
          LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                  + Tfield(ix-1,iy,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy+1,iz,nbuf) 
     +                  - Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tfield(ix,iy,iz+1,nbuf) 
     +                  + Tfield(ix,iy,iz-1,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
          Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
          dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +         *tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qconv = Qconv + Fconv*tstep*dx*dz
          Qrad = Qrad + Frad*tstep*dx*dz
          iy = ndivy
          y = iy*dy+y0
          LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                  + Tfield(ix-1,iy,iz,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +                 (Tfield(ix,iy-1,iz,nbuf) 
     +                  - Tfield(ix,iy,iz,nbuf))/dy2  +
     +                 (Tfield(ix,iy,iz+1,nbuf) 
     +                  + Tfield(ix,iy,iz-1,nbuf)
     +                  - 2*Tfield(ix,iy,iz,nbuf))/dz2
          Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
          Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
          dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +         *tstep/(cheat*rho)
          Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
          Qconv = Qconv + Fconv*tstep*dx*dz
          Qrad = Qrad + Frad*tstep*dx*dz
        enddo
      enddo

c  apply the boundary conditions on the 4 edges (-x,-y),(-x,+y),(+x,+y),(+x,-y)
      do iz=2,ndivz-1
        z = iz*dz+z0
        ix = 1
        x = ix*dx+x0
        iy = 1
        y = iy*dy+y0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tfield(ix,iy,iz-1,nbuf)
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z)+kcond*LaplacianT-(Fconv+Frad)*(1/dx+1/dy))
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
        Qrad = Qrad + Frad*tstep*dz*(dx+dy)
        iy = ndivy
        y = iy*dy+y0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy-1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tfield(ix,iy,iz-1,nbuf)
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z)+kcond*LaplacianT-(Fconv+Frad)*(1/dx+1/dy))
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
        Qrad = Qrad + Frad*tstep*dz*(dx+dy)
        ix = ndivx
        x = ix*dx+x0
        LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy-1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tfield(ix,iy,iz-1,nbuf)
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z)+kcond*LaplacianT-(Fconv+Frad)*(1/dx+1/dy))
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
        Qrad = Qrad + Frad*tstep*dz*(dx+dy)
        iy = 1
        y = iy*dy+y0
        LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tfield(ix,iy,iz-1,nbuf)
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z)+kcond*LaplacianT-(Fconv+Frad)*(1/dx+1/dy))
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
        Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      enddo

c  apply the boundary conditions on the 8 edges at the ends
      do ix=2,ndivx-1
        x = ix*dx+x0
        iy = 1
        y = iy*dy+y0
        iz = 1
        z = iz*dz+z0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                + Tfield(ix-1,iy,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tamb
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dx
        Qrad = Qrad + Frad*tstep*dz*dx
        iz = ndivz
        z = iz*dz+z0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                + Tfield(ix-1,iy,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tamb
     +                + Tfield(ix,iy,iz-1,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dx
        Qrad = Qrad + Frad*tstep*dz*dx
        iy = ndivy
        y = iy*dy+y0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                + Tfield(ix-1,iy,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy-1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tamb
     +                + Tfield(ix,iy,iz-1,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dx
        Qrad = Qrad + Frad*tstep*dz*dx
        iz = 1
        z = iz*dz+z0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                + Tfield(ix-1,iy,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy-1,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tamb
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dy)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dx
        Qrad = Qrad + Frad*tstep*dz*dx
      enddo
      do iy=2,ndivy-1
        y = iy*dy+y0
        ix = 1
        x = ix*dx+x0
        iz = 1
        z = iz*dz+z0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                + Tfield(ix,iy-1,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tamb
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dy
        Qrad = Qrad + Frad*tstep*dz*dy
        iz = ndivz
        z = iz*dz+z0
        LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                + Tfield(ix,iy-1,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tamb
     +                + Tfield(ix,iy,iz-1,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dy
        Qrad = Qrad + Frad*tstep*dz*dy
        ix = ndivx
        x = ix*dx+x0
        LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                + Tfield(ix,iy-1,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tamb
     +                + Tfield(ix,iy,iz-1,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dy
        Qrad = Qrad + Frad*tstep*dz*dy
        iz = 1
        z = iz*dz+z0
        LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +                - Tfield(ix,iy,iz,nbuf))/dx2  +
     +               (Tfield(ix,iy+1,iz,nbuf) 
     +                + Tfield(ix,iy-1,iz,nbuf) 
     +                - 2*Tfield(ix,iy,iz,nbuf))/dy2  +
     +               (Tfield(ix,iy,iz+1,nbuf) 
     +                + Tamb
     +                - 2*Tfield(ix,iy,iz,nbuf))/dz2
        Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
        Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
        dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)/dx)
     +       *tstep/(cheat*rho)
        Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
        Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                  *tstep*dx*dy
        Qconv = Qconv + Fconv*tstep*dz*dy
        Qrad = Qrad + Frad*tstep*dz*dy
      enddo

c  now only the eight corners are left
      ix = 1
      x = ix*dx+x0
      iy = 1
      y = iy*dy+y0
      iz = 1
      z = iz*dz+z0
      LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy+1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tfield(ix,iy,iz+1,nbuf) 
     +              + Tamb
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      ix = ndivx
      x = ix*dx+x0
      LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy+1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tfield(ix,iy,iz+1,nbuf) 
     +              + Tamb
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      iy = ndivy
      y = iy*dy+y0
      LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy-1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tfield(ix,iy,iz+1,nbuf) 
     +              + Tamb
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      ix = 1
      x = ix*dx+x0
      LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy-1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tfield(ix,iy,iz+1,nbuf) 
     +              + Tamb
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      iz = ndivz
      z = iz*dz+z0
      ix = 1
      x = ix*dx+x0
      iy = 1
      y = iy*dy+y0
      LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy+1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tamb
     +              + Tfield(ix,iy,iz-1,nbuf) 
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      ix = ndivx
      x = ix*dx+x0
      LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy+1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tamb
     +              + Tfield(ix,iy,iz-1,nbuf) 
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      iy = ndivy
      y = iy*dy+y0
      LaplacianT = (Tfield(ix-1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy-1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tamb
     +              + Tfield(ix,iy,iz-1,nbuf) 
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      ix = 1
      x = ix*dx+x0
      LaplacianT = (Tfield(ix+1,iy,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dx2  +
     +             (Tfield(ix,iy-1,iz,nbuf) 
     +              - Tfield(ix,iy,iz,nbuf))/dy2  +
     +             (Tamb
     +              + Tfield(ix,iy,iz-1,nbuf) 
     +              - 2*Tfield(ix,iy,iz,nbuf))/dz2
      Fconv = hconv*(Tfield(ix,iy,iz,nbuf)-Tamb)
      Frad = ems*SBconst*(Tfield(ix,iy,iz,nbuf)**4-Tamb**4)
      dT = (Source(x,y,z) + kcond*LaplacianT - (Fconv+Frad)*(1/dx+1/dy))
     +     *tstep/(cheat*rho)
      Tfield(ix,iy,iz,obuf) = Tfield(ix,iy,iz,nbuf) + dT
      Qcond = Qcond + kcond*(Tfield(ix,iy,iz,nbuf)-Tamb)/dz
     +                *tstep*dx*dy
      Qconv = Qconv + Fconv*tstep*dz*(dx+dy)
      Qrad = Qrad + Frad*tstep*dz*(dx+dy)
      end

      real function Source(x,y,z)
*	implicit none
      real x,y,z
      include 'fusion.inc'
      real dist2
      real time0,sigmat
      data time0,sigmat/1.,0.1/
      if (Snorm .eq. 0) then
        Snorm = (2*pi)**1.5 * Ssize(1)*Ssize(2)*Ssize(3)
      endif
      dist2 = ((x-Scent(1))/Ssize(1))**2 +
     +        ((y-Scent(2))/Ssize(2))**2 +
     +        ((z-Scent(3))/Ssize(3))**2
      Source = Spower*exp(-dist2/2)/Snorm
      Source = Source*exp(-0.5*((time-time0)/sigmat)**2)
     +         / (sqrt(2*pi)*sigmat)
      end

      subroutine setPower(watts)
*	implicit none
      real watts
      include 'fusion.inc'
      integer ix,iy,iz
      real x,y,z
      real dx,dy,dz
      real x0,y0,z0
      real totalPower
      if (watts .eq. 0) then
        write(6,*) 'current value of Spower =',Spower
      else
        Spower = watts
      endif
      dx = fiberdim(1)/ndivx
      dy = fiberdim(2)/ndivy
      dz = fiberdim(3)/ndivz
      x0 = -(fiberdim(1)+dx)/2
      y0 = -(fiberdim(2)+dy)/2
      z0 = -(fiberdim(3)+dz)/2
      totalPower = 0
      do ix=1,ndivx
        x = ix*dx+x0
        do iy=1,ndivy
          y = iy*dy+y0
          do iz=1,ndivz
            z = iz*dz+z0
            totalPower = totalPower+Source(x,y,z)
          enddo
        enddo
      enddo
*     print *, 'total power is',totalPower*dx*dy*dz,' W'
      end
      
      subroutine reset(T)
*	implicit none
      real T
      include 'fusion.inc'
      integer ix,iy,iz
      real T0
      if (T .lt. 0) then
        T0 = Tamb
      else
        T0 = T
      endif
      time = 0
      nstep = 0
      do ix=1,ndivx
        do iy=1,ndivy
          do iz=1,ndivz
            Tfield(ix,iy,iz,1) = T0
            Tfield(ix,iy,iz,2) = T0
          enddo
        enddo
      enddo
      Qcond = 0
      Qconv = 0
      Qrad = 0
      end

      subroutine increment_time(step)
*	implicit none
      real step
      include 'fusion.inc'
      time = time + step
      tstep = step
      nstep = nstep+1
      nbuf = 3-nbuf
      obuf = 3-obuf
      end

      subroutine makehist(id)
*	implicit none
      integer id
      include 'fusion.inc'
      integer ix,iy,iz
      integer ix0,iy0,iz0
      character*80 title
      logical hexist
      external hexist
      real cutx(ndivz,ndivy)
      real cuty(ndivz,ndivx)
      real cutz(ndivx,ndivy)
      ix0 = (ndivx+1)/2
      do iz=1,ndivz
        do iy=1,ndivy
          cutx(iz,iy) = Tfield(ix0,iy,iz,obuf)
        enddo
      enddo
      iy0 = (ndivy+1)/2
      do iz=1,ndivz
        do ix=1,ndivx
          cuty(iz,ix) = Tfield(ix,iy0,iz,obuf)
        enddo
      enddo
      iz0 = (ndivz+1)/2
      do ix=1,ndivx
        do iy=1,ndivy
          cutz(ix,iy) = Tfield(ix,iy,iz0,obuf)
        enddo
      enddo

      if (hexist(id)) then
        call hdelet(id)
      endif
      write(title,*) 'slice through x=0 at time',time
      call hbook2(id,title,ndivz,-fiberdim(3)/2,fiberdim(3)/2,
     +                     ndivy,-fiberdim(2)/2,fiberdim(2)/2,0)
      call hpak(id,cutx)

      if (hexist(id+1)) then
        call hdelet(id+1)
      endif
      write(title,*) 'slice through y=0 at time',time
      call hbook2(id+1,title,ndivz,-fiberdim(3)/2,fiberdim(3)/2,
     +                       ndivx,-fiberdim(1)/2,fiberdim(1)/2,0)
      call hpak(id+1,cuty)

      if (hexist(id+2)) then
        call hdelet(id+2)
      endif
      write(title,*) 'slice through z=0 at time',time
      call hbook2(id+2,title,ndivx,-fiberdim(1)/2,fiberdim(1)/2,
     +                       ndivy,-fiberdim(2)/2,fiberdim(2)/2,0)
      call hpak(id+2,cutz)

c     print *, 'the 8 corners are',
c    +         (((Tfield(i,j,k,obuf),
c    +            i=1,ndivx,ndivx-1),
c    +              j=1,ndivy,ndivy-1),
c    +                k=1,ndivz,ndivz-1)
      end

      subroutine getQ(Q)
      real Q(3)
*	implicit none
      include 'fusion.inc'
      Q(1) = Qcond
      Q(2) = Qconv
      Q(3) = Qrad
      end