c------------------------------------------------------
c
c     Forcem user subroutine set up to detect contact.
c     When element is in contact, pressure is not applied.
c     When element is not in contact, the pressure given 
c     by the model is applied.
c
c     Only tested for 4-noded 2D continuum elements.      
c
c------------------------------------------------------
      subroutine forcem(press,th1,th2,nn,n)
      implicit real*8 (a-h,o-z)                                               dp
      dimension n(7)
c* * * * * *
c     defined pressure on an element.
c     press  total pressure            th1    coordinate
c     th2    coordinate                nn     integration point number
c     n(1)   element number            n(2)   load type
c     n(3)   integration point number  n(4)   not used
c     n(5)   dist load index           n(6)   not used
c     n(7)   internal element number
c* * * * * *
      include '../common/form'
      include '../common/array4'
      include '../common/concom'
      include '../common/space'
      include '../common/blnk'
      common /mydata/press_val
      logical stays_off
      save stays_off
      parameter (MAX_ELS=1000, NODES_PER_EL=4)
      dimension nod(NODES_PER_EL), icheck(NODES_PER_EL)
      dimension stays_off(MAX_ELS), press_val(MAX_ELS)
      equivalence (nod(1),lm(1))
c-----------------------------------------------------
      do 5 k=1,NODES_PER_EL
5     icheck(k) = 0
c
c-----------------------------------------------------
c     Contact checks
c
      do 100 i1=1,numdie
      nbct = ints(inbct+i1-1)
      j1 = itouch + (i1 - 1)*nbcn
      j2 = inf    + (i1 - 1)*nbcn
        do 90 i2=1,nbct
        idie = ints(j1+i2-1)
        inod = ints(j2+i2-1)
        if(idie.eq.0) go to 90
c******   don't count def contact
        if(idie.eq.1) go to 90
           do 80 k=1,NODES_PER_EL
              if(inod.eq.nod(k)) icheck(k) = 1
   80      continue
   90   continue
  100 continue
c
      ichk = 0
      do 120 i1=1,NODES_PER_EL
      ichk = ichk + icheck(i1)
  120 continue
c
c     End of contact check
c------------------------------------------------------
c
      if(ichk.eq.2) then
         factor = 0.0
         stays_off(n(1))=.true.
         write(6,1001) inc,n(1),nn,factor,n(2),(icheck(k),k=1,4)
      else 
         if(stays_off(n(1))) then
            factor = 0.0
            write(6,1002) inc,n(1),nn,factor,n(2),(icheck(k),k=1,4)
         else
            factor = 1.0
            write(6,1003) inc,n(1),nn,factor,n(2),(icheck(k),k=1,4)
         end if 
      end if
 1001 format('B1: inc,n(1),nn,factor,n(2),side=',3i5,1pe10.3,5i5)
 1002 format('B2: inc,n(1),nn,factor,n(2),side=',3i5,1pe10.3,5i5)
 1003 format('B3: inc,n(1),nn,factor,n(2),side=',3i5,1pe10.3,5i5)
c
      press = factor * press
      press_val(n(1)) = press
c      write(6,1000) inc,n(1),nn,press,n(2),(icheck(k),k=1,4)
c1000  format(' inc,n(1),nn,press,n(2),side=',3i5,1pe10.3,5i5)
      return
      end
c
c
      subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,
     * nshear,jpltcd)
c* * * * * *
c     select a variable contour plotting (user subroutine).
c     v      variable                    s      stress array
c     sp     stresses in preferred dir.  etot   total strain (generalized)
c     eplas  total plastic strain        ecreep total creep strain
c     t      current temperature         m      element number
c     nn     integration point number    layer  layer number
c     ndi    # of direct stress comp     nshear # of shear stress components
c* * * * * *
      implicit real*8 (a-h,o-z)                                               dp      dimension s(1),etot(1),eplas(1),ecreep(1),sp(1)
      common /mydata/press_val
      parameter (MAX_ELS=1000)
      dimension press_val(MAX_ELS)
      v=press_val(m)
      return
      end