** this is a simplified CALSPAN tire model to be used
** with the 1994 Chevy C15000 pick-up truck model (tr_*.inp)
** it includes only radial forces in the tires 
** lateral forces/moments and self-aligning moments have been stripped off
** and are not included in this version of the UEL
*******************************************************************
      subroutine crossprod(a,b,c)

      implicit real*8(a-h,o-z)
      
      dimension a(3),b(3),c(3)

      c(1) = a(2)*b(3) - a(3)*b(2)
      c(2) = a(3)*b(1) - a(1)*b(3)
      c(3) = a(1)*b(2) - a(2)*b(1)

      return
      end
*******************************************************************
      subroutine dotprod(a,b,c)

      implicit real*8(a-h,o-z)
      
      dimension a(3),b(3)

      c = a(1)*b(1) + a(2)*b(2) + a(3)*b(3)

      return
      end
*******************************************************************
c     matrix multiplication
      subroutine matrixmult(a,b,c,n1,n2,n3)

      implicit real*8(a-h,o-z)
      
      parameter(zero=0.d0)

      dimension a(n1,n2),b(n2,n3),c(n1,n3)

      do i=1,n1
         do j = 1,n3
            c(i,j) = zero
            do k = 1,n2
               c(i,j) = c(i,j) + a(i,k)*b(k,j)
            end do
         end do
      end do

      return
      end

c     compute outer product of two vectors
      subroutine outprod(a,b,aoutb)
      
      implicit real*8(a-h,o-z)

      dimension a(*), b(*), aoutb(3,*)
      
      do i = 1,3
         do j = 1,3
            aoutb(i,j) = a(i)*b(j)
         end do 
      end do

      return 
      end 
      
***************************************************
      subroutine rotatevector(bign,an,rot,nvec,
     $     large)
C
c     input: bign -> vectors to be rotated
c          : rot  -> rotation vector
c          : nvec -> number of vectors to be rotated
c          : large-> 1 if NLGEOM=yes, 0 otherwise
c    output: an   -> rotated vectors

      implicit real*8(a-h,o-z)
      parameter(precis=1.d-16,bigang=15848925.0d0,one=1.d0)

      dimension bign(3,*),an(3,*),rot(3)
C
      dimension dir(3)

      if(large.eq.0)then
         do kvec=1,nvec
            an(1,kvec)=bign(1,kvec)
            an(2,kvec)=bign(2,kvec)
            an(3,kvec)=bign(3,kvec)
         end do
      else if(rot(1)**2+rot(2)**2+rot(3)**2.le.precis)then
         do kvec=1,nvec
            an(1,kvec)=bign(1,kvec)+bign(3,kvec)*rot(2)
     *           - bign(2,kvec)*rot(3)
            an(2,kvec)=bign(2,kvec)+bign(1,kvec)*rot(3)
     *           - bign(3,kvec)*rot(1)
            an(3,kvec)=bign(3,kvec)+bign(2,kvec)*rot(1)
     *           - bign(1,kvec)*rot(2)
         end do
      else
         angle=min(sqrt(rot(1)**2+rot(2)**2+rot(3)**2),bigang)
         ca=cos(angle)
         sa=sin(angle)
         oa=one/angle
         dir(1)=oa*rot(1)
         dir(2)=oa*rot(2)
         dir(3)=oa*rot(3)
         do kvec=1,nvec
            t2=(dir(1)*bign(1,kvec)+dir(2)*bign(2,kvec)
     1           +dir(3)*bign(3,kvec))*(one-ca)
            an(1,kvec)=t2*dir(1)+ca*bign(1,kvec)+sa*(dir(2)*bign(3,kvec)
     1           -dir(3)*bign(2,kvec))
            an(2,kvec)=t2*dir(2)+ca*bign(2,kvec)+sa*(dir(3)*bign(1,kvec)
     1           -dir(1)*bign(3,kvec))
            an(3,kvec)=t2*dir(3)+ca*bign(3,kvec)+sa*(dir(1)*bign(2,kvec)
     1           -dir(2)*bign(1,kvec))
         end do
      end if

      return
      end
***************************************************
c     load stiffness matrix
      subroutine loadstiff(amatrx,ndofel,uu,tt,ut,tu,tt_self,tu_self,
     $           stiffrad,llateral,lselfalign,oalphad,ldebug)
 
      implicit real*8(a-h,o-z)

      dimension uu(3,3),ut(3,3),tu(3,3),tt(3,3),tt_self(3,3),
     $     tu_self(3,3),amatrx(ndofel,ndofel),stiffrad(3,3)

      do i = 1,3
         do j = 1,3
            aux = oalphad*stiffrad(i,j)
            amatrx(i,j)     =  aux
         end do
      end do
      
      if (ldebug.ge.2)then 
         write(7,105)((amatrx(i,j),i=1,6),j=1,6)
 105     format(6('am: ', 6(1E13.6,1x),/))
      end if

      return 
      end
*******************************************************************
C
C  User subroutine for a simplified CALSPAN tire

      subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,
     1     props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,
     2     kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,
     3     lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)
c
* UnusedArgument nsvars stevenso
* UnusedArgument nprops stevenso
* UnusedArgument du stevenso
* UnusedArgument a stevenso
* UnusedArgument jtype stevenso
* UnusedArgument ndload stevenso
* UnusedArgument jdltyp stevenso
* UnusedArgument adlmag stevenso
* UnusedArgument predef stevenso
* UnusedArgument npredf stevenso
* UnusedArgument ddlmag stevenso
* UnusedArgument mdload stevenso
* UnusedArgument pnewdt stevenso
* UnusedArgument njprop fox
* UnusedVariable limpul oancea
C
      implicit real*8(a-h,o-z)
c
      dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),svars(*),props(*),
     1   energy(7),coords(mcrd,nnode),u(ndofel),du(mlvarx,*),v(ndofel),
     2   a(ndofel),time(2),params(*),jdltyp(mdload,*),adlmag(mdload,*),
     4   ddlmag(mdload,*),predef(2,npredf,nnode),lflags(4),jprops(*)
c
**      character*255 tname1,tname2
      character*256 outdir
      character*14 fname_loc,tname1
      character*22 fname_static,tname2

      parameter(zero=0.d0,half=0.5d0,one=1.d0,asmall=1.d-36,abig=1.d36,
     $     pi=3.14159265358979d0,two=2.d0,p25=0.25d0,p1=0.1d0,
     $     asmal=1.d-15)

      parameter(ldebug = 0)

***************************************************************
c
c     parameters that may  be changed for different car models
c
c     distances between cg and front/rear wheels in the x-direction
      parameter(x_cg_init=2252.728d0,
     *          centerfront=1387.532d0,
     *          centerback=1977.2087d0)

c     original orientations of the wheel axis in nominal config
      parameter(xw_axis=0.d0,yw_axis=-1.d0,zw_axis=0.d0)

c     tire element numbers (left = driver side)
      parameter(
     $     lfrontright=701000,lfrontleft=702000,
     $     lbackright=704000,lbackleft=703000)

c     write output to separate files (not thread safe!)
      parameter(lwrite=0)

***************************************************************  
c    common block for uel output only
      common /uelout/time_out,x_cg,
     $ fr_lf,fl_lf,tmz_lf,a_lf,g_lf,b_lf,xb_lf,zb_lf,c_lf,
     $ fr_rf,fl_rf,tmz_rf,a_rf,g_rf,b_rf,xb_rf,zb_rf,c_rf,
     $ fr_rb,fl_rb,tmz_rb,a_rb,g_rb,b_rb,xb_rb,zb_rb,c_rb,
     $ fr_lb,fl_lb,tmz_lb,a_lb,g_lb,b_lb,xb_lb,zb_lb,c_lb,
     $ kinc_old,lstatic_force_written

***************************************************************

      dimension curcoor(6,1),curvec(6),origvec(6),
     $     sresid(100),
     $     vbar(3,3),pv(3,3),e2orig(3),e2(3), eL(3),dnor(3),
     $     ev(3),evp(3),
     $     frad(3),flat(3),flatmom(3),falign(3),stiffrad(3,3),
c
     $     v1(3),v2(3),v9(3),v10(3),v11(3),v12(3),v13(3),v14(3),
     $     v15(3),v16(3),v17(3),et(3),vofb(3),v20(3),v21(3),
c
     $     eLv13(3,3),eLv14(3,3),dnorv15(3,3),
     $     v16v13(3,3),v17et(3,3),v16v14(3,3),
     $     dnorv20(3,3),dnorv21(3,3),vofbet(3,3),
c
     $     pva(3,3),pn(3,3),v100(3),v101(3),pv100(3,3),pv101(3,3),
     $     pv100pn(3,3),aa(3,3),bb(3,3),bbf(3),bb1(3),bb2(3),bb3(3),
     $     cc1(3),cc2(3),cc3(3),cc(3,3),dd(3,3),
c
     $     uu(3,3),ut(3,3),tu(3,3),tt(3,3),tt_self(3,3),tu_self(3,3)
      
      logical lopened

c     set solution flags
      large = 0 
      if (lflags(2).eq.1) large = 1
      limpul = 0
      if (lflags(3).eq.6) limpul = 3
      ldyn = 0
      if (lflags(1).eq.11 .or. lflags(1).eq.12) ldyn = 1
      lstatic = 0
      if (lflags(1).eq.1 .or. lflags(1).eq.2) lstatic = 1

c     assign the original direction of tire axis (nominal config)
      e2orig(1) = xw_axis
      e2orig(2) = yw_axis
      e2orig(3) = zw_axis

c     direction of travel
      et(1) = one
      et(2) = zero
      et(3) = zero

c     read properties
      stiff  = props(1)
      reflen = props(2)

      bumplength = props(3) !5
      bumpheigth = props(4) !0.2
      x_bumpstart= props(5) !arbitrary
      linphase   = jprops(1) !bumps in phase or not
      llateral   = jprops(2) !1-include lateral force effects
      lselfalign = jprops(3) !1-include self aligning torque

c     test if a tire is on the right or on the left side
      llefttire  = 0
      if (jelem.eq.lbackleft .or. jelem.eq.lfrontleft)
     $     llefttire = 1

c     zero stiffness and RHS
      do k1 = 1,ndofel
         sresid(k1) = zero
         do krhs = 1,nrhs
            rhs(k1,krhs) = zero
         end do
         do k2 = 1,ndofel
            amatrx(k2,k1) = zero
         end do
      end do

c     compute current nodal coordinates
      ndim = mcrd
      if (ndim.gt.3) ndim=3
      do i = 1,ndim
         do j=1,nnode
            curcoor(i,j) = coords(i,j)+u((j-1)*6+i)
         end do
      end do

c     corrdinates of the wheel center
      xw = curcoor(1,1)
      yw = curcoor(2,1)
      zw = curcoor(3,1)
      if (ldebug.ne.0) then 
         write(7,*)'*************** BEGIN *******************'
         write(7,*)'element number = ', jelem
         write(7,2005)xw,yw,zw
 2005    format('xw= ',1E13.6, ' yw= ',1E13.6, ' zw= ', 1E13.6)
      end if

c     compute current bump parameters: xb, zb, beta angle     
      xb = xw
      zb = zero
      aux = zero
      aofxa = zero
      if (linphase .eq. -1) then  
c        steering tests: no bumps
      else
         x_bumpend   = x_bumpstart + bumplength
         if (xb.gt.x_bumpstart .and. xb.lt.x_bumpend) then 
            if (ldebug.ne.0) write(7,*)'doing the bump'
            p4 = two/bumplength
            trigarg = p4*pi*(xb - x_bumpstart)
            zb  = (half*bumpheigth)*(one-cos(trigarg))
            aux = (half*bumpheigth)*p4*pi*sin(trigarg)
            aofxa = p4*pi*aux
            if (linphase.eq.0) then 
               if (llefttire.eq.1) then 
                  zb  = - zb
                  aux = -aux 
               end if
            end if
         end if
      end if
      beta  = atan(aux)
      cbeta = cos(beta)
      sbeta = sin(beta)
      ocbeta = one/cbeta

c     current coordinates of the contact point Q
      xq = xw + (zw-zb)*sbeta*cbeta
      yq = yw
      zq = zw*sbeta*sbeta + zb*cbeta*cbeta

c     distance from wheel center to contact point
      distQ = sqrt((xq-xw)**2 + (yq-yw)**2 + (zq-zw)**2)
      if (ldebug.ne.0) write(7,*)'distq = ',distQ

c     compute current orientation of the wheel axis
      call rotatevector(e2orig,e2,u(4),1,large)

c     normal to the contact plane
      dnor(1) =-sbeta   
      dnor(2) = zero
      dnor(3) = cbeta

c     projection of the wheel axis onto ground plane
      e2dotdnor = e2(1)*dnor(1) + e2(2)*dnor(2) +  e2(3)*dnor(3)
      eL(1) = e2(1) - e2dotdnor*dnor(1)
      eL(2) = e2(2) - e2dotdnor*dnor(2)
      eL(3) = e2(3) - e2dotdnor*dnor(3)

c     make the projected vector a unit vector         
c     thus eL is a unit vector contained in the ground plane
c     perpendicular to the wheel plane
      eLmag  = sqrt(eL(1)**2 + eL(2)**2 +  eL(3)**2)
      if (eLmag.gt.zero) then
         eL(1) = eL(1)/eLmag
         eL(2) = eL(2)/eLmag
         eL(3) = eL(3)/eLmag
      else
         write(7,*) 'ERROR: wrong data in uel.f'
         goto 990
      end if

      gamma = zero
      alpha = zero
      lmove = 0
      if (ldebug.ne.0) then 
         write(7,11)(evp(i),i=1,3)
 11      format('evp = ',3(1x,1E13.6))
         write(7,12)(e2(i),i=1,3)
 12      format('e2 = ',3(1x,1E13.6))
         write(7,13)(el(i),i=1,3)
 13      format('eL = ',3(1x,1E13.6))
         write(7,14)(dnor(i),i=1,3)
 14      format('dnor = ',3(1x,1E13.6))
         write(7,*)'gamma = ', gamma, '  alpha = ', alpha
      end if

c     radial force on the tire (simplified Calspan tire) 
      rtu = reflen
      fr_entire = zero
      lcontact = 0              ! assume wheel out of contact with the ground
      disthubground = zw - rtu      
      if (zb .gt. disthubground) lcontact = 1
      if (lcontact.eq.1) fr_entire = stiff*(zb-zw+rtu)*ocbeta
      if (ldebug.ne.0) then 
         if (lcontact.eq.0) then 
            write(7,*)'contact lossed'
         else
            write(7,*)'contact gained'
         end if
      end if
      fr = fr_entire

c     other forces on the tire
      fl  = zero
      tmz = zero
      gfunc = zero
      gprime = zero
      c1 = zero
      c2 = zero
      xmuy = zero
      af1 = zero
      af2 = zero
      af3 = zero
   
**      call calspan(fl,tmz, fr,alpha,gamma,gfunc,gprime,c1,c2,xmuy,
**     $     af1,af2,af3,ldebug)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c            put forces and moments in vector form         
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      do ii = 1,3
         frad(ii) = zero
         flat(ii) = zero
         flatmom(ii) = zero
         falign(ii) = zero
      end do

c     radial force
      frad(1) = fr*dnor(1)    
      frad(2) = fr*dnor(2)
      frad(3) = fr*dnor(3)                  

c     lateral force
      flat(1) = -fl*eL(1)       
      flat(2) = -fl*eL(2)
      flat(3) = -fl*eL(3)         
      
c     moment due to the lateral force  = -dnor x flat
      call crossprod(flat,dnor,flatmom)
      do ii = 1,3
         flatmom(ii) = flatmom(ii)*distQ
      end do

c     self-aligning torque
      falign(1)  = dnor(1)*tmz
      falign(2)  = dnor(2)*tmz
      falign(3)  = dnor(3)*tmz

      if (ldebug.ne.0) then 
         write(7,1999)(u(i),i=1,6)
         write(7,2000)frad(1),frad(2),frad(3)
         write(7,2001)flat(1),flat(2),flat(3)
         write(7,2002)flatmom(1),flatmom(2),flatmom(3)
         write(7,2003)falign(1),falign(2),falign(3)
 1999    format('disp-rot = ', 6(1x,1e13.6))
 2000    format('frad   = ',3(1x,1E13.6))
 2001    format('flat   = ',3(1x,1E13.6))
 2002    format('flatmom= ',3(1x,1E13.6))
 2003    format('falign = ',3(1x,1E13.6))
      end if

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                  stiffness from radial force
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do i = 1,3
         do j = 1,3
            vbar(i,j) = dnor(i)*dnor(j)
            pv (i,j)  = zero
            if (distQ.gt.asmall) then 
               pv(i,j) = -vbar(i,j)
               if (i.eq.j) pv(i,j) = one+pv(i,j)
               pv(i,j) = pv(i,j)/distQ
            end if
            stiffrad(i,j) = stiff*vbar(i,j) + fr*pv(i,j)
         end do 
      end do

      if (ldebug.ge.2) then 
         write(7,100)((stiffrad(i,j),i=1,3),j=1,3)
 100     format(3('sr: ', 3(1E13.6,1x),/))
      end if

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                     assemble residual
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     if not moving forward, do not add contributions from lateral forces
c     and torques
      if (lmove.eq.0) then 
         llateral = 0
         lselfalign = 0
      end if

      do ii = 1,3
         if (llateral.eq.0) then 
            flat(ii) = zero
            flatmom(ii) =zero
         end if
         if (lselfalign.eq.0) falign(ii) = zero
      end do

      do i =1,3
         sresid(i)   = -(frad(i) + flat(i))
         sresid(i+3) = -(falign(i) + flatmom(i))
      end do

      if (ldebug.ne.0) then 
         write(7,130)(sresid(i),i=1,6)
 130     format('sresid = ',6(1x,1E13.6))
      end if

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c              element kinematics - for energy only
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      relpos = zero
      relcon = zero
      constref = zero
      alphad  = zero
      oalphad = zero

c     compute unit vector between nodes and relative position
      do i = 1,ndim
         curvec(i) =  - curcoor(i,1)
         origvec(i)=  - coords(i,1)
      end do
      if (ndim.eq.2) curvec(3) = zero


      dist = zero
      dlen = zero
      do i = 1,ndim
         dlen = dlen +  curvec(i)*curvec(i)
         dist = dist +  origvec(i)*origvec(i)
      end do
      dlen = sqrt(dlen)
      dist = sqrt(dist)

      if (dlen.le.asmall) then  !pick a direction 
         relpos = zero
      else
         relpos = dlen
      end if

c     compute reference length if not specified
      constref= reflen
      if (reflen.eq.-abig) constref = dist
      relcon = relpos - constref

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                analysis type dependent coding
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      nr_entries = 6
      if (lflags(3).eq.1) then 
c        normal incrementation
         if (lstatic.ne.0) then !*STATIC

            oalphad = one
            call loadstiff(amatrx,ndofel,uu,tt,ut,tu,tt_self,tu_self,
     $           stiffrad,llateral,lselfalign,oalphad,ldebug)

c           RHS
            do i=1,ndofel
               rhs(i,1)    = rhs(i,1) - sresid(i)
            end do
            energy(2) = half*fr*relcon
         else if (ldyn.ne.0) then !*DYNAMIC
            alphad = params(1)
            oalphad = one+alphad

            call loadstiff(amatrx,ndofel,uu,tt,ut,tu,tt_self,tu_self,
     $           stiffrad,llateral,lselfalign,oalphad,ldebug)

            do i=1,nr_entries
               rhs(i,1)  = rhs(i,1)-(oalphad*sresid(i)-alphad*svars(i))
            end do

            if (ldebug.ne.0) then 
               write(7,131)(rhs(i,1),i=1,nr_entries)
 131           format('rhs = ',6(1x,1E13.6))
            end if
            
C              SVARS  - 1-nr_entries contains the static residual at time t 
C              upon entering the routine. SRESID is copied to 
C              SVARS(1-nr_entries) after the dynamic residual has been 
C              calculated.
            do k1 = 1,ndofel
               svars(k1+nr_entries) = svars(k1)
               svars(k1)   = sresid(k1)
            end do
            
            energy(1) = zero    !no kinetic energy
            energy(2) = half*fr*relcon
         end if
      else if (lflags(3).eq.2) then 
c        stiffness calculation
         oalphad = one
         call loadstiff(amatrx,ndofel,uu,tt,ut,tu,tt_self,tu_self,
     $        stiffrad,llateral,lselfalign,oalphad,ldebug)
      else if (lflags(3).eq.4) then 
c           mass calculation
c           there is no mass associated with this element         
      else if (lflags(3).eq.5) then 
c           half-step residual calculation
         alphad = params(1)
C           For half-step residual calculations: 
c           nr_entries+1...2*nr_entries
C           contains the static residual at the beginning 
C           of the previous increment. SVARS(1...nr_entries) are copied 
C           into SVARS(nr_entries+1...2*nr_entries) 
c           after the dynamic residual has been calculated.
         do i=1,nr_entries
            rhs(i,1) = rhs(i,1)-oalphad*sresid(i)+
     $           half*alphad*(svars(i)+svars(i+nr_entries))
         end do         
      else if (lflags(3).eq.6) then 
c        initial acceleration calculation
         do i=1,nr_entries
            rhs(i,1)    = rhs(i,1)   - sresid(i)
         end do         
         do k1 = 1,ndofel
            svars(k1) = sresid(k1)
         end do
         energy(1) = zero
         energy(2) = half*fr*relcon
      end if

c     write output 

c     initiliaze common block vars for output and open needed files
      if (kinc.eq.1 .and. kstep.eq.1)  then 
         x_cg   = x_cg_init
         lstatic_force_written = 0

c        get the current output directory
         loutdr = 0
         outdir = ' '
         call getoutdir(outdir, loutdr)

c        assign file name
         fname_loc  = 'truck-uel.data'

c        create full path name
         tname1 = fname_loc
**         tname1= outdir(1:loutdr)//fname_loc

         
         if (ldebug.ge.2) then 
            write(7,*)'tname1 = ', tname1
            write(7,*)'fname_loc = ',fname_loc
**            write(7,*)'outdir = ', outdir(1:loutdr)
         end if

c        open file (only once)
         if (lwrite.eq.1)then 
            inquire(file=tname1, opened=lopened)
            if (.not.lopened) then 
c              if not opened, try to open an  OLD file
               iold = 0
               open(unit=150,file=tname1,status='OLD',iostat=iold)
               if (iold.eq.0) then 
                  rewind(150)
               else
c                  write(7,*)'open 150 new'               
                  open(unit=150,file=tname1,status='NEW')
               end if
            end if

c           open file to write the static forces in the tires at the end of
c           the 1st static step
            fname_static = 'tire_static_force.data'
            tname2= fname_static
            if (ldebug.ge.2) then 
               write(7,*)'tname2 = ', tname2
            end if
            
**            tname2= outdir(1:loutdr)//fname_static
            inquire(file=tname2, opened=lopened)
            if (.not.lopened) then 
c              if not opened, try to open an  OLD file
               iold = 0
c               write(7,*)'open 151'
               open(unit=151,file=tname2,status='OLD',iostat=iold)
               if (iold.eq.0) then 
                  rewind(151)
               else
                  open(unit=151,file=tname2,status='NEW')
               end if
            end if
         end if

      end if

      if (kstep.eq.1 .and. lstatic_force_written.eq.0 .and.
     *     lwrite.eq.1) then 
         write(151,11111)fr_lf,fr_rf,fr_rb,fr_lb
11111    format('RADIAL FORCES IN TIRES AT STATIC EQUILIBRIUM:',/
     $          '   FRONT-LEFT  TIRE: ',1E13.6,/
     $          '   FRONT-RIGHT TIRE: ',1E13.6,/
     $          '   REAR-RIGHT  TIRE: ',1E13.6,/
     $          '   REAR-LEFT   TIRE: ',1E13.6)
         lstatic_force_written = 1
         close(151)
      end if
      
      lstart = 0
      if ((x_cg+centerfront).ge.(x_bumpstart-p1)) lstart = 1

      if (kinc.ne.kinc_old .and. lstart.eq.1 .and. lwrite.eq.1) then

         write(150,10002)time_out,x_cg,
     $        fr_lf,
     $        fr_rf,
     $        fr_lb,
     $        fr_rb
10002    format(6(1x,1e13.6))

      end if

      kinc_old = kinc

      if (ldebug.ge.2) then 
         write(7,*)'time12 = ', time(1),time(2)
         write(7,*)' x-Bstart= ', x_bumpstart
      end if
      time_out = time(2)      

 990  continue

      return
      end

*******************************************************************