C
C     USER SUBROUTINE TO CONTROL INCREMENTATION FOR FREE ROLLING
      subroutine urdfil(lstop,lovrwrt,kstep,kinc,pnewdt,time)
c     
c     
      include 'aba_param.inc'
c     
      dimension array(513),jrray(nprecd,513),time(*),lout(2)
      equivalence (array(1),jrray(1,1))
      parameter (nmany = 999999,zero=0.0d0,small=1.0d-6,cmax=2.d-2,
     *     cmax1=8.d-1,one=1.0d0)
      parameter (omega=10.d0,ftol=10.d0)
      common /cfroll1/ omega1, omega0, temp1, temp0
      common /cfroll2/ kflag
      common /cfroll3/ time0, time1
C
C     ftol  - the tolerance for zero force
C     cmax  - max increment in rolling speed allowed
C     omega - upper bound estimate for the free-rolling velocity 
C     kflag - flag set when free-rolling soln. is obtained
C
c      write(7,*) "Entering urdfil", kstep, kinc, pnewdt, kflag
      if (kflag .eq. 1) return 
      if (kinc .eq. 0) return
C     
C     SPECIFY NODE NUMBER AND DEGREE OF FREEDOM TO MONITOR
C     
      ldofusr = 5
      nodeusr = 199
      otol = cmax*omega
      write(7,1000)
C     
C     INITIALLY INCREMENT SOLUTION WITH SMALL STEPS
C     
      if (kstep.lt.3 .or. kinc .eq. 1) then
         pnewdt = cmax
         otol1  = cmax1*omega  
         time0  = time(1)
         time1  = time(1)+cmax1
         omega0 = otol1
         omega1 = omega0+otol
         return
      end if
C     
C     FIND PREVIOUS INCREMENT
      do k1=kinc-1,1,-1
         call dbfile(2,array,jrcd)
         call posfil(kstep,k1,array,jrcd)
         if (jrcd .eq. 0) goto 10
      end do
C     EXIT IF RESULTS FOR PREVIOUS INCREMENT CANNOT BE FOUND
      call xit
      return
 10   continue
      lproc = jrray(1,7)
      if (lproc .ne. 85) then
C     EXIT IF THE PROCEDURE IS NOT SST
         call xit
         return
      end if
      time0 = array(4)
      kinc0 = jrray(1,9)
      omega0 = temp0
C     
C     FIND THE TORQUE AT PREVIOUS INCREMENT
      torq0 = zero
      lfnd0 = 0
      do k1=1,nmany
         call dbfile(0,array,jrcd)
         if (jrcd .ne. 0) goto 100
         key = jrray(1,2)
         if (key .eq. 104) then
            node = jrray(1,3)
            if (node .eq. nodeusr) then
               torq0 = array(3+ldofusr)
               lfnd0 = 1
               goto 100
            end if
         else if (key .eq. 2001) then
            goto 100
         end if
      end do
C     EXIT IF END-OF-FILE MARKER NOT FOUND
      call xit
      return
 100  continue
      if (lfnd0 .eq. 0) then
         lout(1) = nodeusr
         lout(2) = ldofusr
C     EXIT IF CANNOT FIND FORCES FOR NODE
         call xit
         return
      end if      
C     
C     FIND CURRENT INCREMENT
      call posfil(kstep,kinc,array,jrcd)
      if (jrcd .ne. 0) then
C     EXIT IF RESULTS FOR CURRENT INCREMENT CANNOT BE FOUND
         call xit
         return
      end if      
      time1 = array(4)
      kinc1 = jrray(1,9)
      omega1 = temp1
C
C     FIND THE TORQUE AT THIS INCREMENT
      torq1 = zero
      lfnd1 = 0
      do k1=1,nmany
         call dbfile(0,array,jrcd)
         if (jrcd .ne. 0) goto 200
         key = jrray(1,2)
         if (key .eq. 104) then
            node = jrray(1,3)
            if (node .eq. nodeusr) then
               torq1 = array(3+ldofusr)
               lfnd1 = 1
               goto 200
            end if
         end if
      end do
C     EXIT IF END-OF-FILE MARKER NOT FOUND
      call xit
      return
 200  continue
      if (lfnd1 .eq. 0) then
        lout(1) = nodeusr
        lout(2) = ldofusr
C     EXIT IF CANNOT FIND FORCES FOR NODE
        call xit
        return
      end if      
C
C     ESTIMATE ANGULAR VELOCITY TO OBTAIN FREE ROLLING
      denom = torq1-torq0
      write(7,1001) time0,torq0,time1,torq1
      if (abs(torq1) .le. ftol) then
         write (7,1004)
         kflag = 1
         pnewdt= one - time(1)
         return
      else if (abs(denom) .lt. small*abs(torq1)) then
         write (7,1002) kinc1, kinc0
         domega = otol
      else
         omegafr = (torq1*omega0-torq0*omega1)/denom
         write (7,1003) omegafr
         domega = omegafr-omega1
         if (abs(domega) .gt. abs(otol)) then
            write (7,1005) otol
            domega = otol
         end if
      end if
      omega0 = omega1
      omega1 = omega1+domega
      pnewdt = cmax
      time0 = time(1)
      time1 = time(1)+pnewdt
      write(7,1006) time0, omega1, omega0
C     
 1000 format(//,15x,'FREE ROLLING SOLUTION CONTROL')
 1001 format(/,2x,'TORQUE AT TIME ',1pg12.5,' = ',1pg12.5,
     *     /,2x,'TORQUE AT TIME ',1pg12.5,' = ',1pg12.5)
 1002 format(/,2x,'TORQUE AT CURRENT INCREMENT (',i5,
     *     ') AND PREVIOUS INCREMENT (',i5,') ARE THE SAME.',
     *     /,2x,'CONTINUE WITH REGULAR INCREMENTATION.',/)
 1003 format(2x,'FREE ROLLING ESTIMATED TO OCCUR AT SPEED ',1pg12.5)
 1004 format(/,2x,'FREE ROLLING SOLUTION OBTAINED.',
     *     ' END STEP.',//)
 1005 format(2x,'ESTIMATED INCREMENT TO OBTAIN FREE ROLLING',
     *     ' EXCEEDS ',1pg12.5,
     *     /,2x,'CONTINUE WITH REGULAR INCREMENTATION.',/)
 1006 format(/2x,'END CALCULATIONS AT TIME',1pg12.5,
     *     'NEW ANGULAR VELOCITY ',1pg12.5,
     *     'PREVIOUS ESTIMATE ',1pg12.5,/)
C     
      return
      end
C     
C     User subroutine to define motion for free-rolling
      subroutine umotion(u,kstep,kinc,time,node,jdof)
c     
      include 'aba_param.inc'
c     
      dimension time(2)
      parameter (nmany = 999999,zero=0.0d0,small=1.0d-6,cmax=2.d-2,
     *     cmax1=8.d-1,one=1.0d0)
      parameter (omega=10.d0,ftol=10.d0)
      common /cfroll1/ omega1, omega0, temp1, temp0
      common /cfroll2/ kflag
      common /cfroll3/ time0, time1
C     
      if (kinc .eq. 0) return
      if (kinc .eq. 1) then
         kflag = 0
         u = cmax1*omega
         temp0=omega
         temp1=u
         return
      end if
C     
C     CHECK IF FREE ROLLING SOLUTION IS OBTAINED
C     
      if (kflag .eq. 1) then
         u = omega1
         return
      else
C     
C     CHECK FOR CUTBACK OR INCREASE IN INCREMENT
C     
         if (time(1) .eq. time1) then     
            u = omega1
         else
            u = omega0+(omega1-omega0)/(time1-time0)*(time(1)-time0)
         end if
         temp0=omega0
         temp1=u
      end if
C     
      return
      end