************************************************************************
*
*  subroutine to calculate force acting between segments  ver.1 2003.8
*
************************************************************************
      subroutine bodyforce(dfbx,dfby,dfbz,fxup,fyup,fzup
     &             ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &             ,x0a,y0a,z0a,x0b,y0b,z0b,xb,yb,zb,j
     &             ,xg,yg,zg)
      implicit double precision(a-z)
      integer j,jf,nf,numb,nm
      
      parameter(numb =  21)
      parameter(nf   =  21)
      parameter(nm   =  21)

      double precision  fxup(nf),fxdw(nf)
     &                 ,fyup(nf),fydw(nf)
     &                 ,fzup(nf),fzdw(nf)
     &                 ,mxup(nm),mxdw(nm)
     &                 ,myup(nm),mydw(nm)
     &                 ,mzup(nm),mzdw(nm)
     &       ,x0a(numb),y0a(numb),z0a(numb)
     &       ,x0b(numb),y0b(numb),z0b(numb) 
     &                 ,dfbx,dfby,dfbz
     &                 ,u(3),v(3)

      jf=1
      if((j.ge.8).and.(j.le.15))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=2
      if(((j.ge.8).and.(j.le.15)).or.(j.eq.1))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=3
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.2)))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=4
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.3)))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=5
      if((j.eq.5).or.(j.eq.6).or.(j.eq.7))then
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      else
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      endif

      jf=6
      if((j.eq.6).or.(j.eq.7))then
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      else
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      endif

      jf=7
      if(j.eq.7)then
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      else
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      endif

      jf=8
      if((j.ge.8).and.(j.le.15))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=9
      if((j.ge.9).and.(j.le.15))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=10
      if((j.eq.10).or.(j.eq.12).or.(j.eq.14))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=11
      if((j.eq.11).or.(j.eq.13).or.(j.eq.15))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=12
      if((j.eq.12).or.(j.eq.14))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=13
      if((j.eq.13).or.(j.eq.15))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=14
      if(j.eq.14)then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=15
      if(j.eq.15)then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=16
      if((j.eq.16).or.(j.eq.18).or.(j.eq.20))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=17
      if((j.eq.17).or.(j.eq.19).or.(j.eq.21))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=18
      if((j.eq.18).or.(j.eq.20))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=19
      if((j.eq.19).or.(j.eq.21))then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=20
      if(j.eq.20)then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif

      jf=21
      if(j.eq.21)then
         fxdw(jf) = fxdw(jf)+dfbx
         fydw(jf) = fydw(jf)+dfby
         fzdw(jf) = fzdw(jf)+dfbz
      else
         fxup(jf) = fxup(jf)+dfbx
         fyup(jf) = fyup(jf)+dfby
         fzup(jf) = fzup(jf)+dfbz
      endif


******************
* bending moment *
******************

      u(1) = (xb - xg)
      u(2) = (yb - yg)
      u(3) = (zb - zg)

      v(1) = dfbx
      v(2) = dfby
      v(3) = dfbz

      jf=1
      if((j.ge.8).and.(j.le.15))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=2
      if(((j.ge.8).and.(j.le.15)).or.(j.eq.1))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=3
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.2)))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=4
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.3)))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=5
      if((j.eq.5).or.(j.eq.6).or.(j.eq.7))then
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      else
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      endif

      jf=6
      if((j.eq.6).or.(j.eq.7))then
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      else
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      endif

      jf=7
      if(j.eq.7)then
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      else
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      endif

      jf=8
      if((j.ge.8).and.(j.le.15))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=9
      if((j.ge.9).and.(j.le.15))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=10
      if((j.eq.10).or.(j.eq.12).or.(j.eq.14))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=11
      if((j.eq.11).or.(j.eq.13).or.(j.eq.15))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=12
      if((j.eq.12).or.(j.eq.14))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=13
      if((j.eq.13).or.(j.eq.15))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=14
      if(j.eq.14)then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=15
      if(j.eq.15)then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=16
      if((j.eq.16).or.(j.eq.18).or.(j.eq.20))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=17
      if((j.eq.17).or.(j.eq.19).or.(j.eq.21))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=18
      if((j.eq.18).or.(j.eq.20))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=19
      if((j.eq.19).or.(j.eq.21))then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=20
      if(j.eq.20)then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif

      jf=21
      if(j.eq.21)then
         mxdw(jf) = mxdw(jf)+exprox(u,v)
         mydw(jf) = mydw(jf)+exproy(u,v)
         mzdw(jf) = mzdw(jf)+exproz(u,v)
      else
         mxup(jf) = mxup(jf)+exprox(u,v)
         myup(jf) = myup(jf)+exproy(u,v)
         mzup(jf) = mzup(jf)+exproz(u,v)
      endif


      return
      end

*************************************************************
*
*  subroutine to calculate inertial torque    ver.1 2003.9
*
*************************************************************
      subroutine inertorque(e01x,e01y,e01z
     &                     ,e02x,e02y,e02z
     &                     ,e03x,e03y,e03z
     &                     ,di1,di2,di3,dm
     &                     ,awxg,awyg,awzg
     &                     ,ome1,ome2,ome3
     &                     ,mxup,myup,mzup,mxdw,mydw,mzdw
     &                     ,x0a,y0a,z0a,x0b,y0b,z0b,xb,yb,zb,i,j
     &                     ,xg,yg,zg,inerh,dt)

      implicit double precision(a-z)
      integer i,k,i1,i2,j,k1,k2,jf,numb,nm,div
      
      parameter(numb =  21)
      parameter(nm   =  21)
      parameter(div  =  50)

      double precision e01x,e01y,e01z
     &                ,e02x,e02y,e02z
     &                ,e03x,e03y,e03z
     &                ,mxup(nm),mxdw(nm)
     &                ,myup(nm),mydw(nm)
     &                ,mzup(nm),mzdw(nm)
     &                ,x0a(numb),y0a(numb),z0a(numb)
     &                ,x0b(numb),y0b(numb),z0b(numb) 
     &                ,di1,di2,di3,ome1,ome2,ome3
     &                ,awxg,awyg,awzg,dm
     &                ,u(3),v(3),iner(3,3),inerd(3,3)
     &                ,inerd2(3,3),p(3,3),ome(3),aome(3)
     &                ,inerh(numb,div,nm,3,3),diner(3,3)

      do 10 i1=1,3
      do 10 i2=1,3
         iner (i1,i2) = 0.0d0
         inerd(i1,i2) = 0.0d0
  10  continue

* tensor of principal moment of inertia
      iner(1,1) = di1
      iner(2,2) = di2
      iner(3,3) = di3

* rotate tensor into direction of absolute coordinate
      p(1,1) = e01x
      p(2,1) = e01y
      p(3,1) = e01z
      p(1,2) = e02x
      p(2,2) = e02y
      p(3,2) = e02z
      p(1,3) = e03x
      p(2,3) = e03y
      p(3,3) = e03z

      do 11 i1=1,3
      do 11 i2=1,3
      do 11 k1=1,3
      do 11 k2=1,3
         inerd(i1,i2) = inerd(i1,i2) + p(i1,k1)*iner(k1,k2)*p(i2,k2)
  11  continue

      ome(1)=ome1
      ome(2)=ome2
      ome(3)=ome3

      aome(1)=awxg
      aome(2)=awyg
      aome(3)=awzg

* calculate bending moment

         u(1) = (xb - xg)
         u(2) = (yb - yg)
         u(3) = (zb - zg)

         inerd2(1,1) = inerd(1,1) + (u(2)*u(2)+u(3)*u(3))*dm
         inerd2(2,1) = inerd(2,1) -  u(1)*u(2)*dm
         inerd2(3,1) = inerd(3,1) -  u(1)*u(3)*dm
         inerd2(1,2) = inerd(1,2) -  u(1)*u(2)*dm
         inerd2(2,2) = inerd(2,2) + (u(1)*u(1)+u(3)*u(3))*dm
         inerd2(3,2) = inerd(3,2) -  u(2)*u(3)*dm
         inerd2(1,3) = inerd(1,3) -  u(1)*u(3)*dm
         inerd2(2,3) = inerd(2,3) -  u(2)*u(3)*dm
         inerd2(3,3) = inerd(3,3) + (u(1)*u(1)+u(2)*u(2))*dm


      do jf=1,nm

* finite difference of tensor and save it
      do i1 = 1,3
      do i2 = 1,3
         diner(i1,i2) = (inerd2(i1,i2) - inerh(j,i,jf,i1,i2))/dt
         inerh(j,i,jf,i1,i2) = inerd2(i1,i2)
      enddo
      enddo

      do k=1,3
         v(k)=0.0d0
      enddo

* calculate inertia torque
      do i1=1,3
      do k=1,3
         v(i1) = v(i1) + inerd2(i1,k)*ome(k)
      enddo
      enddo

      if(jf.eq.1)then
      if((j.ge.8).and.(j.le.15))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.2)then
      if(((j.ge.8).and.(j.le.15)).or.(j.eq.1))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.3)then
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.2)))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.4)then
      if(((j.ge.8).and.(j.le.15)).or.((j.ge.1).and.(j.le.3)))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.5)then
      if((j.eq.5).or.(j.eq.6).or.(j.eq.7))then
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.6)then
      if((j.eq.6).or.(j.eq.7))then
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.7)then
      if(j.eq.7)then
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.8)then
      if((j.ge.8).and.(j.le.15))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.9)then
      if((j.ge.9).and.(j.le.15))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.10)then
      if((j.eq.10).or.(j.eq.12).or.(j.eq.14))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.11)then
      if((j.eq.11).or.(j.eq.13).or.(j.eq.15))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.12)then
      if((j.eq.12).or.(j.eq.14))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.13)then
      if((j.eq.13).or.(j.eq.15))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.14)then
      if(j.eq.14)then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.15)then
      if(j.eq.15)then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.16)then
      if((j.eq.16).or.(j.eq.18).or.(j.eq.20))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.17)then
      if((j.eq.17).or.(j.eq.19).or.(j.eq.21))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.18)then
      if((j.eq.18).or.(j.eq.20))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.19)then
      if((j.eq.19).or.(j.eq.21))then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.20)then
      if(j.eq.20)then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else if(jf.eq.21)then
      if(j.eq.21)then
         mxdw(jf) = mxdw(jf) - exprox(ome,v)
         mydw(jf) = mydw(jf) - exproy(ome,v)
         mzdw(jf) = mzdw(jf) - exproz(ome,v)
         do i1 =1,3
         mxdw(jf) = mxdw(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         mydw(jf) = mydw(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzdw(jf) = mzdw(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      else
         mxup(jf) = mxup(jf) - exprox(ome,v)
         myup(jf) = myup(jf) - exproy(ome,v)
         mzup(jf) = mzup(jf) - exproz(ome,v)
         do i1 =1,3
         mxup(jf) = mxup(jf)-inerd2(1,i1)*aome(i1)-diner(1,i1)*ome(i1)
         myup(jf) = myup(jf)-inerd2(2,i1)*aome(i1)-diner(2,i1)*ome(i1)
         mzup(jf) = mzup(jf)-inerd2(3,i1)*aome(i1)-diner(3,i1)*ome(i1)
         enddo
      endif

      else
      endif

      enddo

      return
      end


***********************************************************
*
* subroutine to calculate torque component due to distance 
* between calculation surface and mass center 
* and force acting on the surface
*
*                               ver.1.0 2003.11
***********************************************************
      subroutine rexfa(fxup,fyup,fzup,fxdw,fydw,fzdw
     &     ,mxup,myup,mzup,mxdw,mydw,mzdw
     &     ,x0a,y0a,z0a,x0b,y0b,z0b,xg,yg,zg)

      implicit double precision(a-z)
      integer jf,nf,numb,nm
      
      parameter(numb =  21)
      parameter(nf   =  21)
      parameter(nm   =  21)

      double precision  fxup(nf),fxdw(nf)
     &                 ,fyup(nf),fydw(nf)
     &                 ,fzup(nf),fzdw(nf)
     &                 ,mxup(nm),mxdw(nm)
     &                 ,myup(nm),mydw(nm)
     &                 ,mzup(nm),mzdw(nm)
     &       ,x0a(numb),y0a(numb),z0a(numb)
     &       ,x0b(numb),y0b(numb),z0b(numb)
     &       ,r(3),f(3)

      do 100 jf =1,numb
         r(1) = x0a(jf) - xg
         r(2) = y0a(jf) - yg
         r(3) = z0a(jf) - zg

         f(1) = fxup(jf)
         f(2) = fyup(jf)
         f(3) = fzup(jf)

         mxup(jf) = mxup(jf) - exprox(r,f)
         myup(jf) = myup(jf) - exproy(r,f)
         mzup(jf) = mzup(jf) - exproz(r,f)

         f(1) = fxdw(jf)
         f(2) = fydw(jf)
         f(3) = fzdw(jf)

         mxdw(jf) = mxdw(jf) - exprox(r,f)
         mydw(jf) = mydw(jf) - exproy(r,f)
         mzdw(jf) = mzdw(jf) - exproz(r,f)
 100  continue

      return
      end

***********************************************************
*
*  subroutine to coordinate transform force and moment
*
*                               ver.1.0 2003.11
***********************************************************
      subroutine rfandm(fxup,fyup,fzup,fxdw,fydw,fzdw
     &                 ,mxup,myup,mzup,mxdw,mydw,mzdw
     &     ,fxbup,fybup,fzbup,fxbdw,fybdw,fzbdw
     &     ,mxbup,mybup,mzbup,mxbdw,mybdw,mzbdw
     &     ,e01x,e01y,e01z,e02x,e02y,e02z,e03x,e03y,e03z)

      implicit double precision(a-z)
      integer jf,nf,nm
      
      parameter(numb =  21)
      parameter(nf   =  21)
      parameter(nm   =  21)

      double precision  fxup(nf),fxdw(nf)
     &                 ,fyup(nf),fydw(nf)
     &                 ,fzup(nf),fzdw(nf)
     &                 ,mxup(nm),mxdw(nm)
     &                 ,myup(nm),mydw(nm)
     &                 ,mzup(nm),mzdw(nm)
     &                 ,fxbup(nf),fxbdw(nf)
     &                 ,fybup(nf),fybdw(nf)
     &                 ,fzbup(nf),fzbdw(nf)
     &                 ,mxbup(nm),mxbdw(nm)
     &                 ,mybup(nm),mybdw(nm)
     &                 ,mzbup(nm),mzbdw(nm)

      do 421 jf=1,nf
         fxd = fxup(jf) *e01x + fyup(jf) *e01y + fzup(jf) *e01z
         fyd = fxup(jf) *e02x + fyup(jf) *e02y + fzup(jf) *e02z
         fzd = fxup(jf) *e03x + fyup(jf) *e03y + fzup(jf) *e03z

         fxbup(jf) = fxd 
         fybup(jf) = fyd 
         fzbup(jf) = fzd 

         fxd = fxdw(jf) *e01x + fydw(jf) *e01y + fzdw(jf) *e01z
         fyd = fxdw(jf) *e02x + fydw(jf) *e02y + fzdw(jf) *e02z
         fzd = fxdw(jf) *e03x + fydw(jf) *e03y + fzdw(jf) *e03z

         fxbdw(jf) = fxd 
         fybdw(jf) = fyd 
         fzbdw(jf) = fzd 
 421  continue

      do 422 jf=1,nm
         mxd = mxup(jf) *e01x + myup(jf) *e01y + mzup(jf) *e01z
         myd = mxup(jf) *e02x + myup(jf) *e02y + mzup(jf) *e02z
         mzd = mxup(jf) *e03x + myup(jf) *e03y + mzup(jf) *e03z

         mxbup(jf) = mxd 
         mybup(jf) = myd 
         mzbup(jf) = mzd 

         mxd = mxdw(jf) *e01x + mydw(jf) *e01y + mzdw(jf) *e01z
         myd = mxdw(jf) *e02x + mydw(jf) *e02y + mzdw(jf) *e02z
         mzd = mxdw(jf) *e03x + mydw(jf) *e03y + mzdw(jf) *e03z

         mxbdw(jf) = mxd 
         mybdw(jf) = myd 
         mzbdw(jf) = mzd 
 422  continue

      return
      end
