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)
c 
c             User defined interfacial constitutive behavior.
c             Exponential press-clearance relation, coulomb friction
c
c Local variables
c
      dp = zero
      taucrit = zero
      stiff = zero
      dtau = zero
      tauelas = zero
      dslip = zero
      gbarelas = zero
      taupr = zero
c
c Material properties
c
      xmu = props(1) 
      gcrit = props(2) 
      c_0 = props(3) 
      p_0 = props(4) 

      h = rdisp(1)
      if ( h .gt. -c_0) then
c    
c    Compute normal stress and Tau_crit
c
         stress(1) = (p_0/(exp(one)-one))*
     $        ( ((h/c_0)+one)*(exp(h/c_0+one)-one) )
         taucrit = xmu*stress(1)         
         stiff = taucrit/gcrit
c
c    Compute tau assuming increment is elastic
c
         tauelas = stiff*(statev(1)+drdisp(2))
c
c    Check whether slip is "elastic" or "plastic"
c
         if ( abs(tauelas).lt. taucrit) then
c
c       slip is "elastic"
c
            stress(2) = tauelas
            if (stiff.ge.toler) then
               statev(1) = statev(1)+drdisp(2)
               sed = half*stress(2)*stress(2)/stiff
            else
               sed = zero
            end if
            dslip = zero
            statev(2) = statev(2)+dslip
c       Jacobian
            ddsddr(1,1) = (p_0/(exp(one)-one))*
     $           ( ((h/c_0+two)*exp(h/c_0+one))/c_0 - one/c_0 )
            ddsddr(2,1) = statev(1)*xmu*ddsddr(1,1)/gcrit
            ddsddr(2,2) = stiff
         else
c
c       slip is "plastic"            
c
            gtplusdt = rdisp(2)
            stress(2) = gtplusdt*taucrit/abs(gtplusdt)
c       Jacobian
            ddsddr(1,1) = (p_0/(exp(one)-one))*
     $           ( ((h/c_0+two)*exp(h/c_0+one))/c_0 - one/c_0 )

            ddsddr(2,1) = (gtplusdt/abs(gtplusdt))*xmu*ddsddr(1,1)
c       Plastic slip
            gbarelas = statev(1)
            taupr = stiff*(gbarelas+drdisp(2))
            dslip = (taupr - stress(2))/stiff
            gbarelas = stress(2)/stiff

            statev(1) = gbarelas
            statev(2) = statev(2)+dslip            
            sed = half*stress(2)*stress(2)/stiff
         end if
 900     continue
      else
c
c    Gap is open
c
         if (h.eq.-c_0) then            
            ddsddr(2,2) = 1.D6
         end if
         dslip = zero
         statev(1) = zero
         statev(2) = statev(2)+dslip
         sed = zero
c    
      end if

      
      return
      end