c The following user-tailored subroutines are used in SPF for the following 
c reasons:
c 
c forcem.f .... To control pressure based on the material law and targer edot
c plotv.f ..... To write out the "calculated pressure" for post processing
c uedinc.f .... Stops the run after % nodes come in contact (user control)
c uinstr.f .... Pre stress the membrane for the 1st few increments
c urpflo.f .... Material law 
c stress1.f ... utility routine for urplow            
c 7/15/96 Reza Sadeghi
c=========================================
      subroutine forcem(ptot,x1,x2,nn,nnx)
c
      implicit real*8(a-h,o-z)
c
      include '../common/concom'
      include '../common/nzro1'
      include '../common/dimen'
      include '../common/array2'
      include '../common/arrays'
      include '../common/space'
      include '../common/heat'
      include '../common/creeps'
      include '../common/strtgy'
      include '../common/dyns'
      include '../common/far'
      include '../common/iautcr'
c
      common /spf/ fa,fm,pmin,pmax,target,pst,incres,porn,matl,
     +gstot,emin
      dimension nnx(2)
c
      if(nfff.eq.0) then
        if(inc.gt.0) incres = inc
        pnew    = pmin
        ptot    = porn*pmin
        lastinc = inc    
        nfff = 1
        go to 999
      endif
c
c----------------------------------------------------
c     return the same value for p if the same inc 
c----------------------------------------------------
c
      if(lastinc.eq.inc) then
         ptot = porn*pnew
         go to 999
      end if      
c
      lastinc = inc
      if(inc.eq.1) time0 = timinc
c
c-----------------------------------------------------------------------
c     calculate the average strain rate, the maximum incremental strain,
c     the maximum incremental displacement, and the maximum total
c     displacement
c     this works for out-of-core and for all element types
c
c     eavg = average strain rate
c     emax = maximum strain rate
c     dtot = maximum total displacement
c     dmax = maximum incremental displacement
c-----------------------------------------------------------------------
c
      eavg = 0.0d0
      emax = 0.0d0
      nx   = 1
c
      do 10 i1=1,numel
c
      mmx = i1
      if(ielsto.eq.0) nx = i1
      if(ielsto.eq.1) call wrat3(vars(ielsbn),nelsto,mmx,0,numel)
      lofrx = (nx - 1)*nelstr
      loffx = (nx - 1)*nelsto
c
      kitypx = ints(iitype + loffx)
      if(kitypx.eq.0) go to 10
      itypx  = iabs(kitypx)
c
      ishelx = ints(iishl + itypx - 1)
      ktpx = 1
      if(ishelx.eq.1) ktpx = nseqst
c
        j1 = lofrx + ieprat
        do 5 i2=1,nstres
        j2 = j1 + (i2 - 1)*neqst - 1
          do 3 i3=1,ktpx
          fac1 = vars(j2 + i3)
          etot = etot + fac1
          elavg = elavg + fac1
          if(emax.lt.fac1) emax = fac1
    3     continue 
    5   continue 
c**********************************************
c lines added to include cut-off strain rate...
c   elstrn = element strain-rate
c
        elstrn = elavg/dfloat(nstres*ktpx)
        if(elstrn.ge.emin)then
           actstrn = actstrn + elstrn
           kelact = kelact + 1
        end if
        elavg = 0.0d0
c**********************************************
   10 continue 
c
      if(ielsto.eq.1) call wrat3(vars(ielsbn),nelsto,nnx(1),0,numel)
c
      etot = etot/dfloat(nstres*ktpx*numel)
      if(emax.ge.emin)then
         eavg = actstrn/kelact
      else
         eavg = etot
      end if
c
      dmax = 0.0d0
      do 50 i1=1,numnp
      j1 = idsx  + (i1 - 1)*ndeg - 1
      j2 = idsxt + (i1 - 1)*ndeg - 1
        do 40  i2=1,ncrd
        fac1 = dabs(vars(j1 + i2))
        fac2 = dabs(vars(j2 + i2))
        if(dmax.lt.fac1) then
        dmax = fac1
        dtot = fac2
        endif
   40   continue
   50 continue
c
c----------------------------------------------------------------------
c     calculate a new value for pressure based on desired strain rate
c----------------------------------------------------------------------
c
      call stress1(eavg,starget,stress)
      sra    = starget/stress
      pnew   = sra*pnew
      if(pnew.ge.pmax) pnew = pmax
      if(pnew.le.pmin) pnew = pmin
      if(inc.eq.1    ) pnew = pmin
      ptot = porn*pnew
c
      if (matl .eq. 1)then
      write(0,1001) inc,pnew,target,starget,eavg,stress,emax,dmax,timinc
      write(6,1001) inc,pnew,target,starget,eavg,stress,emax,dmax,timinc
 
 1001 format('     increment number      = ',11x,i5,/,
     1       '     pressure              = ',e16.6,/,
     2       '     target strain rate    = ',e16.6,/,
     3       '     target stress         = ',e16.6,/,
     4       '     average strain rate   = ',e16.6,/,
     5       '     average stress        = ',e16.6,/,
     6       '     maximum strain        = ',e16.6,/,
     7       '     maximum displacement  = ',e16.6,/,
     8       '     time step             = ',e16.6)
      elseif (matl .eq. 2) then
      write(7,1002) inc,pnew,target,starget,eavg,savg,emax,dmax,
     1               timinc,gstot
      write(6,1002) inc,pnew,target,starget,eavg,savg,emax,dmax,
     1               timinc,gstot
c
 1002 format('     increment number      = ',11x,i5,/,
     1       '     pressure              = ',e16.6,/,
     2       '     target strain rate    = ',e16.6,/,
     3       '     target stress         = ',e16.6,/,
     4       '     average strain rate   = ',e16.6,/,
     5       '     average stress        = ',e16.6,/,
     6       '     maximum strain        = ',e16.6,/,
     7       '     maximum displacement  = ',e16.6,/,
     8       '     time step             = ',e16.6,/,
     9       '     grain size            = ',e16.6)
      end if
 999  continue
      return
      end
c================================================================
      subroutine stress1(eavg,starget,stress)
c     ****************************************************
c
c     Calculate stress based on desired material law.
c     Set "matl = 1" for the Power Law
c     Set "matl = 2" for Hamilton's Equation.
c     target = material target strain-rate
c     emin   = cut-off strain-rate(default=10% of target)
c     pmin   = minimum pressure (psi)
c     pmax   = maximum pressure (psi)
c     pst    = element pre-stress (psi)
c     porn   = pressure multiplier 
c              set to -1 for 3d
c              set to +1 for 2d
c
c     ****************************************************
      implicit real*8(a-h,o-z)                                               dp
      common /spf/ fa,fm,pmin,pmax,target,pst,incres,porn,matl,
     +gstot,emin
      include '../common/concom'
      include '../common/creeps'

      if(inc .le. 1)then 
c---------------------------------------
c SPF Paramater Constants...
      target = 0.0005
      emin   = 0.1*target
      pmin   = 0.001
      pmax   = 300.0
      pst    = 50.0
      porn   = -1
      matl   = 1
c---------------------------------------
c Power law material constants:
      fa     = 15300.0
      fm     = 0.5
c-------------------------------------------
c Constants for Hamilton's Equation:
c 
c   gs    = initial grain size (cm)
c   gsm   = incremental grain size (microns)
c   gstot = total grain size (microns)
c
c   *** 5083 Al Alloy ***
c
      gs    = 8.0e-4
      AII   = 8.0e-21
      p     = 4.3d0
      sm    = 0.6d0
      s0    = 10.0d0
      B     = 3.9e19
      q     = 22.0
      degg  = 0.31
      gstot = gs*10000.0
c
c   *** 6Al-4V Ti Alloy ***
c
c      AII = 1.3e-18
c      p  = 3.0d0
c      sm = 0.7d0
c      s0 = 25.0d0
c      B  = 0.0412
c      q  = 3.9
c      degg = 1300.0
c      gs = 6.0e-4
c      gstot = gs*10000.0
c---------------------------------------------
      end if

      if (matl .eq. 1)then
         starget = fa*target**fm
         stress = fa*eavg**fm
         return

      elseif (matl .eq. 2)then
c
c        calculate grain growth (microns) based on the avg strain rate 
c        at the beginning of each increment...
c
         gsm = ((B/(q*gstot**(q-1.0))) + degg*gstot*eavg)*timinc
         if(linc.ne.inc) then
            gstot = gsm + gstot
            gs = gstot * 0.0001
            linc = inc
         endif
c
c        calculate equivalent stress
c
         starget = (((target*(gs**p)/AII)**sm) + s0)*1.45
         stress  = (((eavg*(gs**p)/AII)**sm) + s0)*1.45
         return
      endif
      end
c================================================================
      subroutine uinstr(s,ndi,nshear,m,nn,kc,xintp,ncrd,inc,cptim,
     1                  timinc)
      implicit real*8(a-h,o-z)                                               dp
      common /spf/ fa,fm,pmin,pmax,target,pst,incres,porn,matl,
     +gstot,emin
      dimension s(1),xintp(1)
c
      if(inc.le.5) then
      s(1) = pst
      s(2) = pst
      s(3) = 00.0d0
      endif
c
      return
      end
c
c ******************************************************************************
c
      subroutine urpflo(mdum,nn,layer,mats,inc,ndi,ngens,ncrd,nstat,
     1                  cptim,timinc,ebar,erate,dt,dtdl,stats,dstats,
     2                  coord,yd)
c     
c**************************************************************************
c      user routine for rigid-plastic flow:
c      -----------------------------------
c
c
c      passed into routine:
c      -------------------
c      mdum   = element number
c      nn     = integration point number
c      layer  = layer number
c      mats   = material set number
c      inc    = increment number
c      ndi    = number of direct components
c      ngens  = total number of components
c      nstat  = number of state variables excluding temperature
c               (nstat = one less than number of state variables
c               defined on "state var,X" command)
c      cptim  = time at beginning of increment
c      timinc = incremental time
c      dt     = temperature at beginning of increment
c      dtdl   = incremental temperature
c      ebar   = total equivalent strain at beginning of increment
c      stats  = values of state variables excluding temperature
c               at beginning of increment
c      erate  = equivalent strain rate
c      stats  = array of state variables (excluding temperature)
c               at beginning of increment
c      coord  = integration point coordinates
c
c      to be passed back:
c      -----------------
c
c      yd     = equivalent stress; if not calculated here, program will
c               find the value of yd from the input data
c      dstats = incremental state variables (excluding temperature)
c
c**************************************************************************
c
      implicit real*8(a-h,o-z)
      common /spf/ fa,fm,pmin,pmax,target,pst,incres,porn,matl,
     +gstot,emin
      dimension mdum(2),stats(nstat),dstats(nstat),coord(ncrd)
      call stress1(erate,starget,stress)
      yd = stress
      return
      end
c=================================================================
      subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,nndi,
     * nnshear,jpltcd)
c* * * * * *
c
c     select a variable contour plotting (user subroutine).
c
c     v            variable
c     s (idss)         stress array
c     sp           stresses in preferred direction
c     etot          total strain (generalized)
c     eplas         total plastic strain
c     ecreep        total creep strain
c     t             current temperature
c     m            element number
c     nn           integration point number
c     layer        layer number
c     ndi (3)       number of direct stress components
c     nshear (3)    number of shear stress components
c
c* * * * * *
      implicit real*8 (a-h,o-z) 
      dimension s(1),etot(1),eplas(1),ecreep(1),sp(1)            
      dimension ndint(5000),cnode(4),dnode(4),coor(5000,4)
      dimension et(5000),icon(5000,4),ipair(5000)    
      common /spf/ fa,fm,pmin,pmax,target,pst,incres,porn,matl,
     +gstot,emin
      include '../common/elmcom'
      include '../common/concom'
      include '../common/dimen'
      include '../common/arrays'
      include '../common/array2'
      include '../common/develp'
      include '../common/space'
      include '../common/prepro'
      include '../common/blnk'
      include '../common/again'
     
      goto(100,200) jpltcd

100   v=pnew 
      return

200   if(lclass.ne.5)  then 
        v = 0.
        return
      end if 
      if(inc.eq.0.or.inc.eq.incres) then
        if(lastm.ne.m) then
          icon(m,1) = lm(1)
          icon(m,2) = lm(2)
          icon(m,3) = lm(3)
          icon(m,4) = lm(4)
        end if
        if(m.eq.numel) then
          do i = 1,numel
            ipair(i) = 0
            do j = 1,numel
              j1 = icon(i,4)
              j2 = icon(i,3)
              if(icon(j,1).eq.j1.and.icon(j,2).eq.j2) ipair(i) = j
            end do
          end do
        end if
      end if
c
      if(inc.ne.lastinc) then
        do i = 1,numnp
          ndint(i) = ibsrch(i,int(inoids),numnp,1)
          la2 = inpnum+ndint(i)-1
          if(joptit.ne.0) ndint(i) =igetsh(ints(la2),0)
c
          call vecftc(cnode,vars(ixord),ncrdmx,ncrd,ndint(i),0,2,1)
          call vecftc(dnode,vars(idsxt),ndegmx,ndeg,ndint(i),0,2,5)
            do j = 1,ncrd
              coor(ndint(i),j) = cnode(j) + dnode(j)
            end do
        end do          
        do i=1,numel
          if(ipair(i).ne.0) then
            i1=icon(i,1)
            i2=icon(i,2)
            i3=icon(ipair(i),3)
            i4=icon(ipair(i),4)
            t1=dsqrt((coor(i4,1)-coor(i1,1))**2+
     +               (coor(i4,2)-coor(i1,2))**2)
            t2=dsqrt((coor(i3,1)-coor(i2,1))**2+
     +               (coor(i3,2)-coor(i2,2))**2)
            et(i)=(t1+t2)/2.0
            et(ipair(i))=et(i)      
          end if
        end do
      end if
      lastinc = inc
      lastm = m
      v=et(m) 
      return
      end
c====================================================================
      subroutine uedinc(idum1,idum2)
      implicit real*8 (a-h,o-z)                                               dp
c* * * * * *
c     Stop when a certain fraction of all nodes are in contact.
c* * * * * *
      include '../common/form'
      include '../common/array4'
      include '../common/concom'
      include '../common/space'
      include '../common/blnk'
      include '../common/dimen'
      include '../common/machin'
      dimension nod(4)
      equivalence (nod(1),lm(1))
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c At 99% of all nodes in contact, The run is terminated.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      done=.99
      ndef=numdie-1
      kount=0
      do 100 i1=1,ndef
      nbct = ints(inbct+i1-1)
      j1 = itouch + (i1 - 1)*nbcn
      j2 = inf    + (i1 - 1)*nbcn
c     write(6,300) numdie,nbct,nbcn,j1,j2
300   format(' numdie,nbct,nbcn,j1,j2=',5i10)
        do 90 i2=1,nbct
        idie = ints(j1+i2-1)
        inod = ints(j2+i2-1)
        if(idie.eq.0) go to 90
        kount=kount+1
c          write(6,200) i1,idie,i2,inod,nod(k)
200        format(' i1,idie,i2,inod,nod(k)=',5i10)
   90   continue
  100 continue
      fins=float(kount)/float(numnp)
      if(kount.gt.klast) write(iprtl,400) kount,numnp
400   format(31x,i5,' nodes of ',i5,' in contact')
      klast=kount
      if(fins.ge.done) then
        don100=100*done
        write(6,500) don100
        write(iprtl,500) don100
500     format(39x,'marc exit number 3004',/, 
     +  31x,f5.1,' percent of all nodes in contact',/,
     +  33x,'program stop in subroutine uedinc')
c23456789012345678901234567890123456789012345678901234567890
c                                      marc exit number 3004
c                              xxxxx percent of all nodes in contact
c                                program stop in spf user subroutines
        stop
      end if
      return
      end
c================================================================