subroutine elevar(nm,nnm,layer,gstran,gstres,stress,pstran,
     1cstran,vstran,cauchy,eplas,equivc,swell,krtyp,prang,dt,
     2gsv,ngens1,ngen1,nstats,nstass,thmstr)
c
c CONTACT: ap@marc.de
c      
c* * * * *
c
c User routine to compute the momentum and kinetic energy
c of two impacting bodies.
c It will print:
c   1) Total mass of both bodies.
c   2) The kinetic energy of both bodies.
c   3) The total strain energy of both bodies.
c   4) The total momentum of both bodies.
c   5) The average velocity of both bodies.
c   6) The total momentum of the model (momentum balance).
c
c NOTE: The element numbering must start at 1 and must
c       be consecutive. The first body must have the lowest
c       element numbers and the second body the highest.
c       The element type must either be type 7 (hexa 8)
c       or element type 75 (quad 4 shell)
c       or element type 9 (2 node truss).
c       The user has to define the parameter NELB1 (the highest
c       element number of body one) in the parameter definition
c       below.
c
c* * * * *
      implicit real*8 (a-h,o-z)                                               dp
      dimension gstran(ngens),gstres(ngens),
     1stress(ngen1),pstran(ngen1),cstran(ngen1),vstran(ngen1),
     2cauchy(ngen1),dt(nstats),gsv(1)
     3,thmstr(1),crang(3,3),krtyp(4)
c
      dimension veloc(3),ph(8),rvelo(12)
c* * * * *
c veloc = integration point velocity vector
c ph    = element shape function array
c rvelo = nodal point velocity vector
c* * * * *
      common/impuls_com/xm1,xm2,xpuls1(3),xpuls2(3),xpulst(3),
     *                  veloc1(3),veloc2(3),ekin1,ekin2,epot1,epot2
c* * * * *
c xm1    = total mass of body one.
c xm2    = total mass of body two.
c xpuls1 = total momentum of body one.
c xpuls2 = total momentum of body two.
c xpulst = total momentum of the model.
c veloc1 = average velocity of body one.
c veloc2 = average velocity of body two.
c ekin1  = kinetic energy of body one.
c ekin2  = kinetic energy of body two.
c epot1  = total strain energy of body one.
c epot2  = total strain energy of body two.
c* * * * *
      parameter(NELB1=120)
c* * * * *
c NELB1 = highest element number of body one.
c* * * * *
c
      include '/usr/software/marck62_R10000/common/strvar'
      include '/usr/software/marck62_R10000/common/strnen'
      include '/usr/software/marck62_R10000/common/arrays'
      include '/usr/software/marck62_R10000/common/array4'
      include '/usr/software/marck62_R10000/common/heat'
      include '/usr/software/marck62_R10000/common/dimen'
      include '/usr/software/marck62_R10000/common/elmcom'
      include '/usr/software/marck62_R10000/common/matdat'
      include '/usr/software/marck62_R10000/common/space'
      include '/usr/software/marck62_R10000/common/lass'
      include '/usr/software/marck62_R10000/common/blnk'
      include '/usr/software/marck62_R10000/common/creeps'
c
c... return if layer number not equal to one
      if (layer.gt.1) return
c
c... initialize the mass and momentum values
      if (nm.eq.1.and.nnm.eq.1.and.layer.eq.1) then
       xm1   = 0.0d0
       xm2   = 0.0d0
       ekin1 = 0.0d0
       ekin2 = 0.0d0
       epot1 = 0.0d0
       epot2 = 0.0d0
       do i1=1,3
        xpuls1(i1) = 0.0d0
        xpuls2(i1) = 0.0d0
        veloc1(i1) = 0.0d0
        veloc2(i1) = 0.0d0
       enddo
      endif
c
c... data base offset (here use internal element number)
      lofr = (n-1)*nelstr
c
c... get the integration point coordinates
c     la1  = icrxpt + lofr + (nn-1)*ncrd
c     do 100,i1=1,ncrd
c      xip(i1) = vars(la1+i1-1)
c100  continue
c
c... dv = the jacobian (the integration point volume factor)
      la1   = ijacoc  + lofr + nn - 1
      dv    = vars(la1)
      thick = 1.0d0
c
c... thick = the shell thickness (defaults to 1.0 for non-shells)
      if (jtype.eq.75) then
       la1   = ithick + lofr + nn -1
       thick = vars(la1)
      endif
c
c... scale jacobian with shell thickness
      dv = dv*thick
c
c... determine velocity at integration point
c ... determine isoparametric coordinates of integration point by
c ... binary decoding of the integration point number
      agaus = 1d0/dsqrt(3d0)
      nint  = nnm-1
      nze   = nint/4
      if (nze.eq.0) then
       cze = -agaus
      else
       cze =  agaus
      endif
      if (jtype.ne.7) cze = -1d0
      net = (nint-4*nze)/2
      if (net.eq.0) then
       cet = -agaus
      else
       cet =  agaus
      endif
      if (jtype.eq.9) cet = -1d0
      nxi = nint-4*nze-2*net
      if (nxi.eq.0) then
       cxi = -agaus
      else
       cxi =  agaus
      endif
      if (jtype.eq.9) cxi = 0d0
c ... compute integration point shape function values
      ph(1) = 0.125d0*(1d0-cxi)*(1-cet)*(1d0-cze)
      ph(2) = 0.125d0*(1d0+cxi)*(1-cet)*(1d0-cze)
      ph(3) = 0.125d0*(1d0+cxi)*(1+cet)*(1d0-cze)
      ph(4) = 0.125d0*(1d0-cxi)*(1+cet)*(1d0-cze)
      ph(5) = 0.125d0*(1d0-cxi)*(1-cet)*(1d0+cze)
      ph(6) = 0.125d0*(1d0+cxi)*(1-cet)*(1d0+cze)
      ph(7) = 0.125d0*(1d0+cxi)*(1+cet)*(1d0+cze)
      ph(8) = 0.125d0*(1d0-cxi)*(1+cet)*(1d0+cze)
c ... compute integration point velocity in veloc
      do i1=1,3
       veloc(i1) = 0d0
      enddo
      do i1=1,nnode
       jrdpre=0
       call vecftc(rvelo,vars(idynvs),ncrdmx,ncrd,lm(i1),jrdpre,2,1)
       do i2=1,3
        veloc(i2) = veloc(i2) + ph(i1)*rvelo(i2)
       enddo
      enddo
c
c... sum the integration point contributions
      rhodv = dabs(rho)*dv
      if (nm.le.NELB1) then
       xm1 = xm1 + rhodv
       do i1=1,3
        xpuls1(i1) = xpuls1(i1) + rhodv*veloc(i1) 
        ekin1      = ekin1 + 0.5d0*rhodv*veloc(i1)*veloc(i1)
       enddo
       do i1=1,ngen1
        epot1      = epot1 + 0.5d0*dv*stress(i1)*gstran(i1)
       enddo
      else
       xm2 = xm2 + rhodv
       do i1=1,3
        xpuls2(i1) = xpuls2(i1) + rhodv*veloc(i1) 
        ekin2      = ekin2 + 0.5d0*rhodv*veloc(i1)*veloc(i1)
       enddo
       do i1=1,ngen1
        epot2      = epot2 + 0.5d0*dv*stress(i1)*gstran(i1)
       enddo
      endif
c
c... print the sum when processing last element and integration point
      if (nm.eq.numel.and.nnm.eq.intel) then
       do i1=1,3
        xpulst(i1) = xpuls1(i1) + xpuls2(i1)
        veloc1(i1) = xpuls1(i1)/xm1
        veloc2(i1) = xpuls2(i1)/xm2
       enddo
       write(6,*)
c... print the total mass,kinetic and strain energy
       write(6,1000) xm1
       write(6,1010) xm2
       write(6,1500) ekin1
       write(6,1510) ekin2
       write(6,1600) epot1
       write(6,1610) epot2
c... print the momentum
       write(6,2000) (xpuls1(i1),i1=1,3)
       write(6,2010) (xpuls2(i1),i1=1,3)
       write(6,2020) (xpulst(i1),i1=1,3)
c... print the average velocities
       write(6,3000) (veloc1(i1),i1=1,3)
       write(6,3010) (veloc2(i1),i1=1,3)
c... write the data into special tables
       write(75,4000) cptim,(xpuls1(i1),i1=1,3),(veloc1(i1),i1=1,3)
       write(76,4000) cptim,(xpuls2(i1),i1=1,3),(veloc2(i1),i1=1,3)
       write(77,5000) cptim,(xpulst(i1),i1=1,3)
       write(78,6000) cptim,ekin1,ekin2,epot1,epot2
      endif
c
 1000 format(/,15x,'Total mass of body one : ',e14.6)
 1010 format(  15x,'Total mass of body two : ',e14.6)
 1500 format(/,15x,'Total kinetic energy of body one : ',e14.6)
 1510 format(  15x,'Total kinetic energy of body two : ',e14.6)
 1600 format(/,15x,'Total strain energy of body one : ',e14.6)
 1610 format(  15x,'Total strain energy of body two : ',e14.6)
 2000 format(/,15x,'Momentum of body one     : ', 3e14.6)
 2010 format(  15x,'Momentum of body two     : ', 3e14.6)
 2020 format(  15x,'Total momentum of system : ', 3e14.6)
 3000 format(/,15x,'Average velocity of body one : ', 3e14.6)
 3010 format(  15x,'Average velocity of body two : ', 3e14.6)
 4000 format(7e14.6)
 5000 format(4e14.6)
 6000 format(5e14.6)
c
      return
c
      end