subroutine ubginc(inc,incsub)
      implicit real*8 (a-h,o-z)
c     usersub that would get called at the beginning of each increment
      include '../common/dimen'
      real*8 press_val
      integer*4 itmp
      integer*4 numnp_old, numel_old
      save numnp_old, numel_old
c
c  Dynamic allocation of memory for forcem and plotv usersubs
c
c     itmp        array for temporary workspace in forcem    
c     pitmp       pointer to above
c     numnp       total number of node points (size of itmp array)
c     press_val   array to hold the pressure value computed in forcem
c     ppress_val  pointer to above
c     numel       total number of elements (size of press_val array)
c     marcfalloc  internal marc routine to allocate memory 
c
      common /mydata1/pitmp
      common /mydata2/ppress_val
      pointer(pitmp, itmp(1))
      pointer(ppress_val, press_val(1))
c
c   Allocate memory for temporary storge and pressure storage
c   Update storage if numnp or numel changes
c
c      write(64,100) inc, numnp, numel
c 100  format(1x,/,' Inc =', i5, '  Numnp =', i5, '  Numel =',i5,/)
      if (inc.eq.0) then
        call marcfalloc(numnp,4,1,pitmp)
        call marcfalloc(numel,8,1,ppress_val)
        numnp_old = numnp
        numel_old = numel
      else
        if (numnp.gt.numnp_old) then
          call marcfalloc(numnp,4,-1,pitmp)
          numnp_old = numnp
          do i=1,numnp
            itmp(i) = 0
          enddo
        endif
        if (numel.gt.numel_old) then
          call marcfalloc(numel,8,-1,ppress_val)
          numel_old = numel
          do i=1,numel
            press_val(i) = 0.0
          enddo
        endif
      endif
      return
      end
c
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/blnk'
      include '../common/concom'
      common /mydata1/pitmp
      common /mydata2/ppress_val
      pointer(pitmp, itmp(1))
      pointer(ppress_val, press_val(1))
c
c   flags for keeping pressure off in special cases
c      logical stays_off
c      save stays_off
c
      logical DEBUG
      parameter (DEBUG=.false.)
      integer inc_cur, ncycle_cur
      save inc_cur, ncycle_cur
      data inc_cur, ncycle_cur/ -1, -1 /
c
c   Contact checking variables
c
c    lm()        connectivity array for this element (in blnk common)
c    load_type   element load type (see Vol. B)
c                expect this to be 3,7,9, or 11 for 4-noded quad elements
c                for side 1,2,3, or 4 respectively
c    in_contact  contact flag for element side
c                = 0  no nodes of this side are in contact - pressure on
c                = 1  one node of this side is in contact  - pressure on
c                = 2  two nodes of this side are in contact - pressure off
c                contact means node is a tied node or retained node
c
c    for each node 
c        icntyp = 0  no contact
c               = 1  node is contacting rigid surface only
c               = 2  node is a tied node
c               = 3  node is a retained node
c               >10  node is a retained node in multiple ties
c
c
      if (DEBUG) then
        if (inc.ne.inc_cur.or.ncycle.ne.ncycle_cur) then
          write(65,900) inc,ncycle
 900      format(1x,/,' *** Increment = ',i5,'   Iteration = ',i5,/)
        endif
        inc_cur  =inc
        ncycle_cur = ncycle
      endif

c------------------------------------------------------------
c   TO BE EDITED BY THE USER
c   List the nodes to be checked for contact,
c   based on element type and load type - see Volume B.
c   The following is for 4-noded quad elements (11,80,3,114,10,82)
c
      load_type = n(2)
        if (load_type .eq. 3) then
          ni1 = 1
          ni2 = 2
        else if (load_type .eq. 7) then
          ni1 = 2
          ni2 = 3
        else if (load_type .eq. 9) then
          ni1 = 3
          ni2 = 4
        else if (load_type .eq. 11) then
          ni1 = 4
          ni2 = 1
        endif
c
c------------------------------------------------------------
c   Begin contact check
c
c   zero out the in_contact flag 
      in_contact = 0
c
c   load temporary workspace with contact info
c
      call setwr2(itmp)
c
c   extract contact info for 1st node of this side
c
      icntyp = itmp(lm(ni1))
      if(icntyp.gt.0) in_contact=in_contact+1
      if (DEBUG) then
        write(65,980) ni1, lm(ni1), icntyp
      endif
c
c
c   extract contact info for 2nd node of this side
c
      icntyp = itmp(lm(ni2))
      if(icntyp.gt.0) in_contact=in_contact+1
      if (DEBUG) then
        write(65,990) ni2, lm(ni2), icntyp
      endif
  980 format(1x,'ni1=',i5,'  ln(ni1)=',i5,'  icntyp=', i5)
  990 format(1x,'ni2=',i5,'  ln(ni2)=',i5,'  icntyp=', i5)
c
c   End contact check
c------------------------------------------------------------
c
c   Assume pressure is on (no contact)
c
      factor =1.0
c
c   If both nodes of this side are in contact, turn off pressure
c
      if(in_contact.eq.2) then
         factor = 0.0
c         stays_off(n(1))=.true.
c         if (DEBUG) then
c           write(6,1001) inc,n(1),nn,factor,n(2),in_contact
c         endif
c
c     comment out special stuff for keeping pressure off
c      else 
c         if(stays_off(n(1))) then
c            factor = 0.0
c         if (DEBUG) then
c            write(6,1002) inc,n(1),nn,factor,n(2),in_contact
c         endif
c         else
c            factor = 1.0
c         if (DEBUG) then
c            write(6,1003) inc,n(1),nn,factor,n(2),in_contact
c         endif
c         end if 
      end if
 1000 format(1x,'*** End: el, intpt., n(2), in_contact, pr = ',4i5,2x,
     1       1pe10.3)
 1001 format(1x,'B1: inc,el,nn,factor,n(2),in_contact=',3i5,1pe10.3,2i5)
 1002 format(1x,'B2: inc,el,nn,factor,n(2),in_contact=',3i5,1pe10.3,2i5)
 1003 format(1x,'B3: inc,el,nn,factor,n(2),in_contact=',3i5,1pe10.3,2i5)
c
      press = factor * press
      if (DEBUG) then
        write(65,1000) n(1),nn,n(2),in_contact,press
      endif
      press_val(n(1)) = press
      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)
      include '../common/concom'
      common /mydata2/ppress_val
      pointer(ppress_val, press_val(1))
c
      v=press_val(m)
c      write(66,100) inc, m, v
c 100  format(1x,/,' Inc =', i5, '  El =', i5, '  Press =',1pe10.3,/)
      return
      end