C
C User subroutine vumatht
      subroutine vumatht (
C Read only (unmodifiable) variables -
     *     nblock, nElem, nIntPt, nLayer, nSectPt, 
     *     ntgrad, nstatev, nfieldv, nprops,  
     *     cmname, stepTime, totalTime, dt,  
     *     coordMp, density, props,  
     *     tempOld, fieldOld, stateOld, enerThermOld, 
     *     tempNew, tempgradNew, fieldNew, 
C Write only (modifiable) variables -
     *     stateNew, fluxNew, enerThermNew, dEnerThDTemp, condEff )
C
      include 'vaba_param.inc'
C
      dimension nElem(nblock)
C
      dimension coordMp(nblock,*), density(nblock), props(nprops),
     *    tempOld(nblock), fieldOld(nblock, nfieldv), 
     *    stateOld(nblock, nstatev), enerThermOld(nblock),
     *    tempNew(nblock), tempgradNew(nblock, ntgrad), 
     *    fieldNew(nblock, nfieldv), stateNew(nblock, nstatev),
     *    fluxNew(nblock, ntgrad), enerThermNew(nblock), 
     *    dEnerThDTemp(nblock,2), condEff(nblock)
C
      character*80 cmname
*
      parameter ( zero=0.d0 )
      parameter ( tempZero = 0.d0, tol = 1.d-6)
* properties
* 1. conductivity
* 2. specific heat
* 3. latent heat per unit mass
* 4. solidus temperature 
* 5. liquidus temperature
*
      cond = props(1)
      specHeat = props(2)
      uLatnHeat = props(3)
      tempSol = props(4)
      tempLiq = props(5)
*
      cLatnHeat = uLatnHeat/(tempLiq-tempSol)
*
      do km = 1, nblock
*
*        update heat flux vector
         do i = 1, ntgrad
            fluxNew(km, i) = -cond*tempgradNew(km,i)
         end do
         condEff(km) = cond
*
*        update internal thermal energy and its derivatives
         dtemp = tempNew(km) - tempOld(km)
         if (totalTime .eq. zero) then
            dEnerTh = specHeat*(tempNew(km) - tempZero)
         else
            dEnerTh = specHeat*dtemp
         end if
*
         if (tempOld(km) .ge. tempLiq) then
            uLatnOld = uLatnHeat
         else if (tempOld(km) .ge. tempSol) then
            uLatnOld = cLatnHeat*(tempOld(km) - tempSol)
         else
            uLatnOld = zero
         end if
         if (tempNew(km) .ge. tempLiq) then
            uLatnNew = uLatnHeat
         else if (tempNew(km) .ge. tempSol) then
            uLatnNew = cLatnHeat*(tempNew(km) - tempSol)
         else
            uLatnNew = zero
         end if
         if (abs(dtemp) .gt. tol) then
            dULatnDTemp = (uLatnNew-uLatnOld)/dtemp
         else
            dULatnDTemp = zero
         end if
         dEnerTh = dEnerTh + uLatnNew - uLatnOld
*
         enerThermNew(km) = enerThermOld(km) + dEnerTh
         dEnerThDTemp(km,1) = specHeat + dULatnDTemp
         dEnerThDTemp(km,2) = dULatnDTemp
      end do
*
      return
      end