c     -------------------------------------------------------------------------
c     User subroutine VUSDFLD for obtaining volumetric strain rate and computing porosity
c     -------------------------------------------------------------------------
      subroutine vusdfld(
c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr, 
     *   jElemUid, kIntPt, kLayer, kSecPt, 
     *   stepTime, totalTime, dt, cmname, 
     *   coordMp, direct, T, charLength, props, 
     *   stateOld, 
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension jElemUid(nblock), coordMp(nblock,*), 
     *          direct(nblock,3,3), T(nblock,3,3), 
     *          charLength(nblock), props(nprops), 
     *          stateOld(nblock,nstatev), 
     *          stateNew(nblock,nstatev),
     *          field(nblock,nfieldv)
      character*80 cmname
c
      parameter (zero=0.d0, one=1.d0)
      parameter (l_sdv_Porosity = 1,
     *           l_sdv_VoidR    = 2,
     *           l_sdv_POR      = 3,
     *           l_sdv_EvolR    = 4)
      parameter (l_fv_Porosity = 1,
     *           l_fv_POR     = 2)
c     read properties
      alpha = props(1)
      rKs = props(2)
      rKsInv = one / rKs
c Compute volumetric strain rate
      call getEvolRate( nblock, ndir, stateNew(1,l_sdv_EvolR) )
c Update porosity and set field variable
      do k = 1, nblock
c          
         stateNew(k,l_sdv_POR) = field(k,l_fv_POR) 
c
         stateNew(k,l_sdv_Porosity) = stateOld(k,l_sdv_Porosity) +
     *        (stateOld(k,l_sdv_Porosity)-alpha)*(-dt*stateNew(k,l_sdv_EvolR))
         stateNew(k,l_sdv_Porosity) = max(zero,stateNew(k,l_sdv_Porosity))
c     Contribution due to pore pressure change is one inc behind (added in VUEXPAN)
c     *        + (porosity-alpha)*(three*CTE*TEMP_Inc - ksInv * POR_Inc)
         field(k,l_fv_Porosity) = stateNew(k,l_sdv_Porosity)
      end do
c
      return
      end
c     -------------------------------------------------------------------------
c     User subroutine VUEXPAN for defining thermal and field expansion
c     -------------------------------------------------------------------------
      subroutine vuexpan (
C Read only (unmodifiable)variables -
     * nblock, nDir, nShr, nExpanType,
     * nElem, nIntPt, nLayer, nSectPt,
     * stepTime, totalTime, dt, cmname, 
     * nstatev, nfieldv, nprops, props,
     * tempOld, tempNew, fieldOld, fieldNew,
     * stateOld, 
C Write only (modifiable) variable -
     * stateNew, strainThInc, dStrainTherDT )
C
      include 'vaba_param.inc'
C
      dimension strainThInc(nblock,nDir+nShr),
     * dStrainTherDT(nblock,nDir+nShr),
     * nElem(nblock), props(nprops),
     * tempOld(nblock), tempNew(nblock),
     * fieldOld(nblock,nfieldv),
     * fieldNew(nblock,nfieldv),
     * stateOld(nblock,nstatev),
     * stateNew(nblock,nstatev)
C
      character*80 cmname
C
      parameter (zero=0.d0, one=1.d0, three=3.d0, third=1.d0/3.d0)
      parameter (l_sdv_Porosity = 1,
     *           l_sdv_VoidR    = 2,
     *           l_sdv_POR      = 3,
     *           l_sdv_EvolR    = 4)
      parameter (l_fv_Porosity = 1,
     *           l_fv_POR     = 2)
C
      alpha = props(1)
      rKs = props(2)
      alphaTh = props(3)
      alphaPor = - third / rKs
C
      do k = 1, nblock
         dtemp = tempNew(k) - tempOld(k)
         dPOR = fieldNew(k,l_fv_POR)-fieldOld(k,l_fv_POR)
         deth = alphaTh * dtemp + alphaPor * dPOR
         
         strainThInc(k,1) = deth
         dStrainTherDT(k,1) = alphaTh

         stateNew(k,l_sdv_Porosity) = stateNew(k,l_sdv_Porosity) +
     *        (stateNew(k,l_sdv_Porosity)-alpha)*deth*three
         stateNew(k,l_sdv_Porosity) = max(zero,stateNew(k,l_sdv_Porosity))
         stateNew(k,l_sdv_VoidR) = stateNew(k,l_sdv_Porosity) /
     *        ( one - stateNew(k,l_sdv_Porosity) )
      end do

      return
      end
c     -------------------------------------------------------------------------
c     Utilities
c     -------------------------------------------------------------------------
      subroutine getEvolRate ( nblock, ndir, evolRate )
      include 'vaba_param.inc'
      dimension evolRate(nblock)
c     
      parameter( nrData=6, zero = 0.d0 )
      character*3 cData(maxblk*nrData)
      dimension jData(maxblk*nrData)
      dimension rData(maxblk*nrData)
c     Access strain rate data (ER)
      jStatus = 1
      call vgetvrm( 'ER', rData, jData, cData, jStatus )
      if( jStatus .ne. 0 ) then
         call xplb_abqerr(-2,'Utility routine VGETVRM failed '//
     *        'to get variable.',0,zero,' ')
         call xplb_exit
      end if
c     Compute gammaDot
      if (ndir.eq.3.) then
         call evolRate33( nblock, rData, evolRate )
      end if
c
      return 
      end
c
      subroutine evolRate33 ( nblock, strainRate, evolRate )
      include 'vaba_param.inc'
      dimension strainRate(nblock,6)
      dimension evolRate(nblock)
      do k = 1, nblock
         evolRate(k) = strainRate(k,1) + strainRate(k,2) + strainRate(k,3)
      end do
      return 
      end