subroutine uinter(stress,ddsddr,amki,amski,flux,ddfddt,ddsddt,
     1     ddfddr,statev,sed,sfd,spd,svd,scd,pnewdt,rdisp,drdisp,
     2     temp,dtemp,predef,dpred,time,dtime,freqr,ciname,slname,
     3     msname,props,coords,aLocalDir,drot,area,chrLngth,node,ndir,
     4     nstatv,npred,nprops,mcrd,kstep,kinc,kit,linper,lOpenClose,
     5     lState,lSdi,lPrint)
c 
      include 'aba_param.inc'
c 
      dimension stress(ndir),ddsddr(ndir,ndir),flux(2),ddfddt(2,2),
     $     ddsddt(ndir,2),ddfddr(2,ndir),statev(nstatv),rdisp(ndir),
     $     drdisp(ndir),temp(2),dtemp(2),predef(2,npred),dpred(2,npred),
     $     time(2),props(nprops),coords(mcrd),aLocalDir(mcrd,mcrd),
     $     drot(2,2),amki(ndir,ndir),amski(ndir,ndir)
      character*80 ciname,slname,msname
      parameter(toler = 1.D-12, zero=0.d0, one=1.d0, two=2.d0,
     $     half = one/two, fourth=0.25D0, eighth=0.125D0)
c 
c             User defined interfacial constitutive behavior.
c
c
c Local variables
      p_t = zero
      tau1_t = zero
      tau2_t = zero
      dp = zero
c Material properties
c
      ekbar = props(1) ! Modulus in the normal direction
c
      xk=props(2)+props(3)*(temp(1)+temp(2))/two

      if ( rdisp(1) .ge. toler) then
c
        lOpenClose = 1
c     Compute P_(t+dt) and T_crit
c        
        p_t = stress(1)
        tau1_t = stress(2)
        tau2_t = stress(3)
        dp = ekbar*drdisp(1)
        stress(1) = stress(1) + dp
        stress(2) = drdisp(2)/(xk*dtime)
        stress(3) = drdisp(3)/(xk*dtime)
c
        ddsddr(1,1) = ekbar
        ddsddr(2,2)=one/(xk*dtime)
        ddsddr(3,3)=one/(xk*dtime)
        ddsddr(1,2)=zero
        ddsddr(2,1)=zero
        ddsddr(1,3)=zero
        ddsddr(3,1)=zero
        ddsddr(3,2)=zero
        ddsddr(2,3)=zero  
c
        ddsddt(1,1)=zero
        ddsddt(1,2)=zero
        ddsddt(2,1)=-0.0001*drdisp(2)/(dtime*two*xk**two)
        ddsddt(2,2)=ddsddt(2,1)
        ddsddt(3,1)=-0.0001*drdisp(3)/(dtime*two*xk**two)
        ddsddt(3,2)=ddsddt(3,1)
c
        sfd = stress(2)*drdisp(2) + stress(3)*drdisp(3)
c
        heatGen = half*sfd
        flux(1) = half*heatGen/dtime
        flux(2) = flux(1)
c
        ddfddr(1,1) = zero
        ddfddr(2,1) = zero
        ddfddr(1,2) = fourth*stress(2)/dtime
        ddfddr(2,2) = ddfddr(1,2)
        ddfddr(1,3) = fourth*stress(3)/dtime
        ddfddr(2,3) = ddfddr(1,3)
c
        coeff = eighth*props(3)/(xk*xk*dtime*dtime)
        ddfddt(1,1) = coeff*(drdisp(2)**2 + drdisp(3)**2)
        ddfddt(1,2) = ddfddt(1,1)
        ddfddt(2,1) = ddfddt(1,1)
        ddfddt(2,2) = ddfddt(1,1)
c
      else
c
        lOpenClose = 0
        sfd = zero
c     
      end if
c     
      return
      end