subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,props,
     1 nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,kstep,kinc,
     2 jelem,params,ndload,jdltyp,adlmag,predef,npredf,lflags,
     3 mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)
c
      include 'aba_param.inc'
c
      dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),svars(nsvars),
     1 energy(8),props(*),coords(mcrd,nnode),
     2 u(ndofel),du(mlvarx,*),v(ndofel),a(ndofel),time(2),
     3 params(3),jdltyp(mdload,*),adlmag(mdload,*),
     4 ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),jprops(*)
      dimension sresid(6)

      parameter (zero=0.d0, one=1.d0)
      dimension x(9),y(11),vx(9,11),vy(9,11), coordsu(3), vel(3),
     1          accel(3)

      data x/0.0   ,  21.9 ,  43.8 ,  65.7 , 87.6  ,
     $       109.6 , 131.5 , 153.4 , 175.3  /
      data y/235.16, 232.81, 225.75, 214.0 , 197.53,
     $       176.37, 150.5 , 119.93,  84.66,  44.68, 0.0/
c
c     sresid - stores the static residual at time t+dt
c     svars  - contains the static residual at time t upon entering the
c              routine. sresid is copied to svars after the dynamic
c              residual has been calculated.
c
      area= props(1)
      e   = props(2)
      anu = props(3)
      rho = props(4)
c
      alen = coords(1,2) - coords(1,1)
      ak = area*e/alen
      am = 0.5d0*area*rho*alen
c
      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
      if (lflags(3).eq.1) then
c        normal incrementation
         if (lflags(1).eq.1 .or. lflags(1).eq.2) then
c           *static
            amatrx(1,1) =  ak  
            amatrx(4,4) =  ak  
            amatrx(1,4) = -ak  
            amatrx(4,1) = -ak
         else if (lflags(1).eq.11 .or. lflags(1).eq.12) then
c           *dynamic
            alpha = params(1)
            beta  = params(2)
            gamma = params(3)
c                  
            dadu = 1.0d0/(beta*dtime**2)
            dvdu = gamma/(beta*dtime)
c                  
            do k1=1,ndofel
              amatrx(k1,k1) = am*dadu
            end do
            amatrx(1,1) = amatrx(1,1) + (1.0d0+alpha)*ak  
            amatrx(4,4) = amatrx(4,4) + (1.0d0+alpha)*ak  
            amatrx(1,4) = amatrx(1,4) - (1.0d0+alpha)*ak  
            amatrx(4,1) = amatrx(4,1) - (1.0d0+alpha)*ak
         end if
      else if (lflags(3).eq.2) then
c        stiffness matrix
         amatrx(1,1) =  ak  
         amatrx(4,4) =  ak  
         amatrx(1,4) = -ak  
         amatrx(4,1) = -ak
      else if (lflags(3).eq.4) then
c        mass matrix
         do k1=1,ndofel
           amatrx(k1,k1) = am
         end do
      else if (lflags(3).eq.5) then
c        half-step residual calculation
      else if (lflags(3).eq.6) then
c        initial acceleration calculation
         do k1=1,ndofel
           amatrx(k1,k1) = am
         end do
      else if (lflags(3).eq.100) then
c        output for perturbations
      end if

      if (lflags(3) .gt. 1) return

      ndim = 2
      ij   = 0
      do j=1,11
        coordsu(2) = y(j)
        do i=1,9
          coordsu(1) = x(i)
          iwarn = 0
          call getwavevel (ndim, coordsu, vel, accel, iwarn, 
     1                     noel, exinterm)
          vx(i,j) = vel(1)
          vy(i,j) = vel(2)
          if (iwarn.eq.1) write(7,*) '****** iwarn .eq. 1 *******'

          ij = ij + 1
          svars(ij)     = vel(1)
          svars(ij+ 99) = vel(2)
          svars(ij+198) = accel(1)
          svars(ij+297) = accel(2)
        end do
      end do

      write(7,*)'velocity-vx'
      write(7,222)((vx(k1,k2),k1=1,9),k2=1,11)
      write(7,*)'velocity-vy'
      write(7,222)((vy(k1,k2),k1=1,9),k2=1,11)


 222  format(9f6.2)
 333  format(3f6.2)



      return
      end