************************************************************************
*
*   Swumsuit Analysis Engine Main Program
*
*                                          since 2002.7.19
*                      (c) Motomu Nakashima, Yasufumi Miura, 2002-2006
*
************************************************************************
      subroutine main(multi,numb
     &           ,numl,numl1,numl2,cycle,draggli,cn,cs,ca,lastloop
     &           ,div,divr,divt,inipos,inidirec,inivel,iniavel,direcd
     &           ,flagoutr,flagoutf,flagoutj,flagouta,numsp
     &           ,comsolv,conam
     &           ,vgxsolv,vgysolv,vgzsolv,ome1solv,ome2solv,ome3solv
     &           ,flagpene,flagoutm,pcn)
      implicit double precision(a-z)
      integer i,j,j2,j3,j4,j5,k,kk,div,divr,divt,ite,step,numb,loop
     &       ,numl,numl1,numl2,conam,multi,numbd
      integer ii,ij,imin,imax,imax1,imax2,in
      integer kd,rkloop,lastloop
      integer numsp,nbu,nf,nm,jf,ii1,ii2,igomi
      integer jj,ja,jl,div0
      integer flagoutr,flagoutf,flagoutj,flagouta,flagoutm
      integer comsolv,vgxsolv,vgysolv,vgzsolv,ome1solv,ome2solv,ome3solv
      integer js,je,jp,jjs,jje,jjp,flagpene(0:10),pene(0:10)
      integer goorstop
****** number of joint angle for Anybody**************************
      integer dummy,jcount, dum1,dum2,dum3
*******************************************************************
      character*10 chars

*     numb: number of segment
      parameter(numbd= 21)
      parameter(nbu =  9)
      parameter(nf  = 21)
      parameter(nm  = 21)

      double precision u(3),v(3)
      double precision ra1(numbd),ra2(numbd),rb1(numbd),rb2(numbd)
     &   ,len(numbd),leng(numbd),m(numbd)
     &   ,i1(numbd),i2(numbd),i3(numbd)  
     &   ,x1a(numbd),y1a(numbd),z1a(numbd)
     &   ,x1b(numbd),y1b(numbd),z1b(numbd)  
     &   ,ten(3,3)
     &   ,tent(3,3)
     &   ,iner(3,3),inerall(3,3)
     &   ,e11x(numbd),e11y(numbd),e11z(numbd)
     &   ,e12x(numbd),e12y(numbd),e12z(numbd)
     &   ,e13x(numbd),e13y(numbd),e13z(numbd)
     &   ,x1g(numbd),y1g(numbd),z1g(numbd)
     &   ,mata(3,3),matb(3,3),matc(3,3)
     &   ,e11xall(0:numl2),e11yall(0:numl2),e11zall(0:numl2)
     &   ,iall1(0:numl2)
     &   ,e12xall(0:numl2),e12yall(0:numl2),e12zall(0:numl2)
     &   ,iall2(0:numl2)
     &   ,e13xall(0:numl2),e13yall(0:numl2),e13zall(0:numl2)
     &   ,iall3(0:numl2)
     &   ,idall1(numl1)
     &   ,idall2(numl1)
     &   ,idall3(numl1)
     &   ,x1gall(0:numl2),y1gall(0:numl2),z1gall(0:numl2)
     &   ,xia(numbd,0:numl2),yia(numbd,0:numl2),zia(numbd,0:numl2)
     &   ,xib(numbd,0:numl2),yib(numbd,0:numl2),zib(numbd,0:numl2)  
     &   ,ei1x(numbd,0:numl2),ei1y(numbd,0:numl2),ei1z(numbd,0:numl2)
     &   ,ei2x(numbd,0:numl2),ei2y(numbd,0:numl2),ei2z(numbd,0:numl2)
     &   ,ei3x(numbd,0:numl2),ei3y(numbd,0:numl2),ei3z(numbd,0:numl2)
     &   ,vxia(numbd,numl1),vyia(numbd,numl1),vzia(numbd,numl1)
     &   ,vxib(numbd,numl1),vyib(numbd,numl1),vzib(numbd,numl1)  
     &   ,axia(numbd,numl1),ayia(numbd,numl1),azia(numbd,numl1)
     &   ,axib(numbd,numl1),ayib(numbd,numl1),azib(numbd,numl1)  
     &   ,x0a(numbd),y0a(numbd),z0a(numbd)
     &   ,x0b(numbd),y0b(numbd),z0b(numbd)  
     &   ,x0ac(numbd,0:numl2),y0ac(numbd,0:numl2),z0ac(numbd,0:numl2)
     &   ,x0bc(numbd,0:numl2),y0bc(numbd,0:numl2),z0bc(numbd,0:numl2)  
     &   ,x0ad(numbd),y0ad(numbd),z0ad(numbd)
     &   ,x0bd(numbd),y0bd(numbd),z0bd(numbd)  
     &   ,e01x(numbd),e01y(numbd),e01z(numbd)
     &   ,e02x(numbd),e02y(numbd),e02z(numbd)
     &   ,e03x(numbd),e03y(numbd),e03z(numbd)
     &   ,e01xc(numbd,0:numl2),e01yc(numbd,0:numl2),e01zc(numbd,0:numl2)
     &   ,e02xc(numbd,0:numl2),e02yc(numbd,0:numl2),e02zc(numbd,0:numl2)
     &   ,vx0a(numbd),vy0a(numbd),vz0a(numbd)
     &   ,vx0b(numbd),vy0b(numbd),vz0b(numbd)  
     &   ,ax0a(numbd),ay0a(numbd),az0a(numbd)
     &   ,ax0b(numbd),ay0b(numbd),az0b(numbd)  
     &   ,fbx(numbd),fby(numbd),fbz(numbd)
     &   ,fsx(numbd),fsy(numbd),fsz(numbd)
     &   ,ftx(numbd),fty(numbd),ftz(numbd)
     &   ,fnx(numbd),fny(numbd),fnz(numbd)
     &   ,fix(numbd),fiy(numbd),fiz(numbd)
     &   ,ptx(numbd),pty(numbd),ptz(numbd)
     &   ,pnx(numbd),pny(numbd),pnz(numbd)
     &   ,pix(numbd),piy(numbd),piz(numbd)
     &   ,mb1(numbd),mb2(numbd),mb3(numbd)
     &   ,ms1(numbd),ms2(numbd),ms3(numbd)
     &   ,mt1(numbd),mt2(numbd),mt3(numbd)
     &   ,mn1(numbd),mn2(numbd),mn3(numbd)
     &   ,mi1(numbd),mi2(numbd),mi3(numbd)
     &   ,inera(6,6),ineraorg(6,6),fv(6),fvorg(6),ipvt(6),zmat(6)
     &   ,inerad(6,6,numbd,div)
     &   ,fvd(6,4)
     &   ,rho(numbd)
     &   ,tsp(numsp),usp(numsp)
     &   ,xdist(numbd,div)
     &   ,ydist(numbd,div)
     &   ,zdist(numbd,div)
     &   ,xdistc(numbd,div,0:numl2)
     &   ,ydistc(numbd,div,0:numl2)
     &   ,zdistc(numbd,div,0:numl2)
     &   ,fxdist(numbd,div)
     &   ,fydist(numbd,div)
     &   ,fzdist(numbd,div)
     &   ,fxdistc(numbd,div,0:numl2)
     &   ,fydistc(numbd,div,0:numl2)
     &   ,fzdistc(numbd,div,0:numl2)
     &   ,fbxdist(numbd,div)
     &   ,fbydist(numbd,div)
     &   ,fbzdist(numbd,div)
     &   ,fsxdist(numbd,div)
     &   ,fsydist(numbd,div)
     &   ,fszdist(numbd,div)
     &   ,ftxdist(numbd,div)
     &   ,ftydist(numbd,div)
     &   ,ftzdist(numbd,div)
     &   ,fnxdist(numbd,div)
     &   ,fnydist(numbd,div)
     &   ,fnzdist(numbd,div)
     &   ,fixdist(numbd,div)
     &   ,fiydist(numbd,div)
     &   ,fizdist(numbd,div)
     &   ,bu(nbu)
     &   ,fxup(nf),fxdw(nf)
     &   ,fyup(nf),fydw(nf)
     &   ,fzup(nf),fzdw(nf)
     &   ,fxbup(nf),fxbdw(nf)
     &   ,fybup(nf),fybdw(nf)
     &   ,fzbup(nf),fzbdw(nf)
     &   ,mxup(nm),mxdw(nm)
     &   ,myup(nm),mydw(nm)
     &   ,mzup(nm),mzdw(nm)
     &   ,mxbup(nm),mxbdw(nm)
     &   ,mybup(nm),mybdw(nm)
     &   ,mzbup(nm),mzbdw(nm)
     &   ,axg,ayg,azg,awxg,awyg,awzg
     &   ,axifa(numbd),ayifa(numbd),azifa(numbd)
     &   ,axifb(numbd),ayifb(numbd),azifb(numbd)
     &   ,numen2(numbd,50)
     &   ,inerh(numbd,50,nm,3,3),gomi(nm)
     &   ,inipos(3),inidirec(3,3),inivel(3),iniavel(3)
     &   ,cn(numbd),cs(numbd),ca(numbd)
     &   ,vgxg(numl),vgyg(numl),vgzg(numl)
     &   ,ome1g(numl),ome2g(numl),ome3g(numl)
     &   ,xbm(numbd,div,divt),ybm(numbd,div,divt),zbm(numbd,div,divt)
     &   ,dismax(numbd)
c Below line has been modified from 4.2.0
     &   ,omev1(numl1),omev2(numl1),omev3(numl1)
     &   ,xgorgd(0:numl2),ygorgd(0:numl2),zgorgd(0:numl2)
     &   ,xcob(0:numl2),ycob(0:numl2),zcob(0:numl2)
     &   ,e1xorgd(0:numl2),e1yorgd(0:numl2),e1zorgd(0:numl2)
     &   ,e2xorgd(0:numl2),e2yorgd(0:numl2),e2zorgd(0:numl2)
     &   ,e3xorgd(0:numl2),e3yorgd(0:numl2),e3zorgd(0:numl2)
*********body movement,gravity,back direction for Anybody************
     &   ,rtan1,rtan2,rtan3
     &   ,role,pitch,yaw
     &   ,role0,pitch0,yaw0
     &   ,xori,yori
     &   ,waistx,waisty,waistz
     &   ,wstx0,wsty0,wstz0
     &   ,gravX,bd,jbd
*********fluid force,time step for Anybody***************************
     &   ,tany
     &   ,fxany(numbd,10)
     &   ,fyany(numbd,10)
     &   ,fzany(numbd,10)
*********************************************************************
     &   ,waitdum(100,100)
      double precision pcn

      real ineralls(3,3)
      integer lda,n,ldv,job,info
      complex evalc(3),evecc(3,3)
      real work(6)
      character lno*3

      pi=3.14159265358979d0

      dum1 = 0
      dum2 = 0
      dum3 = 0
      pcn =dble(int(pcn*100000.0d0))/100000.0d0
c     write(*,*) 'pcn=',pcn

* non-dimensionalize gravitational acceleration
* grav = 9.8d0 * (cycle)**2 / (human height)
* note that not non-dimensionalize human height here.
      grav =9.8d0 * cycle **2 / 1.0d0
      ttt  = cycle

* coefficient for fluid force due to momentum change at water surface
* however, this force is not considered (so set 0).
      cm = 0.0d0

* adjust direction of propulsion
* (for non-simmetrical motion)
      direc = - direcd * pi/180.0d0

*****************
* set time step *
*****************
      dt = 1.0/dble(numl)


**********************************************
**********************************************
* set velocity and angular velocity if given *
**********************************************
**********************************************
********
* vgxg *
********
      if (vgxsolv.eq.0) then
          open(unit=23,file='input_vgx.dat',status='unknown')
          do i=1,numl,2
              read(23,*)  vgxg(i)
              vgxg(i+1) = vgxg(i)
          enddo
      end if
      close(unit=23)

********
* vgyg *
********
      if (vgysolv.eq.0) then
          open(unit=23,file='input_vgy.dat',status='unknown')
          do i=1,numl,2
              read(23,*)  vgyg(i)
              vgyg(i+1) = vgyg(i)
          end do
      end if
      close(unit=23)

********
* vgzg *
********
      if (vgzsolv.eq.0) then
          open(unit=23,file='input_vgz.dat',status='unknown')
          do i=1,numl,2
              read(23,*)  vgzg(i)
              vgzg(i+1) = vgzg(i)
          end do
      end if
      close(unit=23)

*********
* ome1g *
*********
      if (ome1solv.eq.0) then
          open(unit=23,file='input_ome1.dat',status='unknown')
          do i=1,numl,2
              read(23,*)   ome1g(i)
              ome1g(i+1) = ome1g(i)
          end do
      end if
      close(unit=23)

*********
* ome2g *
*********
      if (ome2solv.eq.0) then
          open(unit=23,file='input_ome2.dat',status='unknown')
          do i=1,numl,2
              read(23,*)   ome2g(i)
              ome2g(i+1) = ome2g(i)
          end do
      end if
      close(unit=23)

*********
* ome3g *
*********
      if (ome3solv.eq.0) then
          open(unit=23,file='input_ome3.dat',status='unknown')
          do i=1,numl,2
              read(23,*)   ome3g(i)
              ome3g(i+1) = ome3g(i)
          end do
      end if
      close(unit=23)

************************************************************
************************************************************
************************************************************
************************************************************
* Give geometry and motion of each truncated elliptic cone
************************************************************
************************************************************
************************************************************
************************************************************

****************************
* input body_geometry.dat  *
****************************
      open(unit=21,file='body_geometry.dat',status='unknown')
      do 491 j=1,21
      read (21,*) ra1(j),ra2(j),rb1(j),rb2(j),len(j),rho(j)
  491 continue

      do 493 j=1,9
      read(21,*) bu(j)
  493 continue
      close(21)

* non-dimensionalize gravitational acceleration by human height
      grav=grav/bu(8)



***********body geometry for anybody******************************

      if (flagoutm.eq.1) then

       call anygeometry(len,bu)

c      call anygeometry(len(1),len(2),len(3),len(4),len(5),len(6)
c     &                ,len(7),len(10),len(12),len(14),len(16),len(18)
c     &                ,bu(5),bu(6),bu(8),bu(9))

      end if
*************************************************************************



************************************************
* obtain mass center and mass for each segment
* based on coordinate of the segment itself
************************************************
      div0 = 100
      
      do 150 j=1,numb

      m(j) =0.0d0
      z = 0.0d0
      mz = 0.0d0
      do 100 i=1,div0
      d3 = len(j) /dble(div0)
      z = len(j)/dble(div0) * (dble(i)-0.5d0)
      
* obtain whole mass by numerically integrating thin elliptic plates
      r1 = ra1(j) + (rb1(j)-ra1(j))*(dble(i)-0.5d0)/dble(div0)
      r2 = ra2(j) + (rb2(j)-ra2(j))*(dble(i)-0.5d0)/dble(div0)
      dm = pi * r1 * r2 * d3 * rho(j)
      m(j) = m(j) + dm
      mz = mz + dm * z
  100 continue
      leng(j) = mz / m(j)

*************************************************************
* obtain principal moment of inertia about its mass center
* for each segment 
*************************************************************
      
      i1(j) = 0.0d0
      i2(j) = 0.0d0
      i3(j) = 0.0d0

      div0 = 100
      do 110 i=1,div0
      d3 = len(j) /dble(div0)
      z = len(j)/dble(div0) * (dble(i)-0.5d0)

      r1 = ra1(j) + (rb1(j)-ra1(j))*(dble(i)-0.5d0)/dble(div0)
      r2 = ra2(j) + (rb2(j)-ra2(j))*(dble(i)-0.5d0)/dble(div0)

      dm = pi * r1 * r2 * d3 * rho(j)

      di1 = pi * r1 * r2 **3 *d3*rho(j)*0.25d0 + dm * (z-leng(j))**2
      di2 = pi * r2 * r1 **3 *d3*rho(j)*0.25d0 + dm * (z-leng(j))**2
      di3 = di1 + di2 - 2.0d0*dm*(z-leng(j))**2


      i1(j) = i1(j) + di1
      i2(j) = i2(j) + di2
      i3(j) = i3(j) + di3

  110 continue

  150 continue


*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
* Quantities below here vary according to time step
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************
*****************************************************


********************************************
* give position and angle of each segment
* accoding to base segment cooridinate
********************************************

      do 200 k=0,numl2

      open(unit=22,file='joint_motion.dat',status='unknown')
      read(22,*) dummy

      do 195 j=1,numsp
      tsp(j)= dble(j-1)/dble(numsp)
  195 continue

* first, give angle of all elements
      do 201 j=1,numb
      e11x(j) = 1.0d0
      e11y(j) = 0.0d0
      e11z(j) = 0.0d0

      e12x(j) = 0.0d0
      e12y(j) = 1.0d0
      e12z(j) = 0.0d0

      e13x(j) = 0.0d0
      e13y(j) = 0.0d0
      e13z(j) = 1.0d0
  201 continue

* make vector from point-b at upper breast to shoulder joint
      xs2ra = 0.0d0
      ys2ra = -rb2(4)
      zs2ra = -bu(2)

      xs2la = 0.0d0
      ys2la = +rb2(4)
      zs2la = -bu(2)

* make vector from point-b at neck to point-a at head
      xnb2ha = bu(3)
      ynb2ha = 0.0d0
      znb2ha = 0.0d0

* give position of hip joint against point-b at hip
      xh2rl = 0.0d0
      yh2rl = 0.0d0 -bu(4)
      zh2rl = 0.0d0 -bu(5)
      xh2ll = 0.0d0
      yh2ll = 0.0d0 +bu(4)
      zh2ll = 0.0d0 -bu(5)

* make vector from right thigh joint to point-a at thigh
      xtj2tra = 0.0d0
      ytj2tra = 0.0d0
      ztj2tra = +bu(5)

* make vector from left thigh joint to point-a at thigh
      xtj2tla = 0.0d0
      ytj2tla = 0.0d0
      ztj2tla = +bu(5)

* give position of foot joint against point-b at right shank
      xsh2rf = 0.0d0
      ysh2rf = 0.0d0
      zsh2rf = 0.0d0+bu(6)

* give position of foot joint against point-b at left shank
      xsh2lf = 0.0d0
      ysh2lf = 0.0d0
      zsh2lf = 0.0d0 +bu(6)

* make vector from right foot joint to point-a
      xfj2ra = 0.0d0
      yfj2ra = 0.0d0
      zfj2ra = -bu(6)

* make vector from left foot joint to point-a
      xfj2la = 0.0d0
      yfj2la = 0.0d0
      zfj2la = -bu(6)

* make vector from right upper arm joint to point-a
      xuaj2ra = 0.0d0
      yuaj2ra = -bu(1)
      zuaj2ra = -bu(2)

* make vector from left upper arm joint to point-a
      xuaj2la = 0.0d0
      yuaj2la = bu(1)
      zuaj2la = -bu(2)

* make lower half body turn downward
      do 205 j=8,15
      if (j.eq.8) then
          theta = pi+bu(7)
      else if (j.eq.9) then
          theta = pi-bu(7)
          call vecrot(xh2rl,yh2rl,zh2rl,0.0d0,theta,0.0d0)
          call vecrot(xh2ll,yh2ll,zh2ll,0.0d0,theta,0.0d0)
      else if (j.eq.10) then
          theta = pi
          call vecrot(xtj2tra,ytj2tra,ztj2tra,0.0d0,theta,0.0d0)
      else if (j.eq.11) then
          theta = pi
          call vecrot(xtj2tla,ytj2tla,ztj2tla,0.0d0,theta,0.0d0)
      else if (j.eq.12) then
          theta = pi
          call vecrot(xsh2rf,ysh2rf,zsh2rf,0.0d0,theta,0.0d0)
      else if (j.eq.13) then
          theta = pi
          call vecrot(xsh2lf,ysh2lf,zsh2lf,0.0d0,theta,0.0d0)
      else if (j.eq.14) then
          theta = pi
          call vecrot(xfj2ra,yfj2ra,zfj2ra,0.0d0,theta,0.0d0)
      else if (j.eq.15) then
          theta = pi
          call vecrot(xfj2la,yfj2la,zfj2la,0.0d0,theta,0.0d0)
      else
          theta = pi
      end if
      if(j.le.numb) then
      call vecrot(e11x(j),e11y(j),e11z(j),0.0d0,theta,0.0d0)
      call vecrot(e12x(j),e12y(j),e12z(j),0.0d0,theta,0.0d0)
      call vecrot(e13x(j),e13y(j),e13z(j),0.0d0,theta,0.0d0)
      end if
  205 continue


* start data reading
      do 202 jl=1,10000
      read(22,*) j,ja
      if (j.eq.0) then
*********** count number of joint motion data******************
         jcount = jl-1
***************************************************************
          goto 203
      end if

      do 204 jj=1,numsp
      read(22,*) uspval
      usp(jj)= uspval *pi/180.0d0
  204 continue

      call spline(numsp,tsp,usp,k,numl,theta)

      if (ja.eq.1) then
          call vecrot(e11x(j),e11y(j),e11z(j),theta,0.0d0,0.0d0)
          call vecrot(e12x(j),e12y(j),e12z(j),theta,0.0d0,0.0d0)
          call vecrot(e13x(j),e13y(j),e13z(j),theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(e11x(j),e11y(j),e11z(j),0.0d0,theta,0.0d0)
          call vecrot(e12x(j),e12y(j),e12z(j),0.0d0,theta,0.0d0)
          call vecrot(e13x(j),e13y(j),e13z(j),0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(e11x(j),e11y(j),e11z(j),0.0d0,0.0d0,theta)
          call vecrot(e12x(j),e12y(j),e12z(j),0.0d0,0.0d0,theta)
          call vecrot(e13x(j),e13y(j),e13z(j),0.0d0,0.0d0,theta)
      end if

      if (j.eq.4) then
********************
* (4) upper breast *
********************
* rotate vector toward shoulder joint
      if (ja.eq.1) then
          call vecrot(xs2ra,ys2ra,zs2ra,theta,0.0d0,0.0d0)
          call vecrot(xs2la,ys2la,zs2la,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xs2ra,ys2ra,zs2ra,0.0d0,theta,0.0d0)
          call vecrot(xs2la,ys2la,zs2la,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xs2ra,ys2ra,zs2ra,0.0d0,0.0d0,theta)
          call vecrot(xs2la,ys2la,zs2la,0.0d0,0.0d0,theta)
      else if (ja.eq.4) then
          call vecrot(xs2ra,ys2ra,zs2ra,theta,0.0d0,0.0d0)
      else if (ja.eq.5) then
          call vecrot(xs2ra,ys2ra,zs2ra,0.0d0,theta,0.0d0)
      else if (ja.eq.6) then
          call vecrot(xs2ra,ys2ra,zs2ra,0.0d0,0.0d0,theta)
      else if (ja.eq.7) then
          call vecrot(xs2la,ys2la,zs2la,theta,0.0d0,0.0d0)
      else if (ja.eq.8) then
          call vecrot(xs2la,ys2la,zs2la,0.0d0,theta,0.0d0)
      else if (ja.eq.9) then
          call vecrot(xs2la,ys2la,zs2la,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.7) then
************
* (7) head *
************
* rotate vector from joint to point-a together
      if (ja.eq.1) then
          call vecrot(xnb2ha,ynb2ha,znb2ha,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xnb2ha,ynb2ha,znb2ha,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xnb2ha,ynb2ha,znb2ha,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.9) then
*****************
* (9) lower hip *
*****************
      if (ja.eq.1) then
          call vecrot(xh2rl,yh2rl,zh2rl,theta,0.0d0,0.0d0)
          call vecrot(xh2ll,yh2ll,zh2ll,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xh2rl,yh2rl,zh2rl,0.0d0,theta,0.0d0)
          call vecrot(xh2ll,yh2ll,zh2ll,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xh2rl,yh2rl,zh2rl,0.0d0,0.0d0,theta)
          call vecrot(xh2ll,yh2ll,zh2ll,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.10) then
********************
* (10) right thigh *
********************
      if (ja.eq.1) then
          call vecrot(xtj2tra,ytj2tra,ztj2tra,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xtj2tra,ytj2tra,ztj2tra,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xtj2tra,ytj2tra,ztj2tra,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.11) then
*******************
* (11) left thigh *
*******************
      if (ja.eq.1) then
          call vecrot(xtj2tla,ytj2tla,ztj2tla,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xtj2tla,ytj2tla,ztj2tla,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xtj2tla,ytj2tla,ztj2tla,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.12) then
********************
* (12) right shank *
********************
      if (ja.eq.1) then
          call vecrot(xsh2rf,ysh2rf,zsh2rf,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xsh2rf,ysh2rf,zsh2rf,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xsh2rf,ysh2rf,zsh2rf,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.13) then
*******************
* (13) left shank *
*******************
      if (ja.eq.1) then
          call vecrot(xsh2lf,ysh2lf,zsh2lf,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xsh2lf,ysh2lf,zsh2lf,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xsh2lf,ysh2lf,zsh2lf,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.14) then
*******************
* (14) right foot *
*******************
      if (ja.eq.1) then
          call vecrot(xfj2ra,yfj2ra,zfj2ra,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xfj2ra,yfj2ra,zfj2ra,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xfj2ra,yfj2ra,zfj2ra,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.15) then
******************
* (15) left foot *
******************
      if (ja.eq.1) then
          call vecrot(xfj2la,yfj2la,zfj2la,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xfj2la,yfj2la,zfj2la,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xfj2la,yfj2la,zfj2la,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.16) then
************************
* (16) right upper arm *
************************
      if (ja.eq.1) then
          call vecrot(xuaj2ra,yuaj2ra,zuaj2ra,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xuaj2ra,yuaj2ra,zuaj2ra,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xuaj2ra,yuaj2ra,zuaj2ra,0.0d0,0.0d0,theta)
      end if

      else if (j.eq.17) then
***********************
* (17) left upper arm *
***********************
      if (ja.eq.1) then
          call vecrot(xuaj2la,yuaj2la,zuaj2la,theta,0.0d0,0.0d0)
      else if (ja.eq.2) then
          call vecrot(xuaj2la,yuaj2la,zuaj2la,0.0d0,theta,0.0d0)
      else if (ja.eq.3) then
          call vecrot(xuaj2la,yuaj2la,zuaj2la,0.0d0,0.0d0,theta)
      end if

      end if
  202 continue

  203 continue


      do 206 j=1,numb
      if (j.eq.1) then
*******************
* (1) lower waist *
*******************
      x1a(j) = 0.0d0
      y1a(j) = 0.0d0
      z1a(j) = 0.0d0
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if ((j.eq.2).or.(j.eq.3).or.(j.eq.5).or.(j.eq.6)) then
*************************************************************
* (2) upper waist, (3) lower breast, (5) shoulder, (6) neck *
*************************************************************
      x1a(j) = x1b(j-1)
      y1a(j) = y1b(j-1)
      z1a(j) = z1b(j-1)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.4) then
********************
* (4) upper breast *
********************
      x1a(j) = x1b(j-1)
      y1a(j) = y1b(j-1)
      z1a(j) = z1b(j-1)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.7) then
************
* (7) head *
************
      x1a(j) = x1b(j-1) + xnb2ha
      y1a(j) = y1b(j-1) + ynb2ha
      z1a(j) = z1b(j-1) + znb2ha
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.8) then
*****************
* (8) upper hip *
*****************
      x1a(j) = x1a(1)
      y1a(j) = y1a(1)
      z1a(j) = z1a(1)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.9) then
*****************
* (9) lower hip *
*****************
      x1a(j) = x1b(j-1)
      y1a(j) = y1b(j-1)
      z1a(j) = z1b(j-1)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.10) then
********************
* (10) right thigh *
********************
      x1a(j) = x1b(9) + xh2rl + xtj2tra
      y1a(j) = y1b(9) + yh2rl + ytj2tra
      z1a(j) = z1b(9) + zh2rl + ztj2tra
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.11) then
*******************
* (11) left thigh *
*******************
      x1a(j) = x1b(9) + xh2ll + xtj2tla
      y1a(j) = y1b(9) + yh2ll + ytj2tla
      z1a(j) = z1b(9) + zh2ll + ztj2tla
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.12) then
********************
* (12) right shank *
********************
      x1a(j) = x1b(j-2)
      y1a(j) = y1b(j-2)
      z1a(j) = z1b(j-2)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.13) then
*******************
* (13) left shank *
*******************
      x1a(j) = x1b(j-2)
      y1a(j) = y1b(j-2)
      z1a(j) = z1b(j-2)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.14) then
*******************
* (14) right foot *
*******************
      x1a(j) = x1b(j-2)+xfj2ra+xsh2rf
      y1a(j) = y1b(j-2)+yfj2ra+ysh2rf
      z1a(j) = z1b(j-2)+zfj2ra+zsh2rf
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.15) then
******************
* (15) left foot *
******************
      x1a(j) = x1b(j-2)+xfj2la+xsh2lf
      y1a(j) = y1b(j-2)+yfj2la+ysh2lf
      z1a(j) = z1b(j-2)+zfj2la+zsh2lf
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.16) then
************************
* (16) right upper arm *
************************
      x1a(j) = x1b(4) + xs2ra + xuaj2ra
      y1a(j) = y1b(4) + ys2ra + yuaj2ra
      z1a(j) = z1b(4) + zs2ra + zuaj2ra
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if (j.eq.17) then
***********************
* (17) left upper arm *
***********************
      x1a(j) = x1b(4) + xs2la + xuaj2la
      y1a(j) = y1b(4) + ys2la + yuaj2la
      z1a(j) = z1b(4) + zs2la + zuaj2la
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)

      else if ((j.eq.18).or.(j.eq.19).or.(j.eq.20).or.(j.eq.21)) then
***********************************************************************
* (18) right forearm,(19) left forearm,(20) right hand,(21) left hand *
***********************************************************************
      x1a(j) = x1b(j-2)
      y1a(j) = y1b(j-2)
      z1a(j) = z1b(j-2)
      x1b(j) = x1a(j) + len(j) * e13x(j)
      y1b(j) = y1a(j) + len(j) * e13y(j)
      z1b(j) = z1a(j) + len(j) * e13z(j)
      end if
  206 continue


*******************************************************
* calculate position for mass center of each segment
* at base coordinate
*******************************************************
      do 209 j=1,numb
      x1g(j) = x1a(j) + e13x(j) * leng(j) 
      y1g(j) = y1a(j) + e13y(j) * leng(j)
      z1g(j) = z1a(j) + e13z(j) * leng(j)
  209 continue
 
******************************************************
* calculate position and mass for total mass center
* at base coordinate
******************************************************
      mall = 0.0d0
      x1mall = 0.0d0
      y1mall = 0.0d0
      z1mall = 0.0d0

      do 210 j=1,numb
      x1mall = x1mall + x1g(j) * m(j)
      y1mall = y1mall + y1g(j) * m(j)
      z1mall = z1mall + z1g(j) * m(j)

      mall = mall + m(j)
  210 continue
      x1gall(k) = x1mall / mall
      y1gall(k) = y1mall / mall
      z1gall(k) = z1mall / mall

**************************************************************
* calculate inertia tensor and sum up them
* at a coordinate which is parallel to base coordinate and
* whose origin is located at total mass center
**************************************************************

* coordinate transform principal moment of inertia into
* inertia tensor about each mass center of each segment
* and about axes parallel to the coordinate above. 

* first, reset matrix
      call matset0(inerall)

      do 220 j=1,numb
      ten(1,1) = e11x(j)
      ten(2,1) = e11y(j)
      ten(3,1) = e11z(j)
      ten(1,2) = e12x(j)
      ten(2,2) = e12y(j)
      ten(3,2) = e12z(j)
      ten(1,3) = e13x(j)
      ten(2,3) = e13y(j)
      ten(3,3) = e13z(j)

      tent(1,1) = ten(1,1) * i1(j)
      tent(2,1) = ten(1,2) * i2(j)
      tent(3,1) = ten(1,3) * i3(j)
      tent(1,2) = ten(2,1) * i1(j)
      tent(2,2) = ten(2,2) * i2(j)
      tent(3,2) = ten(2,3) * i3(j)
      tent(1,3) = ten(3,1) * i1(j)
      tent(2,3) = ten(3,2) * i2(j)
      tent(3,3) = ten(3,3) * i3(j)

* obtain transformed inertia tensor 
* by multiplying coordinate transformation matrix
      call matpro(ten,tent,iner)

* add up to whole inertia tensor
      call matadd(inerall,iner)

* in addition, sum moment of inertia due to distance of mass center
      inerall(1,1) = inerall(1,1)
     &       + ( (y1g(j)-y1gall(k))**2 + (z1g(j)-z1gall(k))**2 ) * m(j)
      inerall(2,1) = inerall(2,1)
     &       - ( (x1g(j)-x1gall(k)) * (y1g(j)-y1gall(k)) ) * m(j)
      inerall(3,1) = inerall(3,1)
     &       - ( (x1g(j)-x1gall(k)) * (z1g(j)-z1gall(k)) ) * m(j)
      inerall(1,2) = inerall(1,2)
     &       - ( (x1g(j)-x1gall(k)) * (y1g(j)-y1gall(k)) ) * m(j)
      inerall(2,2) = inerall(2,2)
     &       + ( (x1g(j)-x1gall(k))**2 + (z1g(j)-z1gall(k))**2 ) * m(j)
      inerall(3,2) = inerall(3,2)
     &       - ( (y1g(j)-y1gall(k)) * (z1g(j)-z1gall(k)) ) * m(j)
      inerall(1,3) = inerall(1,3)
     &       - ( (x1g(j)-x1gall(k)) * (z1g(j)-z1gall(k)) ) * m(j)
      inerall(2,3) = inerall(2,3)
     &       - ( (y1g(j)-y1gall(k)) * (z1g(j)-z1gall(k)) ) * m(j)
      inerall(3,3) = inerall(3,3)
     &       + ( (x1g(j)-x1gall(k))**2 + (y1g(j)-y1gall(k))**2 ) * m(j)

  220 continue

********************************************************
* obtain principal axes and principal moment of inertia
* by solving eigenvalue problem
********************************************************

* first, make data single precision
      call matsubs(inerall,ineralls)
      
* call Slatec library and solve eigenvalue problem
      lda = 3
      n = 3
      ldv = 3
      job = 1
      call sgeev(ineralls,lda,n,evalc,evecc,ldv,work,job,info)

*     first, make eigenvectors double precision and normalize them
      call matdble(evecc,mata)

      call mattrans(mata,matb)
      call matpro(matb,inerall,matc)
      call matpro(matc,mata,matb)

      if (k.eq.0) then
*     at the first time step, obtain an axis closest to zb
*     and make it axis-1 (new algorithm since 4.0.0)
          evalmax = -100000000000000.0d0
          imax = -1
          do 225 j2=1,3
          inproval = dabs(inpro(mata(1,j2),mata(2,j2),mata(3,j2)
     &                    ,0.0d0, 0.0d0, 1.0d0))
          if (inproval.gt.evalmax) then
              imax = j2
              evalmax = inproval
          end if
  225     continue
          imax1 = imax
          e11xall(k) = mata(1,imax)
          e11yall(k) = mata(2,imax)
          e11zall(k) = mata(3,imax)
          iall1(k) = dble(evalc(imax))

*     at the first time step, obtain minimum of eigenvalue,
*     and make it axis-1 (old algorithm by 4.0.0)
c         evalmin = 100000000000000.0d0
c         imin = -1
c         do 225 j2=1,3
c         if (dble(evalc(j2)).lt.evalmin) then
c             imin = j2
c             evalmin = evalc(j2)
c         end if
c 225     continue
c         e11xall(k) = mata(1,imin)
c         e11yall(k) = mata(2,imin)
c         e11zall(k) = mata(3,imin)
c         iall1(k) = dble(evalc(imin))

*     using inner product of axis-1 and zb, adjust direction
          tmpinpro = inpro(e11xall(k),e11yall(k),e11zall(k)
     &                     ,0.0d0,0.0d0,1.0d0)
          if (tmpinpro.lt.0.0d0) then
              mata(1,1) = -mata(1,1)
              mata(2,1) = -mata(2,1)
              mata(3,1) = -mata(3,1)
              mata(1,2) = -mata(1,2)
              mata(2,2) = -mata(2,2)
              mata(3,2) = -mata(3,2)
              mata(1,3) = -mata(1,3)
              mata(2,3) = -mata(2,3)
              mata(3,3) = -mata(3,3)
              e11xall(k) = -e11xall(k)
              e11yall(k) = -e11yall(k)
              e11zall(k) = -e11zall(k)
          end if

*     set axis close to yb as axis-2 
          iabsmax = 0.0d0
          inpromax = 0.0d0
          imax = -1
          do 227 j2=1,3
          tmpinpro = inpro(mata(1,j2),mata(2,j2),mata(3,j2)
     &                     ,0.0d0,1.0d0,0.0d0)
          tmpabs = dabs(tmpinpro)
           
          if (tmpabs.gt.iabsmax) then
          if (j2.ne.imax1) then
              imax = j2
              iabsmax = tmpabs
              inpromax = tmpinpro
          end if
          end if
  227     continue
          imax2 = imax
          if (inpromax.gt.0.0d0) then
              e12xall(k) = mata(1,imax)
              e12yall(k) = mata(2,imax)
              e12zall(k) = mata(3,imax)
              iall2(k) = dble(evalc(imax))
          else
              e12xall(k) = - mata(1,imax)
              e12yall(k) = - mata(2,imax)
              e12zall(k) = - mata(3,imax)
              iall2(k) = dble(evalc(imax))
          end if

*     set axis close to -xb as axis-3
          iabsmax = 0.0d0
          inpromax = 0.0d0
          imax = -1
          do 228 j2=1,3
          tmpinpro = inpro(mata(1,j2),mata(2,j2),mata(3,j2)
     &                     ,-1.0d0,0.0d0,0.0d0)
          tmpabs = dabs(tmpinpro)
           
          if (tmpabs.gt.iabsmax) then
          if (j2.ne.imax1) then
          if (j2.ne.imax2) then
              imax = j2
              iabsmax = tmpabs
              inpromax = tmpinpro
          end if
          end if
          end if
  228     continue
          if (inpromax.gt.0.0d0) then
              e13xall(k) = mata(1,imax)
              e13yall(k) = mata(2,imax)
              e13zall(k) = mata(3,imax)
              iall3(k) = dble(evalc(imax))
          else
              e13xall(k) = - mata(1,imax)
              e13yall(k) = - mata(2,imax)
              e13zall(k) = - mata(3,imax)
              iall3(k) = dble(evalc(imax))
          end if

      else

*     for the steps except the first step, choose axis close to 
*     the former one

*     first, make them be 0 or -1
          e11xall(k) = 0.0d0
          e11yall(k) = 0.0d0
          e11zall(k) = 0.0d0
          iall1(k)   = -1.0d0
          e12xall(k) = 0.0d0
          e12yall(k) = 0.0d0
          e12zall(k) = 0.0d0
          iall2(k)   = -1.0d0
          e13xall(k) = 0.0d0
          e13yall(k) = 0.0d0
          e13zall(k) = 0.0d0
          iall3(k)   = -1.0d0
          do 226 j2=1,3

          mata1=mata(1,j2)
          mata2=mata(2,j2)
          mata3=mata(3,j2)

          inpro1 = inpro(e11xall(k-1),e11yall(k-1),e11zall(k-1)
     &                     ,mata1,mata2,mata3)
          inpro2 = inpro(e12xall(k-1),e12yall(k-1),e12zall(k-1)
     &                     ,mata1,mata2,mata3)
          inpro3 = inpro(e13xall(k-1),e13yall(k-1),e13zall(k-1)
     &                     ,mata1,mata2,mata3)
          abinpro1 = dabs(inpro1)
          abinpro2 = dabs(inpro2)
          abinpro3 = dabs(inpro3)
          if      ((abinpro1.ge.abinpro2).and.
     &             (abinpro1.ge.abinpro3)) then
*************************************************
          if (inpro1.gt.0.0d0) then
              e11xall(k) =   mata1
              e11yall(k) =   mata2
              e11zall(k) =   mata3
              iall1(k) = dble(evalc(j2))
          else if (inpro1.lt.-0.0d0) then
              e11xall(k) =   -mata1
              e11yall(k) =   -mata2
              e11zall(k) =   -mata3
              iall1(k) = dble(evalc(j2))
          end if
*************************************************
          else if ((abinpro2.ge.abinpro1).and.
     &             (abinpro2.ge.abinpro3)) then
*************************************************
          if (inpro2.gt.0.0d0) then
              e12xall(k) =   mata1
              e12yall(k) =   mata2
              e12zall(k) =   mata3
              iall2(k) = dble(evalc(j2))
          else if (inpro2.lt.-0.0d0) then
              e12xall(k) =   -mata1
              e12yall(k) =   -mata2
              e12zall(k) =   -mata3
              iall2(k) = dble(evalc(j2))
          end if
*************************************************
          else if ((abinpro3.ge.abinpro1).and.
     &             (abinpro3.ge.abinpro2)) then
*************************************************
          if (inpro3.gt.0.0d0) then
              e13xall(k) =   mata1
              e13yall(k) =   mata2
              e13zall(k) =   mata3
              iall3(k) = dble(evalc(j2))
          else if (inpro3.lt.-0.0d0) then
              e13xall(k) =   -mata1
              e13yall(k) =   -mata2
              e13zall(k) =   -mata3
              iall3(k) = dble(evalc(j2))
          end if
*************************************************
          else 
                  write(*,*) 'No way! 6'
                  write(*,*) iall1(k)
                  write(*,*) iall2(k)
                  write(*,*) iall3(k)
                  write(*,*) inpro1
                  write(*,*) inpro2
                  write(*,*) inpro3
                  write(*,*) e11xall(k),e11yall(k),e11zall(k)
                  write(*,*) e12xall(k),e12yall(k),e12zall(k)
                  write(*,*) e13xall(k),e13yall(k),e13zall(k)
                  write(*,*) e11xall(k-1),e11yall(k-1),e11zall(k-1)
                  write(*,*) e12xall(k-1),e12yall(k-1),e12zall(k-1)
                  write(*,*) e13xall(k-1),e13yall(k-1),e13zall(k-1)
                  write(*,*) mata(1,1),mata(2,1),mata(3,1)
                  write(*,*) mata(1,2),mata(2,2),mata(3,2)
                  write(*,*) mata(1,3),mata(2,3),mata(3,3)
                  stop
          end if

c          inprotmp = inpro(e11xall(k-1),e11yall(k-1),e11zall(k-1)
c     &                     ,mata1,mata2,mata3)
cc         write(*,*) ' inprotmp1 ', inprotmp
c          if (inprotmp.gt.0.9d0) then
c              e11xall(k) =   mata1
c              e11yall(k) =   mata2
c              e11zall(k) =   mata3
c              iall1(k) = dble(evalc(j2))
c          else if (inprotmp.lt.-0.9d0) then
c              e11xall(k) =   -mata1
c              e11yall(k) =   -mata2
c              e11zall(k) =   -mata3
c              iall1(k) = dble(evalc(j2))
c          end if
c
c          inprotmp = inpro(e12xall(k-1),e12yall(k-1),e12zall(k-1)
c     &                     ,mata1,mata2,mata3)
cc         write(*,*) ' inprotmp2 ', inprotmp
c          if (inprotmp.gt.0.9d0) then
c              e12xall(k) =   mata1
c              e12yall(k) =   mata2
c              e12zall(k) =   mata3
c              iall2(k) = dble(evalc(j2))
c          else if (inprotmp.lt.-0.9d0) then
c              e12xall(k) =   -mata1
c              e12yall(k) =   -mata2
c              e12zall(k) =   -mata3
c              iall2(k) = dble(evalc(j2))
c          end if
c
c          inprotmp = inpro(e13xall(k-1),e13yall(k-1),e13zall(k-1)
c     &                     ,mata1,mata2,mata3)
cc         write(*,*) ' inprotmp3 ', inprotmp
c          if (inprotmp.gt.0.9d0) then
c              e13xall(k) =   mata1
c              e13yall(k) =   mata2
c              e13zall(k) =   mata3
c              iall3(k) = dble(evalc(j2))
c          else if (inprotmp.lt.-0.9d0) then
c              e13xall(k) =   -mata1
c              e13yall(k) =   -mata2
c              e13zall(k) =   -mata3
c              iall3(k) = dble(evalc(j2))
c          end if

  226 continue

*     check
          if ( (iall1(k).lt.0.0d0).or.
     &         (iall2(k).lt.0.0d0).or.
     &         (iall3(k).lt.0.0d0) ) then
          if   (comsolv.ne.0) then
              write(*,*) 'No way! 5'
              write(*,*) k
              write(*,*) iall1(k)
              write(*,*) iall2(k)
              write(*,*) iall3(k)
              write(*,*) e11xall(k),e11yall(k),e11zall(k)
              write(*,*) e12xall(k),e12yall(k),e12zall(k)
              write(*,*) e13xall(k),e13yall(k),e13zall(k)
              write(*,*) e11xall(k-1),e11yall(k-1),e11zall(k-1)
              write(*,*) e12xall(k-1),e12yall(k-1),e12zall(k-1)
              write(*,*) e13xall(k-1),e13yall(k-1),e13zall(k-1)
              write(*,*) mata(1,1),mata(2,1),mata(3,1)
              write(*,*) mata(1,2),mata(2,2),mata(3,2)
              write(*,*) mata(1,3),mata(2,3),mata(3,3)
              write(*,*) ''
              write(*,*) evecc(1,1), evecc(2,1), evecc(3,1)
              write(*,*) evecc(1,2), evecc(2,2), evecc(3,2)
              write(*,*) evecc(1,3), evecc(2,3), evecc(3,3)
              stop
          end if
          end if

      end if

* in the case when center of mass and principal axes are not computed
      if(comsolv.eq.0) then
          x1gall(k) = 0.0d0
          y1gall(k) = 0.0d0
          z1gall(k) = 0.0d0

          e11xall(k) = 1.0d0
          e11yall(k) = 0.0d0
          e11zall(k) = 0.0d0

          e12xall(k) = 0.0d0
          e12yall(k) = 1.0d0
          e12zall(k) = 0.0d0

          e13xall(k) = 0.0d0
          e13yall(k) = 0.0d0
          e13zall(k) = 1.0d0
      end if

******************************************************************
* transform each coordinate in the base cooridinate 
* into a coordinate whose origin is total mass center and
* whose axes are parallel to principal axes of moment of inertia
******************************************************************

      do 230 j=1,numb
* position
      xia(j,k) = e11xall(k) * (x1a(j)-x1gall(k))
     &         + e11yall(k) * (y1a(j)-y1gall(k))
     &         + e11zall(k) * (z1a(j)-z1gall(k))
      yia(j,k) = e12xall(k) * (x1a(j)-x1gall(k))
     &         + e12yall(k) * (y1a(j)-y1gall(k))
     &         + e12zall(k) * (z1a(j)-z1gall(k))
      zia(j,k) = e13xall(k) * (x1a(j)-x1gall(k))
     &         + e13yall(k) * (y1a(j)-y1gall(k))
     &         + e13zall(k) * (z1a(j)-z1gall(k))

      xib(j,k) = e11xall(k) * (x1b(j)-x1gall(k))
     &         + e11yall(k) * (y1b(j)-y1gall(k))
     &         + e11zall(k) * (z1b(j)-z1gall(k))
      yib(j,k) = e12xall(k) * (x1b(j)-x1gall(k))
     &         + e12yall(k) * (y1b(j)-y1gall(k))
     &         + e12zall(k) * (z1b(j)-z1gall(k))
      zib(j,k) = e13xall(k) * (x1b(j)-x1gall(k))
     &         + e13yall(k) * (y1b(j)-y1gall(k))
     &         + e13zall(k) * (z1b(j)-z1gall(k))

* transform principal axes of moment of inertia of each segment
* into whole principal axes coordinate
      ei1x(j,k) = e11xall(k) * e11x(j)
     &          + e11yall(k) * e11y(j)
     &          + e11zall(k) * e11z(j)
      ei1y(j,k) = e12xall(k) * e11x(j)
     &          + e12yall(k) * e11y(j)
     &          + e12zall(k) * e11z(j)
      ei1z(j,k) = e13xall(k) * e11x(j)
     &          + e13yall(k) * e11y(j)
     &          + e13zall(k) * e11z(j)
       
      ei2x(j,k) = e11xall(k) * e12x(j)
     &          + e11yall(k) * e12y(j)
     &          + e11zall(k) * e12z(j)
      ei2y(j,k) = e12xall(k) * e12x(j)
     &          + e12yall(k) * e12y(j)
     &          + e12zall(k) * e12z(j)
      ei2z(j,k) = e13xall(k) * e12x(j)
     &          + e13yall(k) * e12y(j)
     &          + e13zall(k) * e12z(j)

      ei3x(j,k) = e11xall(k) * e13x(j)
     &          + e11yall(k) * e13y(j)
     &          + e11zall(k) * e13z(j)
      ei3y(j,k) = e12xall(k) * e13x(j)
     &          + e12yall(k) * e13y(j)
     &          + e12zall(k) * e13z(j)
      ei3z(j,k) = e13xall(k) * e13x(j)
     &          + e13yall(k) * e13y(j)
     &          + e13zall(k) * e13z(j)

  230 continue

      close(22)
  200 continue


***********convert joint motion for anybody******************************
      if (flagoutm.eq.1) then

       call outanyjoint(numb,dummy,jcount,pi,ttt)

      end if
**************************************************************************


*********************************************************************
* optain velocity in coordinate whose origin is mass center
* and whose axes are principal axes of moment of inertia in advance
*********************************************************************
      do 250 k=1,numl1
      do 250 j=1,numb
      vxia(j,k) = (xia(j,k+1)-xia(j,k-1))/(2.0d0*dt)
      vyia(j,k) = (yia(j,k+1)-yia(j,k-1))/(2.0d0*dt)
      vzia(j,k) = (zia(j,k+1)-zia(j,k-1))/(2.0d0*dt)
      vxib(j,k) = (xib(j,k+1)-xib(j,k-1))/(2.0d0*dt)
      vyib(j,k) = (yib(j,k+1)-yib(j,k-1))/(2.0d0*dt)
      vzib(j,k) = (zib(j,k+1)-zib(j,k-1))/(2.0d0*dt)
  250 continue

*********************************************************************
* optain acceleration in coordinate whose origin is mass center
* and whose axes are principal axes of moment of inertia in advance
*********************************************************************
      do 260 k=1,numl1
      do 260 j=1,numb
      axia(j,k) = (xia(j,k+1)-2.0d0*xia(j,k)+xia(j,k-1))/(dt**2)
      ayia(j,k) = (yia(j,k+1)-2.0d0*yia(j,k)+yia(j,k-1))/(dt**2)
      azia(j,k) = (zia(j,k+1)-2.0d0*zia(j,k)+zia(j,k-1))/(dt**2)
      axib(j,k) = (xib(j,k+1)-2.0d0*xib(j,k)+xib(j,k-1))/(dt**2)
      ayib(j,k) = (yib(j,k+1)-2.0d0*yib(j,k)+yib(j,k-1))/(dt**2)
      azib(j,k) = (zib(j,k+1)-2.0d0*zib(j,k)+zib(j,k-1))/(dt**2)
  260 continue

*****************************************************
* obtain time change of principal moment of inertia
*****************************************************
      do 270 k=1,numl1
      idall1(k) = (iall1(k+1)-iall1(k-1))/(2.0d0*dt)
      idall2(k) = (iall2(k+1)-iall2(k-1))/(2.0d0*dt)
      idall3(k) = (iall3(k+1)-iall3(k-1))/(2.0d0*dt)
  270 continue

********************************************
* give initial value of total mass center
********************************************
      xg = inipos(1)
      yg = inipos(2)
      zg = inipos(3)

      xgloop = xg
      ygloop = yg
      zgloop = zg

******************************************************
* give initial value of direction of principal axes
******************************************************
      e1x = inidirec(1,1)
      e1y = inidirec(1,2)
      e1z = inidirec(1,3)
      call vecrot(e1x,e1y,e1z,0.0d0,0.0d0,-direc)
       
      e2x = inidirec(2,1)
      e2y = inidirec(2,2)
      e2z = inidirec(2,3)
      call vecrot(e2x,e2y,e2z,0.0d0,0.0d0,-direc)

      e3x = inidirec(3,1)
      e3y = inidirec(3,2)
      e3z = inidirec(3,3)
      call vecrot(e3x,e3y,e3z,0.0d0,0.0d0,-direc)

**************************************************
* give initial value of velocity of mass center
**************************************************
      vgx = inivel(1)
      vgy = inivel(2)
      vgz = inivel(3)

****************************************************************
* give initial value of angular velocity about principal axes
****************************************************************
      ome1 =  iniavel(1)
      ome2 =  iniavel(2)
      ome3 =  iniavel(3)

* for "given" case (added 2006/10/10 for version 2.3.0)
      if (vgxsolv.eq.0) then
          vgx = vgxg(1)
      end if
      if (vgysolv.eq.0) then
          vgy = vgyg(1)
      end if
      if (vgzsolv.eq.0) then
          vgz = vgzg(1)
      end if
      if (ome1solv.eq.0) then
          ome1 =  ome1g(1)
      end if
      if (ome2solv.eq.0) then
          ome2 =  ome2g(1)
      end if
      if (ome3solv.eq.0) then
          ome3 =  ome3g(1)
      end if

*initialize inerh!!!miura
      do 99 j=1,numb
      do 99 i=1,div
      do 99 jf = 1,nm
      do 99 ii1 = 1,3
      do 99 ii2 = 1,3
         inerh(j,i,jf,ii1,ii2) = 0.0d0
 99   continue





************************************************************************
************************************************************************
************************************************************************
************************************************************************
************************************************************************
* Compute correction amount due to rotation of principal axes of inertia
************************************************************************
************************************************************************
************************************************************************
************************************************************************
************************************************************************


***********************************************************************
* Calculate position and velocity in the absolute space
* (set the whole mass center as the origin, translational velocity as 0)
***********************************************************************
      omev1all = 0.0d0
      omev2all = 0.0d0
      omev3all = 0.0d0

* -Since ver.4.2.0- below part was modifed 
c     do k=1,numl
      do k=1,numl1
      h1 = 0.0d0
      h2 = 0.0d0
      h3 = 0.0d0
      do j=1,numb

*************************************************
* Compute angular momentum of each body segment *
*   about principal axes of inertia             *
*************************************************

* (1) Due to the translation of the mass cetner

*     Compute position vector of each body segment's mass center
*     (move for leng(j) in the axis 3 of the body segment
*      from the tip-point "a")
      u(1) = xia(j,k) + ei3x(j,k) * leng(j)
      u(2) = yia(j,k) + ei3y(j,k) * leng(j)
      u(3) = zia(j,k) + ei3z(j,k) * leng(j)

*     Compute velocity (momentum) vector
      v(1) = m(j)*(vxia(j,k)*(len(j)-leng(j))/len(j)
     &            +vxib(j,k)*leng(j)/len(j))
      v(2) = m(j)*(vyia(j,k)*(len(j)-leng(j))/len(j)
     &            +vyib(j,k)*leng(j)/len(j))
      v(3) = m(j)*(vzia(j,k)*(len(j)-leng(j))/len(j)
     &            +vzib(j,k)*leng(j)/len(j))

*     Compute angular momentum vector by taking the outer product
*     of the position vector and velocity (momentum) vector
      mvx = exprox(u,v)
      mvy = exproy(u,v)
      mvz = exproz(u,v)

      hm1 = mvx
      hm2 = mvy
      hm3 = mvz

* (2) Due to the moment of inertia
*     First, compute angular velocity of each segment
*       about principal axes of inertia.
*     The axes 1,2,3 are computed from the variation 
*       about 1,2,3 axes, respectively.

*     Axis 1
      dei2x   = (ei2x(j,k+1)-ei2x(j,k-1))/(2.0d0*dt)
      dei2y   = (ei2y(j,k+1)-ei2y(j,k-1))/(2.0d0*dt)
      dei2z   = (ei2z(j,k+1)-ei2z(j,k-1))/(2.0d0*dt)
      om1 = inpro(dei2x,dei2y,dei2z,ei3x(j,k),ei3y(j,k),ei3z(j,k))

*     Axis 2
      dei3x   = (ei3x(j,k+1)-ei3x(j,k-1))/(2.0d0*dt)
      dei3y   = (ei3y(j,k+1)-ei3y(j,k-1))/(2.0d0*dt)
      dei3z   = (ei3z(j,k+1)-ei3z(j,k-1))/(2.0d0*dt)
      om2 = inpro(dei3x,dei3y,dei3z,ei1x(j,k),ei1y(j,k),ei1z(j,k))

*     Axis 3
      dei1x   = (ei1x(j,k+1)-ei1x(j,k-1))/(2.0d0*dt)
      dei1y   = (ei1y(j,k+1)-ei1y(j,k-1))/(2.0d0*dt)
      dei1z   = (ei1z(j,k+1)-ei1z(j,k-1))/(2.0d0*dt)
      om3 = inpro(dei1x,dei1y,dei1z,ei2x(j,k),ei2y(j,k),ei2z(j,k))

*     Compute angular momentum vector
      hi1 = i1(j) * om1 * ei1x(j,k)
     &    + i2(j) * om2 * ei2x(j,k)
     &    + i3(j) * om3 * ei3x(j,k)
      hi2 = i1(j) * om1 * ei1y(j,k)
     &    + i2(j) * om2 * ei2y(j,k)
     &    + i3(j) * om3 * ei3y(j,k)
      hi3 = i1(j) * om1 * ei1z(j,k)
     &    + i2(j) * om2 * ei2z(j,k)
     &    + i3(j) * om3 * ei3z(j,k)

      h1 = h1 + hm1 + hi1
      h2 = h2 + hm2 + hi2
      h3 = h3 + hm3 + hi3

      enddo
******************************************************
* Compute angular velocity, dividing angular momentum 
*  by whole moment of inertia
******************************************************
      omev1(k) = 0.0d0
      omev2(k) = 0.0d0
      omev3(k) = 0.0d0

      if((conam.eq.1).and.(ome1solv.eq.1)) then
        omev1(k) = h1 / iall1(k)
        omev1all = omev1all + omev1(k)
      end if
      if((conam.eq.1).and.(ome2solv.eq.1)) then
        omev2(k) = h2 / iall2(k)
        omev2all = omev2all + omev2(k)
      end if
      if((conam.eq.1).and.(ome3solv.eq.1)) then
        omev3(k) = h3 / iall3(k)
        omev3all = omev3all + omev3(k)
      end if

c     write(*,*) omev1(k), omev2(k), omev3(k)

      enddo

c     write(*,*) 'all', numl,omev1all, omev2all, omev3all

* Finally, correct the angular velocity
* -Since ver.4.2.0- below part was modifed 
c     omev11 =omev1(1)
c     omev21 =omev2(1)
c     omev31 =omev3(1)
c     do k=1,numl1
c     if((conam.eq.1).and.(ome1solv.eq.1)) then
c       omev1(k) = omev1(k) - omev1all/dble(numl)
c       omev1(k) = omev1(k) - omev11
c     end if

c     if((conam.eq.1).and.(ome2solv.eq.1)) then
c       omev2(k) = omev2(k) - omev2all/dble(numl)
c       omev2(k) = omev2(k) - omev21
c     end if

c     if((conam.eq.1).and.(ome3solv.eq.1)) then
c       omev3(k) = omev3(k) - omev3all/dble(numl)
c       omev3(k) = omev3(k) - omev31
c     end if
c     enddo


*******************
* Open data files *
*******************

      call opendata(cycle,bu(8),flagoutr,flagoutf,flagoutj,flagpene,
     &              multi)

****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************
* start time step          *
****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************
****************************

      t = 0.0d0

      roll1 = 0.0d0
      rollb = 0.0d0
      poweraav = 0.0d0
      powerpav = 0.0d0
      powera = 0.0d0
      powerp = 0.0d0
      powtx   = 0.0d0
      pownx   = 0.0d0
      powix   = 0.0d0
      powty   = 0.0d0
      powny   = 0.0d0
      powiy   = 0.0d0
      powtz   = 0.0d0
      pownz   = 0.0d0
      powiz   = 0.0d0
      powtxav = 0.0d0
      pownxav = 0.0d0
      powixav = 0.0d0
      powtyav = 0.0d0
      pownyav = 0.0d0
      powiyav = 0.0d0
      powtzav = 0.0d0
      pownzav = 0.0d0
      powizav = 0.0d0

      do 1000 loop=1,lastloop

      write(*,*) 'cycle= ',loop
      do 1010 kd=1,numl,2
      do 1100 rkloop=1,4
      if (rkloop.eq.1) then
          k=kd
      else if (rkloop.eq.2) then
          k=kd+1
      else if (rkloop.eq.3) then
          k=kd+1
      else if (rkloop.eq.4) then
          k=kd+2
      else
          write(*,*) 'No way!(rkloop)'
          stop
      end if
****************************************************************
* compute absolute position, angle, velocity and acceleration
* (except rigid body component) of each part 
* from present position and angle
****************************************************************
      do 300 j=1,numb
      
* position
      x0a(j) = xg + e1x * xia(j,k)
     &            + e2x * yia(j,k)
     &            + e3x * zia(j,k)
      y0a(j) = yg + e1y * xia(j,k)
     &            + e2y * yia(j,k)
     &            + e3y * zia(j,k)
      z0a(j) = zg + e1z * xia(j,k)
     &            + e2z * yia(j,k)
     &            + e3z * zia(j,k)
      x0ac(j,kd) = x0a(j)
      y0ac(j,kd) = y0a(j)
      z0ac(j,kd) = z0a(j)

      x0b(j) = xg + e1x * xib(j,k)
     &            + e2x * yib(j,k)
     &            + e3x * zib(j,k)
      y0b(j) = yg + e1y * xib(j,k)
     &            + e2y * yib(j,k)
     &            + e3y * zib(j,k)
      z0b(j) = zg + e1z * xib(j,k)
     &            + e2z * yib(j,k)
     &            + e3z * zib(j,k)
      x0bc(j,kd) = x0b(j)
      y0bc(j,kd) = y0b(j)
      z0bc(j,kd) = z0b(j)

* angle
      e01x(j) = e1x * ei1x(j,k)
     &        + e2x * ei1y(j,k)
     &        + e3x * ei1z(j,k)
      e01y(j) = e1y * ei1x(j,k)
     &        + e2y * ei1y(j,k)
     &        + e3y * ei1z(j,k)
      e01z(j) = e1z * ei1x(j,k)
     &        + e2z * ei1y(j,k)
     &        + e3z * ei1z(j,k)

      e02x(j) = e1x * ei2x(j,k)
     &        + e2x * ei2y(j,k)
     &        + e3x * ei2z(j,k)
      e02y(j) = e1y * ei2x(j,k)
     &        + e2y * ei2y(j,k)
     &        + e3y * ei2z(j,k)
      e02z(j) = e1z * ei2x(j,k)
     &        + e2z * ei2y(j,k)
     &        + e3z * ei2z(j,k)

      e03x(j) = e1x * ei3x(j,k)
     &        + e2x * ei3y(j,k)
     &        + e3x * ei3z(j,k)
      e03y(j) = e1y * ei3x(j,k)
     &        + e2y * ei3y(j,k)
     &        + e3y * ei3z(j,k)
      e03z(j) = e1z * ei3x(j,k)
     &        + e2z * ei3y(j,k)
     &        + e3z * ei3z(j,k)

      e01xc(j,kd) = e01x(j)
      e01yc(j,kd) = e01y(j)
      e01zc(j,kd) = e01z(j)
      e02xc(j,kd) = e02x(j)
      e02yc(j,kd) = e02y(j)
      e02zc(j,kd) = e02z(j)

* velocity
* -Since ver.4.2.0- below part was modifed 
c     u(1) = ome1*e1x + ome2*e2x + ome3*e3x
c     u(2) = ome1*e1y + ome2*e2y + ome3*e3y
c     u(3) = ome1*e1z + ome2*e2z + ome3*e3z
      u(1) = (ome1-omev1(k))*e1x+(ome2-omev2(k))*e2x+(ome3-omev3(k))*e3x
      u(2) = (ome1-omev1(k))*e1y+(ome2-omev2(k))*e2y+(ome3-omev3(k))*e3y
      u(3) = (ome1-omev1(k))*e1z+(ome2-omev2(k))*e2z+(ome3-omev3(k))*e3z

      v(1) = x0a(j) - xg
      v(2) = y0a(j) - yg
      v(3) = z0a(j) - zg

      vx0a(j) = vgx
     &        + e1x * vxia(j,k)
     &        + e2x * vyia(j,k)
     &        + e3x * vzia(j,k)
     &        + exprox(u,v)

      vy0a(j) = vgy
     &        + e1y * vxia(j,k)
     &        + e2y * vyia(j,k)
     &        + e3y * vzia(j,k)
     &        + exproy(u,v)

      vz0a(j) = vgz
     &        + e1z * vxia(j,k)
     &        + e2z * vyia(j,k)
     &        + e3z * vzia(j,k)
     &        + exproz(u,v)

      v(1) = x0b(j) - xg
      v(2) = y0b(j) - yg
      v(3) = z0b(j) - zg

      vx0b(j) = vgx
     &        + e1x * vxib(j,k)
     &        + e2x * vyib(j,k)
     &        + e3x * vzib(j,k)
     &        + exprox(u,v)

      vy0b(j) = vgy
     &        + e1y * vxib(j,k)
     &        + e2y * vyib(j,k)
     &        + e3y * vzib(j,k)
     &        + exproy(u,v)

      vz0b(j) = vgz
     &        + e1z * vxib(j,k)
     &        + e2z * vyib(j,k)
     &        + e3z * vzib(j,k)
     &        + exproz(u,v)
     
* acceleration
* a
      v(1) = vxia(j,k)*e1x + vyia(j,k)*e2x + vzia(j,k)*e3x 
      v(2) = vxia(j,k)*e1y + vyia(j,k)*e2y + vzia(j,k)*e3y 
      v(3) = vxia(j,k)*e1z + vyia(j,k)*e2z + vzia(j,k)*e3z 

      ax0a(j) = 
     &        + e1x * axia(j,k)
     &        + e2x * ayia(j,k)
     &        + e3x * azia(j,k)
     &        + 2.0d0 * exprox(u,v)
      ay0a(j) = 
     &        + e1y * axia(j,k)
     &        + e2y * ayia(j,k)
     &        + e3y * azia(j,k)
     &        + 2.0d0 * exproy(u,v)
      az0a(j) = 
     &        + e1z * axia(j,k)
     &        + e2z * ayia(j,k)
     &        + e3z * azia(j,k)
     &        + 2.0d0 * exproz(u,v)

      v(1) = x0a(j) - xg
      v(2) = y0a(j) - yg
      v(3) = z0a(j) - zg
      omepx = exprox(u,v)
      omepy = exproy(u,v)
      omepz = exproz(u,v)
      v(1) = omepx
      v(2) = omepy
      v(3) = omepz

      ax0a(j) = ax0a(j) + exprox(u,v)
      ay0a(j) = ay0a(j) + exproy(u,v)
      az0a(j) = az0a(j) + exproz(u,v)
     
* b
      v(1) = vxib(j,k)*e1x + vyib(j,k)*e2x + vzib(j,k)*e3x 
      v(2) = vxib(j,k)*e1y + vyib(j,k)*e2y + vzib(j,k)*e3y 
      v(3) = vxib(j,k)*e1z + vyib(j,k)*e2z + vzib(j,k)*e3z 

      ax0b(j) = 
     &        + e1x * axib(j,k)
     &        + e2x * ayib(j,k)
     &        + e3x * azib(j,k)
     &        + 2.0d0 * exprox(u,v)
      ay0b(j) = 
     &        + e1y * axib(j,k)
     &        + e2y * ayib(j,k)
     &        + e3y * azib(j,k)
     &        + 2.0d0 * exproy(u,v)
      az0b(j) = 
     &        + e1z * axib(j,k)
     &        + e2z * ayib(j,k)
     &        + e3z * azib(j,k)
     &        + 2.0d0 * exproz(u,v)

      v(1) = x0b(j) - xg
      v(2) = y0b(j) - yg
      v(3) = z0b(j) - zg
      omepx = exprox(u,v)
      omepy = exproy(u,v)
      omepz = exproz(u,v)
      v(1) = omepx
      v(2) = omepy
      v(3) = omepz

      ax0b(j) = ax0b(j) + exprox(u,v)
      ay0b(j) = ay0b(j) + exproy(u,v)
      az0b(j) = az0b(j) + exproz(u,v)

* Below part was added since version 4.2.0
      if (k.eq.1) then
        omev1km1 = omev1(numl)
        omev2km1 = omev2(numl)
        omev3km1 = omev3(numl)
      else
        omev1km1 = omev1(k-1)
        omev2km1 = omev2(k-1)
        omev3km1 = omev3(k-1)
      endif
      u(1) = (omev1(k+1)-omev1km1)/(2.0d0*dt) *e1x
     &      +(omev2(k+1)-omev2km1)/(2.0d0*dt) *e2x
     &      +(omev3(k+1)-omev3km1)/(2.0d0*dt) *e3x
      u(2) = (omev1(k+1)-omev1km1)/(2.0d0*dt) *e1y
     &      +(omev2(k+1)-omev2km1)/(2.0d0*dt) *e2y
     &      +(omev3(k+1)-omev3km1)/(2.0d0*dt) *e3y
      u(3) = (omev1(k+1)-omev1km1)/(2.0d0*dt) *e1z
     &      +(omev2(k+1)-omev2km1)/(2.0d0*dt) *e2z
     &      +(omev3(k+1)-omev3km1)/(2.0d0*dt) *e3z

      v(1) = x0a(j) - xg
      v(2) = y0a(j) - yg
      v(3) = z0a(j) - zg

      ax0a(j) = ax0a(j) - exprox(u,v)
      ay0a(j) = ay0a(j) - exproy(u,v)
      az0a(j) = az0a(j) - exproz(u,v)

      v(1) = x0b(j) - xg
      v(2) = y0b(j) - yg
      v(3) = z0b(j) - zg

      ax0b(j) = ax0b(j) - exprox(u,v)
      ay0b(j) = ay0b(j) - exproy(u,v)
      az0b(j) = az0b(j) - exproz(u,v)

* obtain mass center fixed coordinate
      x0ad(j) = x0a(j) - xgloop
      y0ad(j) = y0a(j) - ygloop
      z0ad(j) = z0a(j) - zgloop
      x0bd(j) = x0b(j) - xgloop
      y0bd(j) = y0b(j) - ygloop
      z0bd(j) = z0b(j) - zgloop
  300 continue



***************************************************************
***************************************************************
***************************************************************
* Output position, angle, and velocity for multi project mode *
***************************************************************
***************************************************************
***************************************************************
      if (multi.ge.1) then

* positions and velocities of all segments
      open(unit=10,file='Output_data/posvelsh.dat'
     &            ,status='old',access='append')
      open(unit=11,file='Output_data/posvelsc.dat',status='unknown')
  302 format(1p3e18.8e3)
      do j=1,numb

* position
c     write(10,302) x0a(j), y0a(j), z0a(j)
c     write(10,302) x0b(j), y0b(j), z0b(j)
      write(11,302) x0a(j), y0a(j), z0a(j)
      write(11,302) x0b(j), y0b(j), z0b(j)

* angle
c     write(10,302) e01x(j), e01y(j), e01z(j)
c     write(10,302) e02x(j), e02y(j), e02z(j)
c     write(10,302) e03x(j), e03y(j), e03z(j)
      write(11,302) e01x(j), e01y(j), e01z(j)
      write(11,302) e02x(j), e02y(j), e02z(j)
      write(11,302) e03x(j), e03y(j), e03z(j)

* velocity
c     write(10,302) vx0a(j), vy0a(j), vz0a(j)
c     write(10,302) vx0b(j), vy0b(j), vz0b(j)
      write(11,302) vx0a(j), vy0a(j), vz0a(j)
      write(11,302) vx0b(j), vy0b(j), vz0b(j)

      enddo
      close(10)
      close(11)

* position and velocity of center of mass
      open(unit=10,file='Output_data/posvelgh.dat'
     &            ,status='old',access='append')
c     write(10,302) xg, yg, zg
c     write(10,302) e1x, e1y, e1z 
c     write(10,302) e2x, e2y, e2z
c     write(10,302) e3x, e3y, e3z
c     write(10,302) vgx, vgy, vgz
c     write(10,302) ome1, ome2, ome3
      close(10)

      open(unit=11,file='Output_data/posvelgc.dat',status='unknown')
      write(11,302) xg, yg, zg
      write(11,302) e1x, e1y, e1z 
      write(11,302) e2x, e2y, e2z
      write(11,302) e3x, e3y, e3z
      write(11,302) vgx, vgy, vgz
* -Since ver.4.2.0- below part was modifed 
c     write(11,302) ome1, ome2, ome3
      write(11,302) ome1-omev1(k), ome2-omev2(k), ome3-omev3(k)
      close(11)

      open(unit=11,file='Output_data/time_current.dat',status='unknown')
      write(11,*) t
      close(11)

      endif

***************************************
***************************************
***************************************
***************************************
***************************************
* Waiting loop for multi project mode *
***************************************
***************************************
***************************************
***************************************
***************************************
      if (multi.ge.1) then

        goorstop = -2
 1115   continue
        open(unit=10,file='go2.dat',status='replace',err=1115)
 1117   continue
        write(10,*,err=1117) '0 '
        close(10)

      do i=1,4000000

 2220   continue
c       write(*,*) "opening"
        open(unit=10,file='go.dat',action='read'
     &         ,status='unknown',err=2220)
 2222   continue
        close(10)
 2221   continue
        open(unit=10,file='go.dat',action='read'
     &         ,status='unknown',err=2221)
c       write(*,*) "reading"
        read(10,*,end=2222,err=2222) goorstop, dum1, dum2, dum3
        close(10)

c       if (goorstop.eq.1) then
        if ((goorstop.eq.1)
     &      .and.(dum1.eq.loop)
     &      .and.(dum2.eq.((kd-1)/2+1))
     &      .and.(dum3.eq.rkloop)) then
          goto 1111
        endif

 1112   continue
 1113   continue
 1114   continue

        if ((i.ge.500).and.
     &      (mod(i,100).eq.0)) then
          write(*,*) 'waiting iteration=',i, goorstop,dum1,dum2,dum3
          write(*,*) 'waiting iteration=',i, 1,loop,(kd-1)/2+1,rkloop
        endif

*       waiting
        call sleepqq(1)

      enddo

      write(*,*) 'Waiting loop was finished. waited too long.'
      stop

 1111 continue

        goorstop = -2
 1116   continue
        open(unit=10,file='go2.dat',status='replace',err=1116)
 1118   continue
        write(10,*,err=1118) '1 '
        close(10)

      endif




************************
************************
************************
************************
************************
************************
************************
************************
************************
************************
************************
* compute fluid forces *
************************
************************
************************
************************
************************
************************
************************
************************
************************
************************
************************
      areasum = 0.0d0
      fxall=0.0d0

      fxalll=0.0d0
      fxallb =0.0d0
      fxalli =0.0d0
      fxallt =0.0d0
      fxalln =0.0d0

      fyall=0.0d0
      fzall=0.0d0
      m1all=0.0d0

      m1alll=0.0d0
      m1allb=0.0d0
      m1alli=0.0d0
      m1allt=0.0d0
      m1alln=0.0d0

      m2all=0.0d0
      m3all=0.0d0

      powera = 0.0d0
      powerp = 0.0d0

      powtx  = 0.0d0
      pownx  = 0.0d0
      powix  = 0.0d0
      powty  = 0.0d0
      powny  = 0.0d0
      powiy  = 0.0d0
      powtz  = 0.0d0
      pownz  = 0.0d0
      powiz  = 0.0d0

      buoxsum = 0.0d0
      buoysum = 0.0d0
      buozsum = 0.0d0
      buosum  = 0.0d0

      do 390  i=1,6
      do 390 j2=1,6
      inera(i,j2) = 0.0d0
  390 continue

      do 391 jf=1,nf
      fxup(jf)=0.0d0
      fyup(jf)=0.0d0
      fzup(jf)=0.0d0

      fxdw(jf)=0.0d0
      fydw(jf)=0.0d0
      fzdw(jf)=0.0d0
  391 continue

      do 392 jf=1,nm
      mxup(jf)=0.0d0
      myup(jf)=0.0d0
      mzup(jf)=0.0d0

      mxdw(jf)=0.0d0
      mydw(jf)=0.0d0
      mzdw(jf)=0.0d0
 392  continue

      do 400 j=1,numb
      fbx(j)=0.0d0
      fby(j)=0.0d0
      fbz(j)=0.0d0

      fsx(j)=0.0d0
      fsy(j)=0.0d0
      fsz(j)=0.0d0

      ftx(j)=0.0d0
      fty(j)=0.0d0
      ftz(j)=0.0d0

      fnx(j)=0.0d0
      fny(j)=0.0d0
      fnz(j)=0.0d0

      fix(j)=0.0d0
      fiy(j)=0.0d0
      fiz(j)=0.0d0

      ptx(j)=0.0d0
      pty(j)=0.0d0
      ptz(j)=0.0d0

      pnx(j)=0.0d0
      pny(j)=0.0d0
      pnz(j)=0.0d0

      pix(j)=0.0d0
      piy(j)=0.0d0
      piz(j)=0.0d0

      mb1(j)=0.0d0
      mb2(j)=0.0d0
      mb3(j)=0.0d0

      ms1(j)=0.0d0
      ms2(j)=0.0d0
      ms3(j)=0.0d0

      mt1(j)=0.0d0
      mt2(j)=0.0d0
      mt3(j)=0.0d0

      mn1(j)=0.0d0
      mn2(j)=0.0d0
      mn3(j)=0.0d0

      mi1(j)=0.0d0
      mi2(j)=0.0d0
      mi3(j)=0.0d0
     
      dlen = len(j)/dble(div)

**********************************************************
* optain buoyancy and its moment acting on both tip side 
* of truncated elliptic cone
**********************************************************

      do 401 i=1,divr
      do 401 j2=1,divt
* a-side
* obtain coordinate which buoyancy is acting on
      xb = x0a(j)
     &   + e01x(j) * ra1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02x(j) * ra2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
 
      yb = y0a(j)
     &   + e01y(j) * ra1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02y(j) * ra2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

      zb = z0a(j)
     &   + e01z(j) * ra1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02z(j) * ra2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

* if zb<0, buoyancy acts since it is under water surface
      if (zb.le.0.0d0) then

*     obtain area of tiny quadrangle
         area = ra1(j)*(dble(i)-0.5d0)/dble(divr)
     &          * 2.0*pi/dble(divt)
     &          * ra1(j)/dble(divr)
     &          * ra2(j)/ra1(j)

*     obtain unit vector normal to the quadrangle
         nx = (x0b(j)-x0a(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))
         ny = (y0b(j)-y0a(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))
         nz = (z0b(j)-z0a(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))

*     compute buoyancy
         dfbx = + area * (-zb) * grav * nx
         dfby = + area * (-zb) * grav * ny
         dfbz = + area * (-zb) * grav * nz

         fbx(j) = fbx(j) + dfbx
         fby(j) = fby(j) + dfby
         fbz(j) = fbz(j) + dfbz

*     obtain moment due to buoyancy
         u(1) = (xb - xg)
         u(2) = (yb - yg)
         u(3) = (zb - zg)

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

         mtmpx = exprox(u,v)
         mtmpy = exproy(u,v)
         mtmpz = exproz(u,v)

         mb1(j) = mb1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
         mb2(j) = mb2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
         mb3(j) = mb3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

*      for COB calculation
         buoxsum = buoxsum + dfbz * xb
         buoysum = buoysum + dfbz * yb
         buozsum = buozsum + dfbz * zb *0.5d0
         buosum  = buosum  + dfbz

* obtain force and moment between segments
* for calculation of joint torque
      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
         call 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)
      endif
      end if

* b-side
* obtain coordinate which buoyancy is acting on
      xb = x0b(j)
     &   + e01x(j) * rb1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02x(j) * rb2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

      yb = y0b(j)
     &   + e01y(j) * rb1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02y(j) * rb2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

      zb = z0b(j)
     &   + e01z(j) * rb1(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &   + e02z(j) * rb2(j) * (dble(i)-0.5d0)/dble(divr)
     &   * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

* if zb<0, buoyancy acts since it is under water surface
      if (zb.le.0.0d0) then

*     obtain area of tiny quadrangle
         area = rb1(j)*(dble(i)-0.5d0)/dble(divr)
     &          * 2.0*pi/dble(divt)
     &          * rb1(j)/dble(divr)
     &          * rb2(j)/rb1(j)

*     obtain unit vector normal to the quadrangle
         nx = (x0a(j)-x0b(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))
         ny = (y0a(j)-y0b(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))
         nz = (z0a(j)-z0b(j))
     &      / vecabs(x0b(j)-x0a(j),y0b(j)-y0a(j),z0b(j)-z0a(j))

*     compute buoyancy
         dfbx = + area * (-zb) * grav * nx
         dfby = + area * (-zb) * grav * ny
         dfbz = + area * (-zb) * grav * nz

         fbx(j) = fbx(j) + dfbx
         fby(j) = fby(j) + dfby
         fbz(j) = fbz(j) + dfbz

*     obtain moment due to buoyancy
         u(1) = (xb - xg)
         u(2) = (yb - yg)
         u(3) = (zb - zg)

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

         mtmpx = exprox(u,v)
         mtmpy = exproy(u,v)
         mtmpz = exproz(u,v)

         mb1(j) = mb1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
         mb2(j) = mb2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
         mb3(j) = mb3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

*      for COB calculation
         buoxsum = buoxsum + dfbz * xb
         buoysum = buoysum + dfbz * yb
         buozsum = buozsum + dfbz * zb *0.5d0
         buosum  = buosum  + dfbz

* obtain force and moment between segments
* for calculation of joint torque
      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
         call 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)
      endif
      end if

  401 continue

********************************************************************
* optain buoyancy and its moment acting on truncated elliptic cone
********************************************************************
      do 410 i=1,div

* obtain position of the ellipse center
      xc = x0a(j) + (x0b(j)-x0a(j))*(dble(i)-0.5d0)/dble(div)
      yc = y0a(j) + (y0b(j)-y0a(j))*(dble(i)-0.5d0)/dble(div)
      zc = z0a(j) + (z0b(j)-z0a(j))*(dble(i)-0.5d0)/dble(div)

* to view force distribution
      if(rkloop.eq.1) then
      xdist(j,i) = xc
      ydist(j,i) = yc
      zdist(j,i) = zc
      xdistc(j,i,kd) = xc
      ydistc(j,i,kd) = yc
      zdistc(j,i,kd) = zc
      end if

      r1 = ra1(j) + (rb1(j)-ra1(j))*(dble(i)-0.5d0)/dble(div)
      r2 = ra2(j) + (rb2(j)-ra2(j))*(dble(i)-0.5d0)/dble(div)

* obtain velocity of the ellipse center
      vxc = vx0a(j) + (vx0b(j)-vx0a(j))*(dble(i)-0.5d0)/dble(div)
      vyc = vy0a(j) + (vy0b(j)-vy0a(j))*(dble(i)-0.5d0)/dble(div)
      vzc = vz0a(j) + (vz0b(j)-vz0a(j))*(dble(i)-0.5d0)/dble(div)

* obtain acceleration of the ellipse center
      axc = ax0a(j) + (ax0b(j)-ax0a(j))*(dble(i)-0.5d0)/dble(div)
      ayc = ay0a(j) + (ay0b(j)-ay0a(j))*(dble(i)-0.5d0)/dble(div)
      azc = az0a(j) + (az0b(j)-az0a(j))*(dble(i)-0.5d0)/dble(div)

* set dm 0
      dm =0.0d0

* set number of submerged quadrangle 0
      numen = 0.0d0

* set ellipse circumference 0
      enshuu = 0.0d0

* set distribution force 0
      fbxdist(j,i) =0.0d0
      fbydist(j,i) =0.0d0
      fbzdist(j,i) =0.0d0

      do 420 j2=1,divt
* calculate ellipse circumference
      xd1= xc 
     &      + e01x(j) * r1 * dcos(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
     &      + e02x(j) * r2 * dsin(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
      yd1= yc 
     &      + e01y(j) * r1 * dcos(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
     &      + e02y(j) * r2 * dsin(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
      zd1= zc 
     &      + e01z(j) * r1 * dcos(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
     &      + e02z(j) * r2 * dsin(2.0d0*pi*(dble(j2)-1.0d0)/dble(divt))
     
      xd2= xc 
     &      + e01x(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))
     &      + e02x(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))
      yd2= yc 
     &      + e01y(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))
     &      + e02y(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))
      zd2= zc 
     &      + e01z(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))
     &      + e02z(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.0d0)/dble(divt))

      enshuu = enshuu
     &      + dsqrt( (xd1-xd2)**2+(yd1-yd2)**2+(zd1-zd2)**2 )

* obtain coordinate which buoyancy is acting on
      xb = xc 
     &      + e01x(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &      + e02x(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
      yb = yc 
     &      + e01y(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &      + e02y(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
      zb = zc 
     &      + e01z(j) * r1 * dcos(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))
     &      + e02z(j) * r2 * dsin(2.0d0*pi*(dble(j2)-0.5d0)/dble(divt))

* store xb,yb,zb as xbm,ybm,zbm arrays
      xbm(j,i,j2) = xb
      ybm(j,i,j2) = yb
      zbm(j,i,j2) = zb

* if zb<0, buoyancy acts since it is under water surface
      if (zb.le.0.0d0) then

*     obtain each vertex of tiny quadrangle
*     (first 1,2-sides)
          xc12 = x0a(j) + (x0b(j)-x0a(j))*(dble(i-1))/dble(div)
          yc12 = y0a(j) + (y0b(j)-y0a(j))*(dble(i-1))/dble(div)
          zc12 = z0a(j) + (z0b(j)-z0a(j))*(dble(i-1))/dble(div)
          r112 = ra1(j) + (rb1(j)-ra1(j))*(dble(i-1))/dble(div)
          r212 = ra2(j) + (rb2(j)-ra2(j))*(dble(i-1))/dble(div)

          xb1 = xc12
     &         + e01x(j) * r112 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02x(j) * r212 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))
          yb1 = yc12
     &         + e01y(j) * r112 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02y(j) * r212 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))
          zb1 = zc12
     &         + e01z(j) * r112 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02z(j) * r212 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))

          xb2 = xc12
     &         + e01x(j) * r112 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02x(j) * r212 * dsin(2.0d0*pi*(dble(j2))/dble(divt))
          yb2 = yc12
     &         + e01y(j) * r112 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02y(j) * r212 * dsin(2.0d0*pi*(dble(j2))/dble(divt))
          zb2 = zc12
     &         + e01z(j) * r112 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02z(j) * r212 * dsin(2.0d0*pi*(dble(j2))/dble(divt))

*     (next, 3,4-sides)
          xc34 = x0a(j) + (x0b(j)-x0a(j))*(dble(i))/dble(div)
          yc34 = y0a(j) + (y0b(j)-y0a(j))*(dble(i))/dble(div)
          zc34 = z0a(j) + (z0b(j)-z0a(j))*(dble(i))/dble(div)
          r134 = ra1(j) + (rb1(j)-ra1(j))*(dble(i))/dble(div)
          r234 = ra2(j) + (rb2(j)-ra2(j))*(dble(i))/dble(div)

          xb3 = xc34
     &         + e01x(j) * r134 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02x(j) * r234 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))
          yb3 = yc34
     &         + e01y(j) * r134 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02y(j) * r234 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))
          zb3 = zc34
     &         + e01z(j) * r134 * dcos(2.0d0*pi*(dble(j2-1))/dble(divt))
     &         + e02z(j) * r234 * dsin(2.0d0*pi*(dble(j2-1))/dble(divt))

          xb4 = xc34
     &         + e01x(j) * r134 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02x(j) * r234 * dsin(2.0d0*pi*(dble(j2))/dble(divt))
          yb4 = yc34
     &         + e01y(j) * r134 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02y(j) * r234 * dsin(2.0d0*pi*(dble(j2))/dble(divt))
          zb4 = zc34
     &         + e01z(j) * r134 * dcos(2.0d0*pi*(dble(j2))/dble(divt))
     &         + e02z(j) * r234 * dsin(2.0d0*pi*(dble(j2))/dble(divt))

*     obtain area of tiny quadrangle
          s11 = xb2 - xb1
          s12 = yb2 - yb1
          s13 = zb2 - zb1

          s21 = xb3 - xb1
          s22 = yb3 - yb1
          s23 = zb3 - zb1

          s31 = xb3 - xb4
          s32 = yb3 - yb4
          s33 = zb3 - zb4

          s41 = xb2 - xb4
          s42 = yb2 - yb4
          s43 = zb2 - zb4

          area1 = 0.5d0 * dsqrt ( 
     &          (vecabs(s11,s12,s13))**2 
     &        * (vecabs(s21,s22,s23))**2
     &        - (inpro(s11,s12,s13,s21,s22,s23))**2 )
          area2 = 0.5d0 * dsqrt ( 
     &          (vecabs(s31,s32,s33))**2 
     &        * (vecabs(s41,s42,s43))**2
     &        - (inpro(s31,s32,s33,s41,s42,s43))**2 )
          area3 = 0.5d0 * dsqrt ( 
     &          (vecabs(s11,s12,s13))**2 
     &        * (vecabs(s41,s42,s43))**2
     &        - (inpro(-s11,-s12,-s13,-s41,-s42,-s43))**2 )
          area4 = 0.5d0 * dsqrt ( 
     &          (vecabs(s31,s32,s33))**2 
     &        * (vecabs(s21,s22,s23))**2
     &        - (inpro(-s31,-s32,-s33,-s21,-s22,-s23))**2 )

          area = (area1 + area2 + area3 + area4) *0.5d0
c         write(*,*) area
          areasum = areasum + area

*     obtain unit vector normal to the quadrangle
          u(1) = ( (xb4 + xb2)-(xb3 + xb1) )*0.5d0
          u(2) = ( (yb4 + yb2)-(yb3 + yb1) )*0.5d0
          u(3) = ( (zb4 + zb2)-(zb3 + zb1) )*0.5d0

          v(1) = ( (xb4 + xb3)-(xb2 + xb1) )*0.5d0
          v(2) = ( (yb4 + yb3)-(yb2 + yb1) )*0.5d0
          v(3) = ( (zb4 + zb3)-(zb2 + zb1) )*0.5d0

          nx = - exprox(u,v) 
     &                / vecabs(exprox(u,v),exproy(u,v),exproz(u,v))
          ny = - exproy(u,v) 
     &                / vecabs(exprox(u,v),exproy(u,v),exproz(u,v))
          nz = - exproz(u,v) 
     &                / vecabs(exprox(u,v),exproy(u,v),exproz(u,v))

*     compute buoyancy
          dfbx = + area * (-zb) * grav * nx
          dfby = + area * (-zb) * grav * ny
          dfbz = + area * (-zb) * grav * nz

          fbx(j) = fbx(j) + dfbx
          fby(j) = fby(j) + dfby
          fbz(j) = fbz(j) + dfbz

          if(rkloop.eq.4) then
          fbxdist(j,i) = fbxdist(j,i) + dfbx
          fbydist(j,i) = fbydist(j,i) + dfby
          fbzdist(j,i) = fbzdist(j,i) + dfbz
          end if

*     compute dm
          dm = dm + area * vzc * nz

*     obtain moment due to buoyancy
          u(1) = (xb - xg)
          u(2) = (yb - yg)
          u(3) = (zb - zg)

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

          mtmpx = exprox(u,v)
          mtmpy = exproy(u,v)
          mtmpz = exproz(u,v)

          mb1(j) = mb1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
          mb2(j) = mb2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
          mb3(j) = mb3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

*      for COB calculation
         buoxsum = buoxsum + dfbz * xb
         buoysum = buoysum + dfbz * yb
         buozsum = buozsum + dfbz * zb *0.5d0
         buosum  = buosum  + dfbz

* obtain force and moment between segments
* for calculation of joint torque
      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
          call 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)
      endif

* add 1 into number of submerged quadrangle
          numen = numen + 1.0d0
      end if
      
  420 continue

* save numen
      numen2(j,i)=numen

********************************************************
* obtain force due to momentum change at water surface *
********************************************************

* obtain vector from tip-a to tip-b
      xba = x0b(j)-x0a(j)
      yba = y0b(j)-y0a(j)
      zba = z0b(j)-z0a(j)

* make a vector which has horizontal components of above vector
* and make a arbitrary vertical vector 
      xbah = xba
      ybah = yba
      zbah = 0.0d0

      xbav = 0.0d0
      ybav = 0.0d0
      zbav = -1.0d0

* normalize vector from tip-a to tip-b
      baabs = vecabs(xba,yba,zba)
      xba = xba/baabs
      yba = yba/baabs
      zba = zba/baabs
      
* obtain component parallel to ellipse axis 
* using inner product against the normalized vector
      inproh = inpro(xbah,ybah,zbah,xba,yba,zba)
      xbah2 = inproh * xba
      ybah2 = inproh * yba
      zbah2 = inproh * zba
      inprov = inpro(xbav,ybav,zbav,xba,yba,zba)
      xbav2 = inprov * xba
      ybav2 = inprov * yba
      zbav2 = inprov * zba

* calculate component normal to ellipse axis
* by decucting component parallel to ellipse axis
      xbah3 = xbah - xbah2
      ybah3 = ybah - ybah2
      zbah3 = zbah - zbah2

      xbav3 = xbav - xbav2
      ybav3 = ybav - ybav2
      zbav3 = zbav - zbav2

* calculate each absolute value, and take larger one
      h3abs = vecabs(xbah3,ybah3,zbah3)
      v3abs = vecabs(xbav3,ybav3,zbav3)

      if (h3abs.gt.v3abs) then
          xba4 = xbah3
          yba4 = ybah3
          zba4 = zbah3
      else
          xba4 = xbav3
          yba4 = ybav3
          zba4 = zbav3
      end if

* normalize
      ba4abs = vecabs(xba4,yba4,zba4)
      xba5 = xba4/ba4abs
      yba5 = yba4/ba4abs
      zba5 = zba4/ba4abs

* take inner product against velocity vector
      vinpro = inpro(vxc,vyc,vzc,xba5,yba5,zba5)

* obtain vector which is normal to axis and does not have
* horizontal component by multiplying value of the inner product
      vxc2 = vinpro * xba5
      vyc2 = vinpro * yba5
      vzc2 = vinpro * zba5

* obtain fluid force
c     dfsx = - dabs(dm) * vxc2 * cm
c     dfsy = - dabs(dm) * vyc2 * cm
c     dfsz = - dabs(dm) * vzc2 * cm
      dfsx = 0.0d0
      dfsy = 0.0d0
      dfsz = 0.0d0

      fsx(j) = fsx(j) + dfsx
      fsy(j) = fsy(j) + dfsy
      fsz(j) = fsz(j) + dfsz

      fsxdist(j,i) = dfsx
      fsydist(j,i) = dfsy
      fszdist(j,i) = dfsz

* obtain moment due to fluid force
      u(1) = (xc - xg)
      u(2) = (yc - yg)
      u(3) = (zc - zg)

      v(1) = dfsx
      v(2) = dfsy
      v(3) = dfsz

      mtmpx = exprox(u,v)
      mtmpy = exproy(u,v)
      mtmpz = exproz(u,v)

      ms1(j) = ms1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      ms2(j) = ms2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      ms3(j) = ms3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

* obtain force and moment between segments
* for calculation of joint torque
      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
      call bodyforce(dfsx,dfsy,dfsz,fxup,fyup,fzup
     &              ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &              ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &              ,xg,yg,zg)
      endif

*********************************************************
* force proportional to square of velocity (tangential)
*********************************************************

* optain tangential component of velocity
* using inner product against unit vector from tip-a to tip-b
      inproh = inpro(vxc,vyc,vzc,xba,yba,zba)
      vxct = inproh * xba
      vyct = inproh * yba
      vzct = inproh * zba

* obtain magnitude
      vctabs = dsqrt(vxct**2+vyct**2+vzct**2)

* obtain fluid force
      dftx = -0.5d0 *dlen*enshuu *cs(j) *vctabs *(numen/dble(divt))*vxct
      dfty = -0.5d0 *dlen*enshuu *cs(j) *vctabs *(numen/dble(divt))*vyct
      dftz = -0.5d0 *dlen*enshuu *cs(j) *vctabs *(numen/dble(divt))*vzct

      ftx(j) = ftx(j) + dftx
      fty(j) = fty(j) + dfty
      ftz(j) = ftz(j) + dftz

* calculate power
      ptx(j) = ptx(j) + dftx * vxc
      pty(j) = pty(j) + dfty * vyc
      ptz(j) = ptz(j) + dftz * vzc

      ftxdist(j,i) = dftx
      ftydist(j,i) = dfty
      ftzdist(j,i) = dftz

* obtain moment due to fluid force
      u(1) = (xc - xg)
      u(2) = (yc - yg)
      u(3) = (zc - zg)

      v(1) = dftx
      v(2) = dfty
      v(3) = dftz

      mtmpx = exprox(u,v)
      mtmpy = exproy(u,v)
      mtmpz = exproz(u,v)

      mt1(j) = mt1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      mt2(j) = mt2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      mt3(j) = mt3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
* obtain force and moment between segments
* for calculation of joint torque
      call bodyforce(dftx,dfty,dftz,fxup,fyup,fzup
     &              ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &              ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &              ,xg,yg,zg)
      endif
*****************************************************
* force proportional to square of velocity (normal)
*****************************************************

* optain normal component by deducting tangential component
      vxcn = vxc - vxct
      vycn = vyc - vyct
      vzcn = vzc - vzct

* make magnitude square by multiplying magnitude itself
      vcnabs = dsqrt(vxcn**2+vycn**2+vzcn**2)
      vxcns = vxcn * vcnabs
      vycns = vycn * vcnabs
      vzcns = vzcn * vcnabs

* obtain two components parallel to ellipse axes
* by inner product
* axis-1
      inpron1 = inpro(vxcns,vycns,vzcns,e01x(j),e01y(j),e01z(j))
      vxcn1 = inpron1 * e01x(j)
      vycn1 = inpron1 * e01y(j)
      vzcn1 = inpron1 * e01z(j)
* axis-2
      inpron2 = inpro(vxcns,vycns,vzcns,e02x(j),e02y(j),e02z(j))
      vxcn2 = inpron2 * e02x(j)
      vycn2 = inpron2 * e02y(j)
      vzcn2 = inpron2 * e02z(j)

* obtain fluid force
      dfnx =
     &   - 0.5d0 *dlen*2.0d0*r2 *(r2/r1) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vxcn1
     &   - 0.5d0 *dlen*2.0d0*r1 *(r1/r2) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vxcn2
      dfny =
     &   - 0.5d0 *dlen*2.0d0*r2 *(r2/r1) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vycn1
     &   - 0.5d0 *dlen*2.0d0*r1 *(r1/r2) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vycn2
      dfnz =
     &   - 0.5d0 *dlen*2.0d0*r2 *(r2/r1) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vzcn1
     &   - 0.5d0 *dlen*2.0d0*r1 *(r1/r2) ** pcn
     &           *cn(j)*(numen/dble(divt)) * vzcn2

      fnx(j) = fnx(j) + dfnx
      fny(j) = fny(j) + dfny
      fnz(j) = fnz(j) + dfnz

* calculate power
      pnx(j) = pnx(j) + dfnx * vxc
      pny(j) = pny(j) + dfny * vyc
      pnz(j) = pnz(j) + dfnz * vzc

      fnxdist(j,i) = dfnx
      fnydist(j,i) = dfny
      fnzdist(j,i) = dfnz

* obtain moment due to fluid force
      u(1) = (xc - xg)
      u(2) = (yc - yg)
      u(3) = (zc - zg)

      v(1) = dfnx
      v(2) = dfny
      v(3) = dfnz

      mtmpx = exprox(u,v)
      mtmpy = exproy(u,v)
      mtmpz = exproz(u,v)

      mn1(j) = mn1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      mn2(j) = mn2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      mn3(j) = mn3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
* obtain force and moment between segments
* for calculation of joint torque
      call bodyforce(dfnx,dfny,dfnz,fxup,fyup,fzup
     &              ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &              ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &              ,xg,yg,zg)
      endif
**************************************************
* inertial force due to added mass 
* (component which can be calculated in advance)
**************************************************

* obtain tangential component of acceleration
* using inner product against unit vector from tip-a to tip-b
      inproh = inpro(axc,ayc,azc,xba,yba,zba)
      axct = inproh * xba
      ayct = inproh * yba
      azct = inproh * zba

* obtain normal component by deducting tangential component
      axcn = axc - axct
      aycn = ayc - ayct
      azcn = azc - azct

* obtain two components parallel to ellipse axes
* by inner product
* axis-1
      inproi1 = inpro(axcn,aycn,azcn,e01x(j),e01y(j),e01z(j))
      axcn1 = inproi1 * e01x(j)
      aycn1 = inproi1 * e01y(j)
      azcn1 = inproi1 * e01z(j)
* axis-2
      inproi2 = inpro(axcn,aycn,azcn,e02x(j),e02y(j),e02z(j))
      axcn2 = inproi2 * e02x(j)
      aycn2 = inproi2 * e02y(j)
      azcn2 = inproi2 * e02z(j)

* obtain fluid force
      dfix =
     &      - dlen*pi*r2**2 * (numen/dble(divt)) * axcn1 * ca(j)
     &      - dlen*pi*r1**2 * (numen/dble(divt)) * axcn2 * ca(j)
      dfiy =
     &      - dlen*pi*r2**2 * (numen/dble(divt)) * aycn1 * ca(j)
     &      - dlen*pi*r1**2 * (numen/dble(divt)) * aycn2 * ca(j)
      dfiz =
     &      - dlen*pi*r2**2 * (numen/dble(divt)) * azcn1 * ca(j)
     &      - dlen*pi*r1**2 * (numen/dble(divt)) * azcn2 * ca(j)

      fix(j) = fix(j) + dfix
      fiy(j) = fiy(j) + dfiy
      fiz(j) = fiz(j) + dfiz

* calculate power
      pix(j) = pix(j) + dfix * vxc
      piy(j) = piy(j) + dfiy * vyc
      piz(j) = piz(j) + dfiz * vzc

      fixdist(j,i) = dfix
      fiydist(j,i) = dfiy
      fizdist(j,i) = dfiz

* obtain moment due to fluid force
      u(1) = (xc - xg)
      u(2) = (yc - yg)
      u(3) = (zc - zg)

      v(1) = dfix
      v(2) = dfiy
      v(3) = dfiz

      mtmpx = exprox(u,v)
      mtmpy = exproy(u,v)
      mtmpz = exproz(u,v)

      mi1(j) = mi1(j) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      mi2(j) = mi2(j) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      mi3(j) = mi3(j) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

* obtain force and moment between segments
* for calculation of joint torque
c           call bodyforce(dfix,dfiy,dfiz,fxup,fyup,fzup
c    &                    ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
c    &                    ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j)

************************************
* inertial force due to added mass 
* (rigid body component)
************************************
      do 450 ij=1,6
      axc = 0.0d0
      ayc = 0.0d0
      azc = 0.0d0
      if(ij.eq.1) then
          axc = 1.0d0
      else if(ij.eq.2) then
          ayc = 1.0d0
      else if(ij.eq.3) then
          azc = 1.0d0
      else if(ij.eq.4) then
          u(1) = e1x
          u(2) = e1y
          u(3) = e1z
          v(1) = xc - xg
          v(2) = yc - yg
          v(3) = zc - zg
          axc = exprox(u,v)
          ayc = exproy(u,v)
          azc = exproz(u,v)
      else if(ij.eq.5) then
          u(1) = e2x
          u(2) = e2y
          u(3) = e2z
          v(1) = xc - xg
          v(2) = yc - yg
          v(3) = zc - zg
          axc = exprox(u,v)
          ayc = exproy(u,v)
          azc = exproz(u,v)
      else if(ij.eq.6) then
          u(1) = e3x
          u(2) = e3y
          u(3) = e3z
          v(1) = xc - xg
          v(2) = yc - yg
          v(3) = zc - zg
          axc = exprox(u,v)
          ayc = exproy(u,v)
          azc = exproz(u,v)
      else
          write(*,*) 'No way! 6'
          stop
      end if

* obtain tangential component of accecleration
* by inner product against unit vector from tip-a to tip-b
      inproh = inpro(axc,ayc,azc,xba,yba,zba)
      axct = inproh * xba
      ayct = inproh * yba
      azct = inproh * zba

* obtain normal component by deducting tangential component
      axcn = axc - axct
      aycn = ayc - ayct
      azcn = azc - azct

* obtain two components parallel to ellipse axes
* by inner product
* axis-1
      inproi1 = inpro(axcn,aycn,azcn,e01x(j),e01y(j),e01z(j))
      axcn1 = inproi1 * e01x(j)
      aycn1 = inproi1 * e01y(j)
      azcn1 = inproi1 * e01z(j)
* axis-2
      inproi2 = inpro(axcn,aycn,azcn,e02x(j),e02y(j),e02z(j))
      axcn2 = inproi2 * e02x(j)
      aycn2 = inproi2 * e02y(j)
      azcn2 = inproi2 * e02z(j)

* obtain fluid force
      dfix =
     &      + dlen*pi*r2**2 * (numen/dble(divt)) * axcn1 * ca(j)
     &      + dlen*pi*r1**2 * (numen/dble(divt)) * axcn2 * ca(j)
      dfiy =
     &      + dlen*pi*r2**2 * (numen/dble(divt)) * aycn1 * ca(j)
     &      + dlen*pi*r1**2 * (numen/dble(divt)) * aycn2 * ca(j)
      dfiz =
     &      + dlen*pi*r2**2 * (numen/dble(divt)) * azcn1 * ca(j)
     &      + dlen*pi*r1**2 * (numen/dble(divt)) * azcn2 * ca(j)

      inera(1,ij) = inera(1,ij) + dfix
      inera(2,ij) = inera(2,ij) + dfiy
      inera(3,ij) = inera(3,ij) + dfiz

      ineraorg(1,ij) = inera(1,ij)
      ineraorg(2,ij) = inera(2,ij)
      ineraorg(3,ij) = inera(3,ij)

      inerad(1,ij,j,i) = dfix
      inerad(2,ij,j,i) = dfiy
      inerad(3,ij,j,i) = dfiz

* obtain moment due to fluid force
      u(1) = (xc - xg)
      u(2) = (yc - yg)
      u(3) = (zc - zg)

      v(1) = dfix
      v(2) = dfiy
      v(3) = dfiz

      mtmpx = exprox(u,v)
      mtmpy = exproy(u,v)
      mtmpz = exproz(u,v)

      inera(4,ij) = inera(4,ij) + inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      inera(5,ij) = inera(5,ij) + inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      inera(6,ij) = inera(6,ij) + inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

      ineraorg(4,ij) = inera(4,ij)
      ineraorg(5,ij) = inera(5,ij)
      ineraorg(6,ij) = inera(6,ij)

      inerad(4,ij,j,i) = inpro(mtmpx,mtmpy,mtmpz,e1x,e1y,e1z)
      inerad(5,ij,j,i) = inpro(mtmpx,mtmpy,mtmpz,e2x,e2y,e2z)
      inerad(6,ij,j,i) = inpro(mtmpx,mtmpy,mtmpz,e3x,e3y,e3z)

  450 continue

      if (rkloop.eq.1) then
      bou = 0.2d0
      fxdist(j,i) = 0.0d0
c    &              + fbxdist(j,i)/len(j)*bou
     &              + fsxdist(j,i)/len(j)*bou
     &              + ftxdist(j,i)/len(j)*bou
     &              + fnxdist(j,i)/len(j)*bou
     &              + fixdist(j,i)/len(j)*bou

      fydist(j,i) = 0.0d0
c    &              + fbydist(j,i)/len(j)*bou
     &              + fsydist(j,i)/len(j)*bou
     &              + ftydist(j,i)/len(j)*bou
     &              + fnydist(j,i)/len(j)*bou
     &              + fiydist(j,i)/len(j)*bou

      fzdist(j,i) = 0.0d0
c    &              + fbzdist(j,i)/len(j)*bou
     &              + fszdist(j,i)/len(j)*bou
     &              + ftzdist(j,i)/len(j)*bou
     &              + fnzdist(j,i)/len(j)*bou
     &              + fizdist(j,i)/len(j)*bou

      fxdistc(j,i,kd) = fxdist(j,i)
      fydistc(j,i,kd) = fydist(j,i)
      fzdistc(j,i,kd) = fzdist(j,i)
      end if

****************************************************************************
****************************
* judge of back direction  *
****************************

       bd = (x0a(10)-x0a(11))*(y0a(7)-y0a(11))
     &      -(y0a(10)-y0a(11))*(x0a(7)-x0a(11))

      if(bd.ge.0.0d0) then
      jbd = 1.0d0

      else 
      jbd = -1.0d0

      end if

*******************************************
*TimeStep for BodyMovement and FluidForce *
*******************************************
      tany = 0.0d0 + (t-(lastloop-1.0d0))*ttt

      if(tany.le.0.001d0) then
      tany = 0.0d0
      end if

*********************************
*compute fluid force for Anybody*
*********************************

c      fxany(j,i) = 0.0d0
c      fyany(j,i) = 0.0d0
c      fzany(j,i) = 0.0d0


      fxany(j,i) = 0.0d0 - (fbzdist(j,i)+fszdist(j,i)+ftzdist(j,i)
     &    +fnzdist(j,i)+fizdist(j,i))
     &    *bu(8)*bu(8)*bu(8)*bu(8)*1000.0d0/ttt/ttt
     &    *jbd
c    &    *0.0001d0
c     &    *998.0d0*height**4/cycle**2
     

      fyany(j,i) = 0.0d0 - (fbxdist(j,i)+fsxdist(j,i)+ftxdist(j,i)
     &    +fnxdist(j,i)+fixdist(j,i))
     &    *bu(8)*bu(8)*bu(8)*bu(8)*1000.0d0/ttt/ttt
c    &    *0.0001d0
c     &    *998.0d0*height**4/cycle**2
     

      fzany(j,i) = 0.0d0 + (fbydist(j,i)+fsydist(j,i)+ftydist(j,i)
     &    +fnydist(j,i)+fiydist(j,i))
     &    *bu(8)*bu(8)*bu(8)*bu(8)*1000.0d0/ttt/ttt
     &    *jbd
c    &    *0.0001d0
c     &    *998.0d0*height**4/cycle**2
     

************************************
*compute body movement for Anybody *
************************************

      rtan1 =   0.0d0 + (z0a(10)-z0a(11))
     &        /(dsqrt((y0a(10)-y0a(11))*(y0a(10)-y0a(11))
     &        +(x0a(10)-x0a(11))*(x0a(10)-x0a(11))))
      role =    0.0d0 - datan(rtan1) * 180.0d0/pi * jbd

      rtan2 =   0.0d0 + (z0b(1)-z0a(1))/len(1)
      pitch =   0.0d0 + dasin(rtan2) * 180.0d0/pi * jbd

      rtan3 =   0.0d0 + (y0b(1)-y0a(1))/(x0b(1)-x0a(1))
      yaw =     0.0d0 - datan(rtan3) * 180.0d0/pi * jbd


      if(tany.eq.0.0d0) then
      role0 = role
      pitch0= pitch
      yaw0  = yaw
      end if


      if(tany.eq.0.0d0) then
      xori=x0a(1)
      yori=y0a(1)
      end if

      waistx = -z0a(1)*bu(8)*jbd
      waisty = -((x0a(1)-xori) - xst*tany/ttt)*bu(8)
      waistz = ((y0a(1)-yori) - yst*tany/ttt)*bu(8)*jbd


      if(tany.eq.0.0d0) then
      wstx0=waistx
      wsty0=waisty
      wstz0=waistz
      end if


************************
* gravity for Anybody  *
************************

      gravX = 9.8d0*jbd

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


**************************
* compute comsumed power *
**************************
      if (rkloop.eq.1) then

* total power
      powera = powera
     &       + ftxdist(j,i) * vxc
     &       + fnxdist(j,i) * vxc
     &       + fixdist(j,i) * vxc

     &       + ftydist(j,i) * vyc
     &       + fnydist(j,i) * vyc
     &       + fiydist(j,i) * vyc

     &       + ftzdist(j,i) * vzc
     &       + fnzdist(j,i) * vzc
     &       + fizdist(j,i) * vzc

      powtx = powtx + ftxdist(j,i) * vxc
      pownx = pownx + fnxdist(j,i) * vxc
      powix = powix + fixdist(j,i) * vxc

      powty = powty + ftydist(j,i) * vyc
      powny = powny + fnydist(j,i) * vyc
      powiy = powiy + fiydist(j,i) * vyc

      powtz = powtz + ftzdist(j,i) * vzc
      pownz = pownz + fnzdist(j,i) * vzc
      powiz = powiz + fizdist(j,i) * vzc

* efficient power used as propulsion
      powerp = powerp
     &       + ftxdist(j,i) * vxc

      end if
  410 continue

********************************************
* sum up all fluid force and its moment
********************************************
      fxall = fxall
     &          + fbx(j)
     &          + fsx(j)
     &          + ftx(j)
     &          + fnx(j)
     &          + fix(j)

      fxalll = fxalll + fbx(j)
     &                + fsx(j)
     &                + ftx(j)
     &                + fnx(j)
     &                + fix(j)
      fxallb = fxallb + fbx(j)
      fxalli = fxalli + fix(j)
      fxallt = fxallt + ftx(j)
      fxalln = fxalln + fnx(j)

      fyall = fyall
     &          + fby(j)
     &          + fsy(j)
     &          + fty(j)
     &          + fny(j)
     &          + fiy(j)
      fzall = fzall
     &          + fbz(j)
     &          + fsz(j)
     &          + ftz(j)
     &          + fnz(j)
     &          + fiz(j)

      m1all = m1all
     &          + mb1(j)
     &          + ms1(j)
     &          + mt1(j)
     &          + mn1(j)
     &          + mi1(j)

      m1alll = m1alll
     &          + mb1(j)
c    &          + ms1(j)
     &          + mt1(j)
     &          + mn1(j)
     &          + mi1(j)
      m1allb = m1allb + mb1(j)
      m1alli = m1alli + mi1(j)
      m1allt = m1allt + mt1(j)
      m1alln = m1alln + mn1(j)

      m2all = m2all
     &          + mb2(j)
     &          + ms2(j)
     &          + mt2(j)
     &          + mn2(j)
     &          + mi2(j)
      m3all = m3all
     &          + mb3(j)
     &          + ms3(j)
     &          + mt3(j)
     &          + mn3(j)
     &          + mi3(j)
  400 continue

**************************************
* calculate COB (center of buoyancy) *
**************************************
      if (rkloop.eq.1) then
      xcob(kd) = buoxsum / buosum
      ycob(kd) = buoysum / buosum
      zcob(kd) = buozsum / buosum
      endif

*********************************************************
* calculate average of comsumed power during one cycle
*********************************************************
      if (rkloop.eq.1) then
      poweraav = poweraav + powera/(dble(numl)/2.0d0)
      powerpav = powerpav + powerp/(dble(numl)/2.0d0)

      powtxav  = powtxav  + powtx /(dble(numl)/2.0d0)
      pownxav  = pownxav  + pownx /(dble(numl)/2.0d0)
      powixav  = powixav  + powix /(dble(numl)/2.0d0)

      powtyav  = powtyav  + powty /(dble(numl)/2.0d0)
      pownyav  = pownyav  + powny /(dble(numl)/2.0d0)
      powiyav  = powiyav  + powiy /(dble(numl)/2.0d0)

      powtzav  = powtzav  + powtz /(dble(numl)/2.0d0)
      pownzav  = pownzav  + pownz /(dble(numl)/2.0d0)
      powizav  = powizav  + powiz /(dble(numl)/2.0d0)

      end if


**********************
* compute roll angle *
**********************

      if (rkloop.eq.1) then
* calculte component about principal axis (axis-1)
* by time-integrating ome1
* -Since ver.4.2.0- below part was modifed 
c     roll1 = roll1 + 2.0*dt * ome1
      roll1 = roll1 + 2.0*dt * (ome1-omev1(k))

* calculate compoent about trunk axis
* from direction of vector

* first, normalize trunk axis vector just in case
      e01abs = vecabs(e01x(1),e01y(1),e01z(1))
      e01x1 = e01x(1)/e01abs
      e01y1 = e01y(1)/e01abs
      e01z1 = e01z(1)/e01abs

      e02abs = vecabs(e02x(1),e02y(1),e02z(1))
      e02x1 = e02x(1)/e02abs
      e02y1 = e02y(1)/e02abs
      e02z1 = e02z(1)/e02abs

      e03abs = vecabs(e03x(1),e03y(1),e03z(1))
      e03x1 = e03x(1)/e03abs
      e03y1 = e03y(1)/e03abs
      e03z1 = e03z(1)/e03abs

c     rollb = dasin( e01z1/dcos(dasin(e03z1)) )
      rollb = dasin( e02z1/dcos(dasin(e03z1)) )
      end if


***********************************************************
* Output force, moment ,power ,fluidforce and bodymovement*
***********************************************************
      if (rkloop.eq.1) then

c     if((flagoutf.eq.2).or.
c    &  ((flagoutf.eq.1).and.(loop.eq.lastloop))) then
c     call outforce(numb,t,fix,fnx,ftx,fbx
c    &                    ,fiy,fny,fty,fby
c    &                    ,fiz,fnz,ftz,fbz)
c     call outmomen(numb,t,mi1,mn1,mt1,mb1
c    &                    ,mi2,mn2,mt2,mb2
c    &                    ,mi3,mn3,mt3,mb3)
c     call outpower(numb,t,pix,pnx,ptx
c    &                    ,piy,pny,pty
c    &                    ,piz,pnz,ptz)
c     endif


************************************************************************
      if ((flagoutm.eq.1).and.(loop.eq.lastloop)) then
************************************************************************

**********************
* Output fluid force *
**********************

 9001 format(1p6e18.8e3)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,1),fyany(1,1),fzany(1,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_1-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,2),fyany(1,2),fzany(1,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,3),fyany(1,3),fzany(1,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,4),fyany(1,4),fzany(1,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,5),fyany(1,5),fzany(1,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,6),fyany(1,6),fzany(1,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,7),fyany(1,7),fzany(1,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,8),fyany(1,8),fzany(1,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,9),fyany(1,9),fzany(1,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_1-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(1,10),fyany(1,10),fzany(1,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,1),fyany(2,1),fzany(2,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,2),fyany(2,2),fzany(2,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,3),fyany(2,3),fzany(2,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,4),fyany(2,4),fzany(2,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,5),fyany(2,5),fzany(2,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,6),fyany(2,6),fzany(2,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,7),fyany(2,7),fzany(2,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,8),fyany(2,8),fzany(2,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,9),fyany(2,9),fzany(2,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_2-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(2,10),fyany(2,10),fzany(2,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,1),fyany(3,1),fzany(3,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,2),fyany(3,2),fzany(3,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,3),fyany(3,3),fzany(3,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,4),fyany(3,4),fzany(3,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,5),fyany(3,5),fzany(3,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,6),fyany(3,6),fzany(3,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,7),fyany(3,7),fzany(3,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,8),fyany(3,8),fzany(3,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,9),fyany(3,9),fzany(3,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_3-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(3,10),fyany(3,10),fzany(3,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,1),fyany(4,1),fzany(4,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,2),fyany(4,2),fzany(4,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,3),fyany(4,3),fzany(4,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,4),fyany(4,4),fzany(4,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,5),fyany(4,5),fzany(4,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,6),fyany(4,6),fzany(4,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,7),fyany(4,7),fzany(4,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,8),fyany(4,8),fzany(4,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,9),fyany(4,9),fzany(4,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_4-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(4,10),fyany(4,10),fzany(4,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,1),fyany(5,1),fzany(5,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,2),fyany(5,2),fzany(5,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,3),fyany(5,3),fzany(5,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,4),fyany(5,4),fzany(5,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,5),fyany(5,5),fzany(5,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,6),fyany(5,6),fzany(5,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,7),fyany(5,7),fzany(5,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,8),fyany(5,8),fzany(5,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,9),fyany(5,9),fzany(5,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_5-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(5,10),fyany(5,10),fzany(5,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,1),fyany(6,1),fzany(6,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,2),fyany(6,2),fzany(6,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,3),fyany(6,3),fzany(6,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,4),fyany(6,4),fzany(6,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,5),fyany(6,5),fzany(6,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,6),fyany(6,6),fzany(6,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,7),fyany(6,7),fzany(6,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,8),fyany(6,8),fzany(6,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,9),fyany(6,9),fzany(6,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_6-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(6,10),fyany(6,10),fzany(6,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,1),fyany(7,1),fzany(7,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,2),fyany(7,2),fzany(7,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,3),fyany(7,3),fzany(7,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,4),fyany(7,4),fzany(7,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,5),fyany(7,5),fzany(7,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,6),fyany(7,6),fzany(7,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,7),fyany(7,7),fzany(7,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,8),fyany(7,8),fzany(7,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,9),fyany(7,9),fzany(7,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_7-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(7,10),fyany(7,10),fzany(7,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,1),fyany(8,1),fzany(8,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,2),fyany(8,2),fzany(8,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,3),fyany(8,3),fzany(8,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,4),fyany(8,4),fzany(8,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,5),fyany(8,5),fzany(8,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,6),fyany(8,6),fzany(8,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,7),fyany(8,7),fzany(8,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,8),fyany(8,8),fzany(8,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,9),fyany(8,9),fzany(8,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_8-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(8,10),fyany(8,10),fzany(8,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,1),fyany(9,1),fzany(9,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,2),fyany(9,2),fzany(9,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,3),fyany(9,3),fzany(9,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,4),fyany(9,4),fzany(9,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,5),fyany(9,5),fzany(9,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,6),fyany(9,6),fzany(9,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,7),fyany(9,7),fzany(9,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,8),fyany(9,8),fzany(9,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,9),fyany(9,9),fzany(9,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_9-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(9,10),fyany(9,10),fzany(9,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,1),fyany(10,1),fzany(10,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_10-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,2),fyany(10,2),fzany(10,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,3),fyany(10,3),fzany(10,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,4),fyany(10,4),fzany(10,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,5),fyany(10,5),fzany(10,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,6),fyany(10,6),fzany(10,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,7),fyany(10,7),fzany(10,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,8),fyany(10,8),fzany(10,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,9),fyany(10,9),fzany(10,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_10-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(10,10),fyany(10,10),fzany(10,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,1),fyany(11,1),fzany(11,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_11-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,2),fyany(11,2),fzany(11,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,3),fyany(11,3),fzany(11,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,4),fyany(11,4),fzany(11,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,5),fyany(11,5),fzany(11,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,6),fyany(11,6),fzany(11,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,7),fyany(11,7),fzany(11,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,8),fyany(11,8),fzany(11,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,9),fyany(11,9),fzany(11,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_11-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(11,10),fyany(11,10),fzany(11,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,1),fyany(12,1),fzany(12,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_12-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,2),fyany(12,2),fzany(12,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,3),fyany(12,3),fzany(12,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,4),fyany(12,4),fzany(12,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,5),fyany(12,5),fzany(12,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,6),fyany(12,6),fzany(12,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,7),fyany(12,7),fzany(12,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,8),fyany(12,8),fzany(12,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,9),fyany(12,9),fzany(12,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_12-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(12,10),fyany(12,10),fzany(12,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,1),fyany(13,1),fzany(13,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_13-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,2),fyany(13,2),fzany(13,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,3),fyany(13,3),fzany(13,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,4),fyany(13,4),fzany(13,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,5),fyany(13,5),fzany(13,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,6),fyany(13,6),fzany(13,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,7),fyany(13,7),fzany(13,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,8),fyany(13,8),fzany(13,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,9),fyany(13,9),fzany(13,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_13-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(13,10),fyany(13,10),fzany(13,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,1),fyany(14,1),fzany(14,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_14-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,2),fyany(14,2),fzany(14,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,3),fyany(14,3),fzany(14,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,4),fyany(14,4),fzany(14,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,5),fyany(14,5),fzany(14,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,6),fyany(14,6),fzany(14,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,7),fyany(14,7),fzany(14,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,8),fyany(14,8),fzany(14,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,9),fyany(14,9),fzany(14,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_14-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(14,10),fyany(14,10),fzany(14,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,1),fyany(15,1),fzany(15,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_15-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,2),fyany(15,2),fzany(15,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,3),fyany(15,3),fzany(15,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,4),fyany(15,4),fzany(15,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,5),fyany(15,5),fzany(15,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,6),fyany(15,6),fzany(15,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,7),fyany(15,7),fzany(15,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,8),fyany(15,8),fzany(15,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,9),fyany(15,9),fzany(15,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_15-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(15,10),fyany(15,10),fzany(15,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,1),fyany(16,1),fzany(16,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_16-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,2),fyany(16,2),fzany(16,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,3),fyany(16,3),fzany(16,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,4),fyany(16,4),fzany(16,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,5),fyany(16,5),fzany(16,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,6),fyany(16,6),fzany(16,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,7),fyany(16,7),fzany(16,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,8),fyany(16,8),fzany(16,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,9),fyany(16,9),fzany(16,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_16-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(16,10),fyany(16,10),fzany(16,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,1),fyany(17,1),fzany(17,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_17-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,2),fyany(17,2),fzany(17,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,3),fyany(17,3),fzany(17,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,4),fyany(17,4),fzany(17,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,5),fyany(17,5),fzany(17,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,6),fyany(17,6),fzany(17,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,7),fyany(17,7),fzany(17,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,8),fyany(17,8),fzany(17,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,9),fyany(17,9),fzany(17,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_17-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(17,10),fyany(17,10),fzany(17,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,1),fyany(18,1),fzany(18,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_18-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,2),fyany(18,2),fzany(18,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,3),fyany(18,3),fzany(18,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,4),fyany(18,4),fzany(18,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,5),fyany(18,5),fzany(18,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,6),fyany(18,6),fzany(18,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,7),fyany(18,7),fzany(18,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,8),fyany(18,8),fzany(18,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,9),fyany(18,9),fzany(18,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_18-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(18,10),fyany(18,10),fzany(18,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,1),fyany(19,1),fzany(19,1)
      close(10)

       open(unit=10,file='AnyBody/FluidForce/AnyFF_19-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,2),fyany(19,2),fzany(19,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,3),fyany(19,3),fzany(19,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,4),fyany(19,4),fzany(19,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,5),fyany(19,5),fzany(19,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,6),fyany(19,6),fzany(19,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,7),fyany(19,7),fzany(19,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,8),fyany(19,8),fzany(19,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,9),fyany(19,9),fzany(19,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_19-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(19,10),fyany(19,10),fzany(19,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,1),fyany(20,1),fzany(20,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,2),fyany(20,2),fzany(20,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,3),fyany(20,3),fzany(20,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,4),fyany(20,4),fzany(20,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,5),fyany(20,5),fzany(20,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,6),fyany(20,6),fzany(20,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,7),fyany(20,7),fzany(20,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,8),fyany(20,8),fzany(20,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,9),fyany(20,9),fzany(20,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_20-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(20,10),fyany(20,10),fzany(20,10)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-1.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,1),fyany(21,1),fzany(21,1)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-2.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,2),fyany(21,2),fzany(21,2)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-3.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,3),fyany(21,3),fzany(21,3)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-4.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,4),fyany(21,4),fzany(21,4)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-5.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,5),fyany(21,5),fzany(21,5)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-6.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,6),fyany(21,6),fzany(21,6)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-7.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,7),fyany(21,7),fzany(21,7)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-8.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,8),fyany(21,8),fzany(21,8)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-9.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,9),fyany(21,9),fzany(21,9)
      close(10)

      open(unit=10,file='AnyBody/FluidForce/AnyFF_21-10.txt'
     &          ,status='unknown',access='append')
      write(10,9001) tany,fxany(21,10),fyany(21,10),fzany(21,10)
      close(10)

*********************
* Output Time Step  *
*********************

      open(unit=10,file='AnyBody/timeStep.any'
     &          ,status='unknown',access='append')

      write(10,*) tany,','

      if (kd.eq.numl-1) then
      write(10,*) ttt
      end if

      close(10)

***********************
* Output bodymovement *
***********************

      open(unit=10,file='AnyBody/BodyMovement/role.any'
     &          ,status='unknown',access='append')
      write(10,*) role,','

      if (kd.eq.numl-1) then
      write(10,*) role0
      end if

      close(10)


      open(unit=10,file='AnyBody/BodyMovement/pitch.any'
     &          ,status='unknown',access='append')
      write(10,*) pitch,','

      if (kd.eq.numl-1) then
      write(10,*) pitch0
      end if

      close(10)


      open(unit=10,file='AnyBody/BodyMovement/yaw.any'
     &          ,status='unknown',access='append')
      write(10,*) yaw,','

      if (kd.eq.numl-1) then
      write(10,*) yaw0
      end if

      close(10)


      open(unit=10,file='AnyBody/BodyMovement/ComX.any'
     &               ,status='unknown',access='append')
      write(10,*) waistx,','

      if (kd.eq.numl-1) then
      write(10,*) wstx0
      end if

      close(10)


      open(unit=10,file='AnyBody/BodyMovement/ComY.any'
     &               ,status='unknown',access='append')
      write(10,*) waisty,','

      if (kd.eq.numl-1) then
      write(10,*) waisty
      end if

      close(10)


      open(unit=10,file='AnyBody/BodyMovement/ComZ.any'
     &               ,status='unknown',access='append')
      write(10,*) waistz,','

      if (kd.eq.numl-1) then
      write(10,*) wstz0
      end if

      close(10)


******************
* Output gravity *
******************

      if (kd.eq.1) then
      open(unit=10,file='AnyBody/Gravity.any'
     &               ,status='unknown',access='append')
      write(10,*) '{',gravX,',0.0,0.0}'
      close(10)
      end if

************************************************************************
      end if
************************************************************************

      end if

**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
* From equation of motion, obtain velocity and angular velocity
* at next time step (Runge-Kutta method)
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************
**********************************************************************

**********************************************************************
* first, add mass and moment of inertia itself 
* to added mass matrix
**********************************************************************
      inera(1,1) = inera(1,1) + mall
      inera(2,2) = inera(2,2) + mall
      inera(3,3) = inera(3,3) + mall
      inera(4,4) = inera(4,4) + iall1(k)
      inera(5,5) = inera(5,5) + iall2(k)
      inera(6,6) = inera(6,6) + iall3(k)

*******************************************
* load external force from other projects *
*******************************************
      fxex = 0.0d0
      fyex = 0.0d0
      fzex = 0.0d0
      m1ex = 0.0d0
      m2ex = 0.0d0
      m3ex = 0.0d0
      if (multi.ge.1) then
      open(unit=10,file='external_force.dat',action='read'
     &         ,status='unknown',err=7772)
      read(10,*,end=7772,err=7772) fxex
      read(10,*,end=7772,err=7772) fyex
      read(10,*,end=7772,err=7772) fzex
      read(10,*,end=7772,err=7772) mxex
      read(10,*,end=7772,err=7772) myex
      read(10,*,end=7772,err=7772) mzex
      close(10)
      goto 7771
 7772 continue
* if failed, all values are set as zero
      fxex = 0.0d0
      fyex = 0.0d0
      fzex = 0.0d0
      mxex = 0.0d0
      myex = 0.0d0
      mzex = 0.0d0
 7771 continue

* Compute moments about 1,2,3 axes from moments about x,y,z axes
      m1ex = inpro(mxex,myex,mzex,e1x,e1y,e1z)
      m2ex = inpro(mxex,myex,mzex,e2x,e2y,e2z)
      m3ex = inpro(mxex,myex,mzex,e3x,e3y,e3z)

      write(*,7773) fxex,fyex,fzex,m1ex,m2ex,m3ex
 7773 format('    ',1p6e13.3e3)
      end if


************************************************
* put force and moment as vector form in order
************************************************
      fv(1) = fxall
     &              + fxex
      fv(2) = fyall
     &              + fyex
      fv(3) = fzall
     &         -mall * grav
     &              + fzex
      fv(4) = 
     &              +(iall2(k)-iall3(k))*ome2*ome3 
     &              + m1all
     &              - idall1(k) * ome1
     &              + m1ex
      fv(5) = 
     &              +(iall3(k)-iall1(k))*ome3*ome1 
     &              + m2all
     &              - idall2(k) * ome2
     &              + m2ex
      fv(6) = 
     &              +(iall1(k)-iall2(k))*ome1*ome2
     &              + m3all
     &              - idall3(k) * ome3
     &              + m3ex
     
      fvorg(1) = fv(1)
      fvorg(2) = fv(2)
      fvorg(3) = fv(3)
      fvorg(4) = fv(4)
      fvorg(5) = fv(5)
      fvorg(6) = fv(6)

*****************************************************
* solve simultaneous equation using SLATEC library
* (returned fv is the solution)
*****************************************************
      call dgeco(inera,6,6,ipvt,rcond,zmat)
      call dgesl(inera,6,6,ipvt,fv,0)

      fvd(1,rkloop)=fv(1)
      fvd(2,rkloop)=fv(2)
      fvd(3,rkloop)=fv(3)
      fvd(4,rkloop)=fv(4)
      fvd(5,rkloop)=fv(5)
      fvd(6,rkloop)=fv(6)

* make vector of principal axes in moving coordinate
* for calculation of rotation
      e11 = 1.0d0
      e12 = 0.0d0
      e13 = 0.0d0

      e21 = 0.0d0
      e22 = 1.0d0
      e23 = 0.0d0

      e31 = 0.0d0
      e32 = 0.0d0
      e33 = 1.0d0

**********************
* Runge-Kutta loop 1 *
**********************
      if (rkloop.eq.1) then
*     first, save present step values
          xgorg = xg
          ygorg = yg
          zgorg = zg
          xgorgd(kd) = xg
          ygorgd(kd) = yg
          zgorgd(kd) = zg

          e1xorg = e1x
          e1yorg = e1y
          e1zorg = e1z
          e2xorg = e2x
          e2yorg = e2y
          e2zorg = e2z
          e3xorg = e3x
          e3yorg = e3y
          e3zorg = e3z

          e1xorgd(kd) = e1x
          e1yorgd(kd) = e1y
          e1zorgd(kd) = e1z
          e2xorgd(kd) = e2x
          e2yorgd(kd) = e2y
          e2zorgd(kd) = e2z
          e3xorgd(kd) = e3x
          e3yorgd(kd) = e3y
          e3zorgd(kd) = e3z

          vgxorg = vgx
          vgyorg = vgy
          vgzorg = vgz

          ome1org = ome1
          ome2org = ome2
          ome3org = ome3

*     save velocity
          vgx1= vgx
          vgy1= vgy
          vgz1= vgz
          ome11 = ome1
          ome21 = ome2
          ome31 = ome3

*     change present value for the case rkloop=2
*     (translational displacement)
          xg = xgorg + dt*vgx
          yg = ygorg + dt*vgy
          zg = zgorg + dt*vgz

*     (rotational displacement)
c         e11new = -ome3*dt *e12 +ome2*dt *e13  + e11
c         e12new =  ome3*dt *e11 -ome1*dt *e13  + e12
c         e13new = -ome2*dt *e11 +ome1*dt *e12  + e13
c         call vecnor(e11new,e12new,e13new)
c      
c         e21new = -ome3*dt *e22 +ome2*dt *e23  + e21
c         e22new =  ome3*dt *e21 -ome1*dt *e23  + e22
c         e23new = -ome2*dt *e21 +ome1*dt *e22  + e23
c         call vecnor(e21new,e22new,e23new)
c      
c         e31new = -ome3*dt *e32 +ome2*dt *e33  + e31
c         e32new =  ome3*dt *e31 -ome1*dt *e33  + e32
c         e33new = -ome2*dt *e31 +ome1*dt *e32  + e33
c         call vecnor(e31new,e32new,e33new)
c    
          e11new = -(ome3-omev3(k))*dt*e12+(ome2-omev2(k))*dt*e13 + e11
          e12new =  (ome3-omev3(k))*dt*e11-(ome1-omev1(k))*dt*e13 + e12
          e13new = -(ome2-omev2(k))*dt*e11+(ome1-omev1(k))*dt*e12 + e13
          call vecnor(e11new,e12new,e13new)
       
          e21new = -(ome3-omev3(k))*dt*e22+(ome2-omev2(k))*dt*e23 + e21
          e22new =  (ome3-omev3(k))*dt*e21-(ome1-omev1(k))*dt*e23 + e22
          e23new = -(ome2-omev2(k))*dt*e21+(ome1-omev1(k))*dt*e22 + e23
          call vecnor(e21new,e22new,e23new)
       
          e31new = -(ome3-omev3(k))*dt*e32+(ome2-omev2(k))*dt*e33 + e31
          e32new =  (ome3-omev3(k))*dt*e31-(ome1-omev1(k))*dt*e33 + e32
          e33new = -(ome2-omev2(k))*dt*e31+(ome1-omev1(k))*dt*e32 + e33
          call vecnor(e31new,e32new,e33new)
     
          e1x = e11new * e1xorg + e12new * e2xorg + e13new * e3xorg
          e1y = e11new * e1yorg + e12new * e2yorg + e13new * e3yorg
          e1z = e11new * e1zorg + e12new * e2zorg + e13new * e3zorg

          e2x = e21new * e1xorg + e22new * e2xorg + e23new * e3xorg
          e2y = e21new * e1yorg + e22new * e2yorg + e23new * e3yorg
          e2z = e21new * e1zorg + e22new * e2zorg + e23new * e3zorg

          e3x = e31new * e1xorg + e32new * e2xorg + e33new * e3xorg
          e3y = e31new * e1yorg + e32new * e2yorg + e33new * e3yorg
          e3z = e31new * e1zorg + e32new * e2zorg + e33new * e3zorg

*     (translational velocity)
          vgx = vgxorg + dt * fv(1)
          vgy = vgyorg + dt * fv(2)
          vgz = vgzorg + dt * fv(3)

          if(vgxsolv.eq.0) then
              vgx = vgxg(kd)
          end if
          if(vgysolv.eq.0) then
              vgy = vgyg(kd)
          end if
          if(vgzsolv.eq.0) then
              vgz = vgzg(kd)
          end if

          vgx2= vgx
          vgy2= vgy
          vgz2= vgz

*     (rotational velocity)
          ome1 = ome1org + dt * fv(4)
          ome2 = ome2org + dt * fv(5)
          ome3 = ome3org + dt * fv(6)

          if(ome1solv.eq.0) then
              ome1 = ome1g(kd)
          end if
          if(ome2solv.eq.0) then
              ome2 = ome2g(kd)
          end if
          if(ome3solv.eq.0) then
              ome3 = ome3g(kd)
          end if

          ome12 = ome1
          ome22 = ome2
          ome32 = ome3

*******************************************************************
* Compute inertial force due to added mass (rigid body component) *
*******************************************************************
            fam1 = - (fvorg(1) - mall*fv(1))
            fam2 = - (fvorg(2) - mall*fv(2))
            fam3 = - (fvorg(3) - mall*fv(3))
            mam1 = - (fvorg(4) - iall1(k)*fv(4))
            mam2 = - (fvorg(5) - iall2(k)*fv(5))
            mam3 = - (fvorg(6) - iall3(k)*fv(6))

          if(vgxsolv.eq.0) then
            if(kd.eq.numl-1) then
                vgxtmp = vgxg(1)
            else
                vgxtmp = vgxg(kd+2)
            end if
            famd1 =    -(
     &               ineraorg(1,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(1,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(1,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(1,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(1,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(1,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if

          if(vgysolv.eq.0) then
            if(kd.eq.numl-1) then
                vgytmp = vgyg(1)
            else
                vgytmp = vgyg(kd+2)
            end if
            famd2 =    -(
     &               ineraorg(2,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(2,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(2,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(2,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(2,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(2,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if

          if(vgzsolv.eq.0) then
            if(kd.eq.numl-1) then
                vgztmp = vgzg(1)
            else
                vgztmp = vgzg(kd+2)
            end if
            famd3 =    -(
     &               ineraorg(3,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(3,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(3,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(3,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(3,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(3,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if

          if(ome1solv.eq.0) then
            if(kd.eq.numl-1) then
                ome1tmp = ome1g(1)
            else
                ome1tmp = ome1g(kd+2)
            end if
            mamd1 =    -(
     &               ineraorg(4,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(4,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(4,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(4,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(4,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(4,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if

          if(ome2solv.eq.0) then
            if(kd.eq.numl-1) then
                ome2tmp = ome2g(1)
            else
                ome2tmp = ome2g(kd+2)
            end if
            mamd2 =    -(
     &               ineraorg(5,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(5,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(5,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(5,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(5,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(5,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if

          if(ome3solv.eq.0) then
            if(kd.eq.numl-1) then
                ome3tmp = ome3g(1)
            else
                ome3tmp = ome3g(kd+2)
            end if
            mamd3 =    -(
     &               ineraorg(6,1) * (vgxtmp - vgx)/(2.0d0*dt)
     &             + ineraorg(6,2) * (vgytmp - vgy)/(2.0d0*dt)
     &             + ineraorg(6,3) * (vgztmp - vgz)/(2.0d0*dt)
     &             + ineraorg(6,4) * (ome1tmp-ome1)/(2.0d0*dt)
     &             + ineraorg(6,5) * (ome2tmp-ome2)/(2.0d0*dt)
     &             + ineraorg(6,6) * (ome3tmp-ome3)/(2.0d0*dt)
     &                  )
          end if


***********************************************************
* Output force, moment ,power ,fluidforce and bodymovement*
***********************************************************
c     if (rkloop.eq.1) then

      if((flagoutf.eq.2).or.
     &  ((flagoutf.eq.1).and.(loop.eq.lastloop))) then
      call outforce(numb,t,fix,fnx,ftx,fbx
     &                    ,fiy,fny,fty,fby
     &                    ,fiz,fnz,ftz,fbz,fam1,fam2,fam3
     &                    ,famd1,famd2,famd3)
      call outmomen(numb,t,mi1,mn1,mt1,mb1
     &                    ,mi2,mn2,mt2,mb2
     &                    ,mi3,mn3,mt3,mb3,mam1,mam2,mam3
     &                    ,mamd1,mamd2,mamd3)
      call outpower(numb,t,pix,pnx,ptx
     &                    ,piy,pny,pty
     &                    ,piz,pnz,ptz)
      endif

**********************
* Runge-Kutta loop 2 *
**********************
      else if (rkloop.eq.2) then

*     change present value for the case rkloop=3
*     (translational displacement)
          xg = xgorg + dt*vgx
          yg = ygorg + dt*vgy
          zg = zgorg + dt*vgz

*     (rotational displacement)
c         e11new = -ome3*dt *e12 +ome2*dt *e13  + e11
c         e12new =  ome3*dt *e11 -ome1*dt *e13  + e12
c         e13new = -ome2*dt *e11 +ome1*dt *e12  + e13
c         call vecnor(e11new,e12new,e13new)
c      
c         e21new = -ome3*dt *e22 +ome2*dt *e23  + e21
c         e22new =  ome3*dt *e21 -ome1*dt *e23  + e22
c         e23new = -ome2*dt *e21 +ome1*dt *e22  + e23
c         call vecnor(e21new,e22new,e23new)
c      
c         e31new = -ome3*dt *e32 +ome2*dt *e33  + e31
c         e32new =  ome3*dt *e31 -ome1*dt *e33  + e32
c         e33new = -ome2*dt *e31 +ome1*dt *e32  + e33
c         call vecnor(e31new,e32new,e33new)
     
          e11new = -(ome3-omev3(k))*dt*e12+(ome2-omev2(k))*dt*e13 + e11
          e12new =  (ome3-omev3(k))*dt*e11-(ome1-omev1(k))*dt*e13 + e12
          e13new = -(ome2-omev2(k))*dt*e11+(ome1-omev1(k))*dt*e12 + e13
          call vecnor(e11new,e12new,e13new)
       
          e21new = -(ome3-omev3(k))*dt*e22+(ome2-omev2(k))*dt*e23 + e21
          e22new =  (ome3-omev3(k))*dt*e21-(ome1-omev1(k))*dt*e23 + e22
          e23new = -(ome2-omev2(k))*dt*e21+(ome1-omev1(k))*dt*e22 + e23
          call vecnor(e21new,e22new,e23new)
       
          e31new = -(ome3-omev3(k))*dt*e32+(ome2-omev2(k))*dt*e33 + e31
          e32new =  (ome3-omev3(k))*dt*e31-(ome1-omev1(k))*dt*e33 + e32
          e33new = -(ome2-omev2(k))*dt*e31+(ome1-omev1(k))*dt*e32 + e33
          call vecnor(e31new,e32new,e33new)
     
          e1x = e11new * e1xorg + e12new * e2xorg + e13new * e3xorg
          e1y = e11new * e1yorg + e12new * e2yorg + e13new * e3yorg
          e1z = e11new * e1zorg + e12new * e2zorg + e13new * e3zorg

          e2x = e21new * e1xorg + e22new * e2xorg + e23new * e3xorg
          e2y = e21new * e1yorg + e22new * e2yorg + e23new * e3yorg
          e2z = e21new * e1zorg + e22new * e2zorg + e23new * e3zorg

          e3x = e31new * e1xorg + e32new * e2xorg + e33new * e3xorg
          e3y = e31new * e1yorg + e32new * e2yorg + e33new * e3yorg
          e3z = e31new * e1zorg + e32new * e2zorg + e33new * e3zorg

*     (translational velocity)
          vgx = vgxorg + dt * fv(1)
          vgy = vgyorg + dt * fv(2)
          vgz = vgzorg + dt * fv(3)

          if(vgxsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgx = vgxg(1)
              else
                  vgx = vgxg(kd+2)
              end if
          end if
          if(vgysolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgy = vgyg(1)
              else
                  vgy = vgyg(kd+2)
              end if
          end if
          if(vgzsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgz = vgzg(1)
              else
                  vgz = vgzg(kd+2)
              end if
          end if

          vgx3= vgx
          vgy3= vgy
          vgz3= vgz

*     (rotational velocity)
          ome1 = ome1org + dt * fv(4)
          ome2 = ome2org + dt * fv(5)
          ome3 = ome3org + dt * fv(6)

          if(ome1solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome1 = ome1g(1)
              else
                  ome1 = ome1g(kd+2)
              end if
          end if
          if(ome2solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome2 = ome2g(1)
              else
                  ome2 = ome2g(kd+2)
              end if
          end if
          if(ome3solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome3 = ome3g(1)
              else
                  ome3 = ome3g(kd+2)
              end if
          end if

          ome13 = ome1
          ome23 = ome2
          ome33 = ome3

**********************
* Runge-Kutta loop 3 *
**********************
      else if (rkloop.eq.3) then
*     change present value for the case rkloop=3
*     (translational displacement)
          xg = xgorg + 2.0d0*dt*vgx
          yg = ygorg + 2.0d0*dt*vgy
          zg = zgorg + 2.0d0*dt*vgz

*     (rotational displacement)
c         e11new = -ome3*2.0*dt *e12 +ome2*2.0*dt *e13  + e11
c         e12new =  ome3*2.0*dt *e11 -ome1*2.0*dt *e13  + e12
c         e13new = -ome2*2.0*dt *e11 +ome1*2.0*dt *e12  + e13
c         call vecnor(e11new,e12new,e13new)
c      
c         e21new = -ome3*2.0*dt *e22 +ome2*2.0*dt *e23  + e21
c         e22new =  ome3*2.0*dt *e21 -ome1*2.0*dt *e23  + e22
c         e23new = -ome2*2.0*dt *e21 +ome1*2.0*dt *e22  + e23
c         call vecnor(e21new,e22new,e23new)
c      
c         e31new = -ome3*2.0*dt *e32 +ome2*2.0*dt *e33  + e31
c         e32new =  ome3*2.0*dt *e31 -ome1*2.0*dt *e33  + e32
c         e33new = -ome2*2.0*dt *e31 +ome1*2.0*dt *e32  + e33
c         call vecnor(e31new,e32new,e33new)
     
          e11new = -(ome3-omev3(k))*2.0*dt*e12
     &             +(ome2-omev2(k))*2.0*dt*e13 + e11
          e12new =  (ome3-omev3(k))*2.0*dt*e11
     &             -(ome1-omev1(k))*2.0*dt*e13 + e12
          e13new = -(ome2-omev2(k))*2.0*dt*e11
     &             +(ome1-omev1(k))*2.0*dt*e12 + e13
          call vecnor(e11new,e12new,e13new)
       
          e21new = -(ome3-omev3(k))*2.0*dt*e22
     &             +(ome2-omev2(k))*2.0*dt*e23 + e21
          e22new =  (ome3-omev3(k))*2.0*dt*e21
     &             -(ome1-omev1(k))*2.0*dt*e23 + e22
          e23new = -(ome2-omev2(k))*2.0*dt*e21
     &             +(ome1-omev1(k))*2.0*dt*e22 + e23
          call vecnor(e21new,e22new,e23new)
       
          e31new = -(ome3-omev3(k))*2.0*dt*e32
     &             +(ome2-omev2(k))*2.0*dt*e33 + e31
          e32new =  (ome3-omev3(k))*2.0*dt*e31
     &             -(ome1-omev1(k))*2.0*dt*e33 + e32
          e33new = -(ome2-omev2(k))*2.0*dt*e31
     &             +(ome1-omev1(k))*2.0*dt*e32 + e33
          call vecnor(e31new,e32new,e33new)
     
          e1x = e11new * e1xorg + e12new * e2xorg + e13new * e3xorg
          e1y = e11new * e1yorg + e12new * e2yorg + e13new * e3yorg
          e1z = e11new * e1zorg + e12new * e2zorg + e13new * e3zorg

          e2x = e21new * e1xorg + e22new * e2xorg + e23new * e3xorg
          e2y = e21new * e1yorg + e22new * e2yorg + e23new * e3yorg
          e2z = e21new * e1zorg + e22new * e2zorg + e23new * e3zorg

          e3x = e31new * e1xorg + e32new * e2xorg + e33new * e3xorg
          e3y = e31new * e1yorg + e32new * e2yorg + e33new * e3yorg
          e3z = e31new * e1zorg + e32new * e2zorg + e33new * e3zorg

*     (translational velocity)
          vgx = vgxorg + 2.0* dt * fv(1)
          vgy = vgyorg + 2.0* dt * fv(2)
          vgz = vgzorg + 2.0* dt * fv(3)

          if(vgxsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgx = vgxg(1)
              else
                  vgx = vgxg(kd+2)
              end if
          end if
          if(vgysolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgy = vgyg(1)
              else
                  vgy = vgyg(kd+2)
              end if
          end if
          if(vgzsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgz = vgzg(1)
              else
                  vgz = vgzg(kd+2)
              end if
          end if

          vgx4= vgx
          vgy4= vgy
          vgz4= vgz

*     (rotational velocity)
          ome1 = ome1org + 2.0* dt * fv(4)
          ome2 = ome2org + 2.0* dt * fv(5)
          ome3 = ome3org + 2.0* dt * fv(6)

          if(ome1solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome1 = ome1g(1)
              else
                  ome1 = ome1g(kd+2)
              end if
          end if
          if(ome2solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome2 = ome2g(1)
              else
                  ome2 = ome2g(kd+2)
              end if
          end if
          if(ome3solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome3 = ome3g(1)
              else
                  ome3 = ome3g(kd+2)
              end if
          end if

          ome14 = ome1
          ome24 = ome2
          ome34 = ome3

**********************
* Runge-Kutta loop 4 *
**********************
      else if (rkloop.eq.4) then
*     calculate new time step values using values so far
*     (translational displacement)
          xgnew = xgorg + 2.0d0*dt*(vgx1/6.0d0 + vgx2/3.0d0
     &                            + vgx3/3.0d0 + vgx4/6.0d0)
          ygnew = ygorg + 2.0d0*dt*(vgy1/6.0d0 + vgy2/3.0d0
     &                            + vgy3/3.0d0 + vgy4/6.0d0)
          zgnew = zgorg + 2.0d0*dt*(vgz1/6.0d0 + vgz2/3.0d0
     &                            + vgz3/3.0d0 + vgz4/6.0d0)

*     (rotational displacement)
          ome1 = ome11/6.0d0 + ome12/3.0d0 + ome13/3.0d0 + ome14/6.0d0
          ome2 = ome21/6.0d0 + ome22/3.0d0 + ome23/3.0d0 + ome24/6.0d0
          ome3 = ome31/6.0d0 + ome32/3.0d0 + ome33/3.0d0 + ome34/6.0d0

c         e11new = -ome3*2.0*dt *e12 +ome2*2.0*dt *e13  + e11
c         e12new =  ome3*2.0*dt *e11 -ome1*2.0*dt *e13  + e12
c         e13new = -ome2*2.0*dt *e11 +ome1*2.0*dt *e12  + e13
c         call vecnor(e11new,e12new,e13new)
c      
c         e21new = -ome3*2.0*dt *e22 +ome2*2.0*dt *e23  + e21
c         e22new =  ome3*2.0*dt *e21 -ome1*2.0*dt *e23  + e22
c         e23new = -ome2*2.0*dt *e21 +ome1*2.0*dt *e22  + e23
c         call vecnor(e21new,e22new,e23new)
c      
c         e31new = -ome3*2.0*dt *e32 +ome2*2.0*dt *e33  + e31
c         e32new =  ome3*2.0*dt *e31 -ome1*2.0*dt *e33  + e32
c         e33new = -ome2*2.0*dt *e31 +ome1*2.0*dt *e32  + e33
c         call vecnor(e31new,e32new,e33new)
     
          e11new = -(ome3-omev3(k))*2.0*dt*e12
     &             +(ome2-omev2(k))*2.0*dt*e13 + e11
          e12new =  (ome3-omev3(k))*2.0*dt*e11
     &             -(ome1-omev1(k))*2.0*dt*e13 + e12
          e13new = -(ome2-omev2(k))*2.0*dt*e11
     &             +(ome1-omev1(k))*2.0*dt*e12 + e13
          call vecnor(e11new,e12new,e13new)
       
          e21new = -(ome3-omev3(k))*2.0*dt*e22
     &             +(ome2-omev2(k))*2.0*dt*e23 + e21
          e22new =  (ome3-omev3(k))*2.0*dt*e21
     &             -(ome1-omev1(k))*2.0*dt*e23 + e22
          e23new = -(ome2-omev2(k))*2.0*dt*e21
     &             +(ome1-omev1(k))*2.0*dt*e22 + e23
          call vecnor(e21new,e22new,e23new)
       
          e31new = -(ome3-omev3(k))*2.0*dt*e32
     &             +(ome2-omev2(k))*2.0*dt*e33 + e31
          e32new =  (ome3-omev3(k))*2.0*dt*e31
     &             -(ome1-omev1(k))*2.0*dt*e33 + e32
          e33new = -(ome2-omev2(k))*2.0*dt*e31
     &             +(ome1-omev1(k))*2.0*dt*e32 + e33
          call vecnor(e31new,e32new,e33new)
     
          e1xnew = e11new * e1xorg + e12new * e2xorg + e13new * e3xorg
          e1ynew = e11new * e1yorg + e12new * e2yorg + e13new * e3yorg
          e1znew = e11new * e1zorg + e12new * e2zorg + e13new * e3zorg

          e2xnew = e21new * e1xorg + e22new * e2xorg + e23new * e3xorg
          e2ynew = e21new * e1yorg + e22new * e2yorg + e23new * e3yorg
          e2znew = e21new * e1zorg + e22new * e2zorg + e23new * e3zorg

          e3xnew = e31new * e1xorg + e32new * e2xorg + e33new * e3xorg
          e3ynew = e31new * e1yorg + e32new * e2yorg + e33new * e3yorg
          e3znew = e31new * e1zorg + e32new * e2zorg + e33new * e3zorg

*     orthonormalize obtained principal axes of inertia
          call vecnoror(e1xnew,e1ynew,e1znew
     &                 ,e2xnew,e2ynew,e2znew
     &                 ,e3xnew,e3ynew,e3znew)

*     (translational velocity)
          vgxnew = vgxorg   + 2.0*dt * (fvd(1,1)/6.0d0 + fvd(1,2)/3.0d0
     &                                + fvd(1,3)/3.0d0 + fvd(1,4)/6.0d0)
          vgynew = vgyorg   + 2.0*dt * (fvd(2,1)/6.0d0 + fvd(2,2)/3.0d0
     &                                + fvd(2,3)/3.0d0 + fvd(2,4)/6.0d0)
          vgznew = vgzorg   + 2.0*dt * (fvd(3,1)/6.0d0 + fvd(3,2)/3.0d0
     &                                + fvd(3,3)/3.0d0 + fvd(3,4)/6.0d0)

          if(vgxsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgxnew = vgxg(1)
              else
                  vgxnew = vgxg(kd+2)
              end if
          end if
          if(vgysolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgynew = vgyg(1)
              else
                  vgynew = vgyg(kd+2)
              end if
          end if
          if(vgzsolv.eq.0) then
              if(kd.eq.numl-1) then
                  vgznew = vgzg(1)
              else
                  vgznew = vgzg(kd+2)
              end if
          end if

*     (rotational velocity)
          ome1new = ome1org + 2.0*dt * (fvd(4,1)/6.0d0 + fvd(4,2)/3.0d0
     &                                + fvd(4,3)/3.0d0 + fvd(4,4)/6.0d0)
          ome2new = ome2org + 2.0*dt * (fvd(5,1)/6.0d0 + fvd(5,2)/3.0d0
     &                                + fvd(5,3)/3.0d0 + fvd(5,4)/6.0d0)
          ome3new = ome3org + 2.0*dt * (fvd(6,1)/6.0d0 + fvd(6,2)/3.0d0
     &                                + fvd(6,3)/3.0d0 + fvd(6,4)/6.0d0)

          if(ome1solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome1new = ome1g(1)
              else
                  ome1new = ome1g(kd+2)
              end if
          end if
          if(ome2solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome2new = ome2g(1)
              else
                  ome2new = ome2g(kd+2)
              end if
          end if
          if(ome3solv.eq.0) then
              if(kd.eq.numl-1) then
                  ome3new = ome3g(1)
              else
                  ome3new = ome3g(kd+2)
              end if
          end if

      if (multi.ge.1) then
      write(*,8888) t, xgnew, ygnew, zgnew
 8888 format('t=',1pe16.6e3,'  (xg,yg,zg)=',1p3e16.6e3)
      end if

******************************************************************
* Compute inertial forces due to added mass (rigid body segment) *
* and add them to distributed forces for animation               *
******************************************************************
      do j = 1,numb
      do i = 1,div
        ineradx = 0.0d0
        inerady = 0.0d0
        ineradz = 0.0d0
        do ij=1,6
          acc = fv(ij)
          if((ij.eq.1).and.(vgxsolv.eq.0)) then
            acc = (vgxnew - vgxorg)/(2.0d0*dt)
          endif
          if((ij.eq.2).and.(vgysolv.eq.0)) then
            acc = (vgynew - vgyorg)/(2.0d0*dt)
          endif
          if((ij.eq.3).and.(vgzsolv.eq.0)) then
            acc = (vgznew - vgzorg)/(2.0d0*dt)
          endif
          if((ij.eq.4).and.(ome1solv.eq.0)) then
            acc = (ome1new - ome1org)/(2.0d0*dt)
          endif
          if((ij.eq.5).and.(ome2solv.eq.0)) then
            acc = (ome2new - ome2org)/(2.0d0*dt)
          endif
          if((ij.eq.6).and.(ome3solv.eq.0)) then
            acc = (ome3new - ome3org)/(2.0d0*dt)
          endif
          ineradx = ineradx + inerad(1,ij,j,i)*acc
          inerady = inerady + inerad(2,ij,j,i)*acc 
          ineradz = ineradz + inerad(3,ij,j,i)*acc 
        enddo
          fxdistc(j,i,kd) = fxdistc(j,i,kd) - ineradx/len(j)*bou
          fydistc(j,i,kd) = fydistc(j,i,kd) - inerady/len(j)*bou
          fzdistc(j,i,kd) = fzdistc(j,i,kd) - ineradz/len(j)*bou
      enddo
      enddo


*****************
* time marching *
*****************
      vgx = vgxnew
      vgy = vgynew
      vgz = vgznew

      ome1 = ome1new
      ome2 = ome2new
      ome3 = ome3new

      xg = xgnew
      yg = ygnew
      zg = zgnew

      e1x = e1xnew
      e1y = e1ynew
      e1z = e1znew

      e2x = e2xnew
      e2y = e2ynew
      e2z = e2znew

      e3x = e3xnew
      e3y = e3ynew
      e3z = e3znew

*     in the case no time-marching
c     xg = xgorg
c     yg = ygorg
c     zg = zgorg

c     e1x = e1xorg
c     e1y = e1yorg
c     e1z = e1zorg
c     e2x = e2xorg
c     e2y = e2yorg
c     e2z = e2zorg
c     e3x = e3xorg
c     e3y = e3yorg
c     e3z = e3zorg

c     vgx = vgxorg
c     vgy = vgyorg
c     vgz = vgzorg

c     ome1 = ome1org
c     ome2 = ome2org
c     ome3 = ome3org

      t = t + 2.0d0*dt

      if((flagoutr.eq.2).or.
     &  ((flagoutr.eq.1).and.(loop.eq.lastloop))) then
      open(unit=10,file='Output_data/position_x.dat'
     &            ,status='old',access='append')
      write(10,*) t, xg
      close(10)

      open(unit=10,file='Output_data/position_y.dat'
     &            ,status='old',access='append')
      write(10,*) t, yg
      close(10)

      open(unit=10,file='Output_data/position_z.dat'
     &            ,status='old',access='append')
      write(10,*) t, zg
      close(10)

      open(unit=10,file='Output_data/principal_axis_direction.dat'
     &            ,status='old',access='append')
 9011 format(1p4e18.8e3)
      write(10,9011) t, e1x, e1y, e1z 
      write(10,9011) t, e2x, e2y, e2z 
      write(10,9011) t, e3x, e3y, e3z 
      close(10)

      open(unit=10,file='Output_data/x-y_locus.dat'
     &            ,status='old',access='append')
      write(10,*) xg, yg
      close(10)

      open(unit=10,file='Output_data/velocity_x.dat'
     &            ,status='old',access='append')
      write(10,*) t, vgx
      close(10)

      open(unit=10,file='Output_data/velocity_y.dat'
     &            ,status='old',access='append')
      write(10,*) t, vgy
      close(10)

      open(unit=10,file='Output_data/velocity_z.dat'
     &            ,status='old',access='append')
      write(10,*) t, vgz
      close(10)

      open(unit=10,file='Output_data/angular_velocity_1.dat'
     &            ,status='old',access='append')
      write(10,*) t, ome1
      close(10)

      open(unit=10,file='Output_data/angular_velocity_2.dat'
     &            ,status='old',access='append')
      write(10,*) t, ome2
      close(10)

      open(unit=10,file='Output_data/angular_velocity_3.dat'
     &            ,status='old',access='append')
      write(10,*) t, ome3
      close(10)

      open(unit=10,file='Output_data/acceleration_x.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(1,1)/6.0d0 + fvd(1,2)/3.0d0
     &              + fvd(1,3)/3.0d0 + fvd(1,4)/6.0d0)
      close(10)

      open(unit=10,file='Output_data/acceleration_y.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(2,1)/6.0d0 + fvd(2,2)/3.0d0
     &              + fvd(2,3)/3.0d0 + fvd(2,4)/6.0d0)
      close(10)

      open(unit=10,file='Output_data/acceleration_z.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(3,1)/6.0d0 + fvd(3,2)/3.0d0
     &              + fvd(3,3)/3.0d0 + fvd(3,4)/6.0d0)
      close(10)

      open(unit=10,file='Output_data/angular_acceleration_1.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(4,1)/6.0d0 + fvd(4,2)/3.0d0
     &              + fvd(4,3)/3.0d0 + fvd(4,4)/6.0d0)
      close(10)

      open(unit=10,file='Output_data/angular_acceleration_2.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(5,1)/6.0d0 + fvd(5,2)/3.0d0
     &              + fvd(5,3)/3.0d0 + fvd(5,4)/6.0d0)
      close(10)

      open(unit=10,file='Output_data/angular_acceleration_3.dat'
     &            ,status='old',access='append')
      write(10,*) t, (fvd(6,1)/6.0d0 + fvd(6,2)/3.0d0
     &              + fvd(6,3)/3.0d0 + fvd(6,4)/6.0d0)
      close(10)
      end if

*************************************************
* judgement of penetration
*************************************************

      do j=0,10
          pene(j) = 0
      enddo

      if((flagpene(0).eq.1).and.(loop.eq.1))then

* compute maximum distance in each body segment
      do j=1,numb
          dismax(j) = 0.0d0
          i=1
          do k=1,divt-1
              dis = (xbm(j,i,k)-xbm(j,i,k+1))**2
     &             +(ybm(j,i,k)-ybm(j,i,k+1))**2
     &             +(zbm(j,i,k)-zbm(j,i,k+1))**2
              if(dis.gt.dismax(j))then
                  dismax(j) = dis
              endif
          enddo

          i=div
          do k=1,divt-1
              dis = (xbm(j,i,k)-xbm(j,i,k+1))**2
     &             +(ybm(j,i,k)-ybm(j,i,k+1))**2
     &             +(zbm(j,i,k)-zbm(j,i,k+1))**2
              if(dis.gt.dismax(j))then
                  dismax(j) = dis
              endif
          enddo

          dis = (xbm(j,1,1)-xbm(j,2,1))**2
     &         +(ybm(j,1,1)-ybm(j,2,1))**2
     &         +(zbm(j,1,1)-zbm(j,2,1))**2
          if(dis.gt.dismax(j))then
              dismax(j) = dis
          endif

          dis = (xbm(j,1,divt)-xbm(j,2,divt))**2
     &         +(ybm(j,1,divt)-ybm(j,2,divt))**2
     &         +(zbm(j,1,divt)-zbm(j,2,divt))**2
          if(dis.gt.dismax(j))then
              dismax(j) = dis
          endif

      enddo

* main judging loop start
      do j3=1,10

          if(j3.eq.1)then
* (1) right upper limb - trunk
          js = 18
          je = 20
          jp =  2
          jjs=  1
          jje=  9
          jjp=  1
      elseif(j3.eq.2)then
* (2) left upper limb - trunk
          js = 19
          je = 21
          jp =  2
          jjs=  1
          jje=  9
          jjp=  1
      elseif(j3.eq.3)then
* (3) right upper limb - left upper limb
          js = 18
          je = 20
          jp =  2
          jjs= 19
          jje= 21
          jjp=  2
      elseif(j3.eq.4)then
* (4) right lower limb - trunk
          js = 12
          je = 14
          jp =  2
          jjs=  1
          jje=  9
          jjp=  1
      elseif(j3.eq.5)then
* (5) left lower limb - trunk
          js = 13
          je = 15
          jp =  2
          jjs=  1
          jje=  9
          jjp=  1
      elseif(j3.eq.6)then
* (6) right lower limb - left lower limb
          js = 12
          je = 14
          jp =  2
          jjs= 13
          jje= 15
          jjp=  2
      elseif(j3.eq.7)then
* (7) right upper limb - right lower limb
          js = 16
          je = 20
          jp =  2
          jjs= 10
          jje= 14
          jjp=  2
      elseif(j3.eq.8)then
* (8) left upper limb - left lower limb
          js = 17
          je = 21
          jp =  2
          jjs= 11
          jje= 15
          jjp=  2
      elseif(j3.eq.9)then
* (9) right upper limb - left lower limb
          js = 16
          je = 20
          jp =  2
          jjs= 11
          jje= 15
          jjp=  2
      elseif(j3.eq.10)then
* (10) left upper limb - right lower limb
          js = 17
          je = 21
          jp =  2
          jjs= 10
          jje= 14
          jjp=  2
      else
      endif

      if(flagpene(j3).eq.1)then
          do  j= js, je, jp
          do jj=jjs,jje,jjp
              if(dismax(j).gt.dismax(jj))then
                  dismax0 = dismax(j)
              else
                  dismax0 = dismax(jj)
              endif
              do i=1,div
              do k=1,divt
              do ii=1,div
              do kk=1,divt
                  dis = (xbm(j,i,k)-xbm(jj,ii,kk))**2
     &                 +(ybm(j,i,k)-ybm(jj,ii,kk))**2
     &                 +(zbm(j,i,k)-zbm(jj,ii,kk))**2
                  if(dis.lt.dismax0*0.5d0)then
c                     pene(j3) = pene(j3) + dsqrt(dismax0)
                      pene(j3) = pene(j3) + 1
                  endif
              enddo
              enddo
              enddo
              enddo

          enddo
          enddo
      else
          pene(j3) = -1
      endif

      enddo


      do j=1,10
          if(pene(j).ge.0)then
              pene(0) = pene(0) + pene(j)
          endif
      enddo

      open(unit=10,file='Output_data/penetration.dat'
     &            ,status='old',access='append')
      write(10,6666) pene(0),pene(1),pene(2),pene(3),pene(4),pene(5)
     &            ,pene(6),pene(7),pene(8),pene(9),pene(10)
 6666 format(11i8)
      close(10)

      endif





      else
          write(*,*) 'No way!'
          stop
      end if

*************************************************
* calculate force and moment 
* due to inertia force and gravitational force
* (!!!miura)
*************************************************
      if(((flagoutj.eq.2).and.(rkloop.eq.4)).or.
     & ((flagoutj.eq.1).and.((loop.eq.lastloop).or.(loop.eq.lastloop-1))
     &                   .and.(rkloop.eq.4)))then
      axg  = fv(1)
      ayg  = fv(2)
      azg  = fv(3)
      aw1g = fv(4)
      aw2g = fv(5)
      aw3g = fv(6)

      wxg = ome1org*e1xorg + ome2org*e2xorg + ome3org*e3xorg
      wyg = ome1org*e1yorg + ome2org*e2yorg + ome3org*e3yorg
      wzg = ome1org*e1zorg + ome2org*e2zorg + ome3org*e3zorg

* calculate acceleration by finite differenciation
* (this part is added for version 2.3.0 2006/10/10 by Motomu Nakashima.)
      if (vgxsolv.eq.0) then
          if (kd.eq.numl-1) then
              axg = (vgxg(1)-vgxg(kd))/(2.0d0*dt)
          else
              axg = (vgxg(kd+2)-vgxg(kd))/(2.0d0*dt)
          end if
      end if
      if (vgysolv.eq.0) then
          if (kd.eq.numl-1) then
              ayg = (vgyg(1)-vgyg(kd))/(2.0d0*dt)
          else
              ayg = (vgyg(kd+2)-vgyg(kd))/(2.0d0*dt)
          end if
      end if
      if (vgzsolv.eq.0) then
          if (kd.eq.numl-1) then
              azg = (vgzg(1)-vgzg(kd))/(2.0d0*dt)
          else
              azg = (vgzg(kd+2)-vgzg(kd))/(2.0d0*dt)
          end if
      end if

      if (ome1solv.eq.0) then
          if (kd.eq.numl-1) then
              aw1g = (ome1g(1)-ome1g(kd))/(2.0d0*dt)
          else
              aw1g = (ome1g(kd+2)-ome1g(kd))/(2.0d0*dt)
          end if
      end if
      if (ome2solv.eq.0) then
          if (kd.eq.numl-1) then
              aw2g = (ome2g(1)-ome2g(kd))/(2.0d0*dt)
          else
              aw2g = (ome2g(kd+2)-ome2g(kd))/(2.0d0*dt)
          end if
      end if
      if (ome3solv.eq.0) then
          if (kd.eq.numl-1) then
              aw3g = (ome3g(1)-ome3g(kd))/(2.0d0*dt)
          else
              aw3g = (ome3g(kd+2)-ome3g(kd))/(2.0d0*dt)
          end if
      end if

*inertia -> absolute
      awxg = aw1g*e1xorg + aw2g*e2xorg + aw3g*e3xorg 
      awyg = aw1g*e1yorg + aw2g*e2yorg + aw3g*e3yorg 
      awzg = aw1g*e1zorg + aw2g*e2zorg + aw3g*e3zorg 

      do 500 j=1,numb
         u(1) = awxg
         u(2) = awyg
         u(3) = awzg

*tip-a
         v(1) = x0a(j) - xgorg
         v(2) = y0a(j) - ygorg
         v(3) = z0a(j) - zgorg

         axifa(j) = ax0a(j) + axg + exprox(u,v)
         ayifa(j) = ay0a(j) + ayg + exproy(u,v)
         azifa(j) = az0a(j) + azg + exproz(u,v)

*tip-b
         v(1) = x0b(j) - xgorg
         v(2) = y0b(j) - ygorg
         v(3) = z0b(j) - zgorg

         axifb(j) = ax0b(j) + axg + exprox(u,v)
         ayifb(j) = ay0b(j) + ayg + exproy(u,v)
         azifb(j) = az0b(j) + azg + exproz(u,v)

         do 501 i = 1,div
* obtain mass of thin elliptic plate
            dlen = len(j) /dble(div)
            r1 = ra1(j) + (rb1(j)-ra1(j))*(dble(i)-0.5d0)/dble(div)
            r2 = ra2(j) + (rb2(j)-ra2(j))*(dble(i)-0.5d0)/dble(div)
            dm = pi * r1 * r2 * dlen * rho(j)

* obtain center position of ellipse
            xc = x0a(j) + (x0b(j)-x0a(j))*(dble(i)-0.5d0)/dble(div)
            yc = y0a(j) + (y0b(j)-y0a(j))*(dble(i)-0.5d0)/dble(div)
            zc = z0a(j) + (z0b(j)-z0a(j))*(dble(i)-0.5d0)/dble(div)

* obtain center acceleration of ellipse
            axc = axifa(j) 
     &          + (axifb(j)-axifa(j))*(dble(i)-0.5d0)/dble(div)
            ayc = ayifa(j)
     &          + (ayifb(j)-ayifa(j))*(dble(i)-0.5d0)/dble(div)
            azc = azifa(j)
     &          + (azifb(j)-azifa(j))*(dble(i)-0.5d0)/dble(div)

* calculate inertia force
            dfix = -dm*axc
            dfiy = -dm*ayc
            dfiz = -dm*azc

            do igomi=1,nm
            gomi(igomi)=0.0d0
            enddo

* obtain force and moment between segments
* for calculation of joint torque
            call bodyforce(dfix,dfiy,dfiz,fxup,fyup,fzup
     &                    ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &                    ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &                    ,xgorg,ygorg,zgorg)

* calculte gravitational force
            dfix = 0.0d0
            dfiy = 0.0d0
            dfiz =  - dm*grav

* obtain force and moment between segments
* for calculation of joint torque
            call bodyforce(dfix,dfiy,dfiz,fxup,fyup,fzup
     &                    ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &                    ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &                    ,xgorg,ygorg,zgorg)

************************************************************
* inertial torque about mass center of thin elliptic plate
************************************************************
c           di1 = pi*r1*(r2**3)*dlen*rho(j)*0.25d0
c           di2 = pi*r2*(r1**3)*dlen*rho(j)*0.25d0
c           di3 = di1+di2

c           call inertorque(e01x(j),e01y(j),e01z(j)
c    &                     ,e02x(j),e02y(j),e02z(j)
c    &                     ,e03x(j),e03y(j),e03z(j)
c    &                     ,di1,di2,di3,dm
c    &                     ,awxg,awyg,awzg
c    &                     ,wxg,wyg,wzg
c    &                     ,mxup,myup,mzup,mxdw,mydw,mzdw
c    &                     ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,i,j
c    &                     ,xgorg,ygorg,zgorg,inerh,dt*2.0d0)

************************************************
* inertial force due to added mass(all) 
************************************************
* obtain vector from tip-a to tip-b
            xba = x0b(j)-x0a(j)
            yba = y0b(j)-y0a(j)
            zba = z0b(j)-z0a(j)

* normalize vector from tip-a to tip-b
            baabs = vecabs(xba,yba,zba)
            xba = xba/baabs
            yba = yba/baabs
            zba = zba/baabs

* obtain tangential component of acceleration 
* using inner product against unit vector from tip-a to tip-b
            inproh = inpro(axc,ayc,azc,xba,yba,zba)
            axct = inproh * xba
            ayct = inproh * yba
            azct = inproh * zba

* obtain normal component by deducting tangential component
            axcn = axc - axct
            aycn = ayc - ayct
            azcn = azc - azct

* obtain two components parallel to two axes of ellipse
* from inner product
* axis-1
            inproi1 = inpro(axcn,aycn,azcn,e01x(j),e01y(j),e01z(j))
            axcn1 = inproi1 * e01x(j)
            aycn1 = inproi1 * e01y(j)
            azcn1 = inproi1 * e01z(j)
* axis-2
            inproi2 = inpro(axcn,aycn,azcn,e02x(j),e02y(j),e02z(j))
            axcn2 = inproi2 * e02x(j)
            aycn2 = inproi2 * e02y(j)
            azcn2 = inproi2 * e02z(j)

* obtain fluid force
            dfix =
     &        - dlen*pi*r2**2 * (numen2(j,i)/dble(divt)) * axcn1 * ca(j)
     &        - dlen*pi*r1**2 * (numen2(j,i)/dble(divt)) * axcn2 * ca(j)
            dfiy =
     &        - dlen*pi*r2**2 * (numen2(j,i)/dble(divt)) * aycn1 * ca(j)
     &        - dlen*pi*r1**2 * (numen2(j,i)/dble(divt)) * aycn2 * ca(j)
            dfiz =
     &        - dlen*pi*r2**2 * (numen2(j,i)/dble(divt)) * azcn1 * ca(j)
     &        - dlen*pi*r1**2 * (numen2(j,i)/dble(divt)) * azcn2 * ca(j)

* obtain force and moment between segments
* for calculation of joint torque
            call bodyforce(dfix,dfiy,dfiz,fxup,fyup,fzup
     &                    ,fxdw,fydw,fzdw,mxup,myup,mzup,mxdw,mydw,mzdw
     &                    ,x0a,y0a,z0a,x0b,y0b,z0b,xc,yc,zc,j
     &                    ,xgorg,ygorg,zgorg)

  501    continue
  500 continue

* torque component due to distance between calculation surface
* and mass center and force acting on the surface
      call rexfa(fxup,fyup,fzup,fxdw,fydw,fzdw
     &          ,mxup,myup,mzup,mxdw,mydw,mzdw
     &          ,x0a,y0a,z0a,x0b,y0b,z0b,xgorg,ygorg,zgorg)

      endif
**********************************************************
* transform cooridinate into direction of lower waist(1)
**********************************************************
      call 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(1),e01y(1),e01z(1)
     &     ,e02x(1),e02y(1),e02z(1),e03x(1),e03y(1),e03z(1))

**********
* output
**********
      if((flagoutj.eq.2).or.
     &  ((flagoutj.eq.1).and.(loop.eq.lastloop))) then
      if (rkloop.eq.4) then
      call outjoint(numb,t-2.0d0*dt,mxup,myup,mzup,mxdw,mydw,mzdw
     &              ,mxbup,mybup,mzbup,mxbdw,mybdw,mzbdw)
      endif
      endif

********************
* finish time step 
********************
 1100 continue
 1010 continue

**************************************************
* calculate propulsive distance during one cycle
**************************************************
      stlen = dsqrt(
     &         (xg-xgloop)**2
     &        +(yg-ygloop)**2
     &        +(zg-zgloop)**2
     &             )
      xst = xg - xgloop
      yst = yg - ygloop
      direc2 = datan(yst/xst)
      muki = direc2 * 180.0d0/pi

      if (xst.lt.0.0d0) then
          muki = muki + 180.0d0
      end if

* output camera information for multi project mode
      open(unit=10,file='Output_data/camera.dat',status='unknown')
      write(10,*) xgloop,ygloop
      write(10,*) xst,yst
      write(10,*) muki
      close(10)

***************************************
***************************************
* Waiting loop for multi project mode *
***************************************
***************************************
      if (multi.ge.1) then

        goorstop = -2
 5555   continue
        open(unit=10,file='go2.dat',status='replace',err=5555)
 5557   continue
        write(10,*,err=5557) '0 '
        close(10)

      do i=1,4000000

        open(unit=10,file='go.dat',action='read'
     &         ,status='unknown',err=5552)
        read(10,*,end=5553,err=5554) goorstop, dum1, dum2, dum3
        close(10)
        if ((goorstop.eq.1)
     &      .and.(dum1.eq.loop)
     &      .and.(dum2.eq.-1)
     &      .and.(dum3.eq.rkloop)) then
          goto 5551
        endif

 5552   continue
 5553   continue
 5554   continue

        if ((i.ge.500).and.
     &      (mod(i,100).eq.0)) then
          write(*,*) 'waiting iteration=',i, goorstop,dum1,dum2,dum3
          write(*,*) 'waiting iteration=',i, 1,loop,-1,rkloop
        endif

*       waiting
        call sleepqq(1)

      enddo
      write(*,*) 'Waiting loop was finished. waited too long.'
      stop

 5551 continue

        goorstop = -2
 5556   continue
        open(unit=10,file='go2.dat',status='replace',err=5556)
 5558   continue
        write(10,*,err=5558) '1 '
        close(10)

      endif



*********************
* output motion.dat *
*********************
      if((flagouta.eq.2).or.
     &  ((flagouta.eq.1).and.(loop.eq.lastloop))) then

      open(unit=20,file='Output_data/motion.dat',status='unknown')
      open(unit=25,file='Output_data/cogcob.dat',status='unknown')
      write(20,*) numl/2, loop, lastloop
      write(20,*) numb
      write(20,*) div
      write(20,*) stlen
      write(20,*) muki


* only sub project
      if (multi.eq.1) then
* read camera information of main project
      open(unit=10,file='../main/Output_data/camera.dat'
     &            ,status='unknown',err=5562)
      read(10,*,end=5562,err=5562) xgloop,ygloop
      read(10,*,end=5562,err=5562) xst,yst
      read(10,*,end=5562,err=5562) muki
      close(10)
      goto 5561
 5562 continue
* if failed, read original camera information again
      open(unit=10,file='Output_data/camera.dat'
     &            ,status='unknown')
      read(10,*) xgloop,ygloop
      read(10,*) xst,yst
      read(10,*) muki
      close(10)
 5561 continue
      end if

      do kd=1,numl,2
      xcamera = xgloop + xst * dble(kd-1)/dble(numl)
      ycamera = ygloop + yst * dble(kd-1)/dble(numl)

      do 490 j=1,numb
      x0at= (x0ac(j,kd)-xcamera)*dcos(pi-muki*pi/180.0d0)
     &     -(y0ac(j,kd)-ycamera)*dsin(pi-muki*pi/180.0d0)
      y0at= (x0ac(j,kd)-xcamera)*dsin(pi-muki*pi/180.0d0)
     &     +(y0ac(j,kd)-ycamera)*dcos(pi-muki*pi/180.0d0)
      x0bt= (x0bc(j,kd)-xcamera)*dcos(pi-muki*pi/180.0d0)
     &     -(y0bc(j,kd)-ycamera)*dsin(pi-muki*pi/180.0d0)
      y0bt= (x0bc(j,kd)-xcamera)*dsin(pi-muki*pi/180.0d0)
     &     +(y0bc(j,kd)-ycamera)*dcos(pi-muki*pi/180.0d0)

      e01xt= e01xc(j,kd)*dcos(pi-muki*pi/180.0d0)
     &      -e01yc(j,kd)*dsin(pi-muki*pi/180.0d0)
      e01yt= e01xc(j,kd)*dsin(pi-muki*pi/180.0d0)
     &      +e01yc(j,kd)*dcos(pi-muki*pi/180.0d0)
      e02xt= e02xc(j,kd)*dcos(pi-muki*pi/180.0d0)
     &      -e02yc(j,kd)*dsin(pi-muki*pi/180.0d0)
      e02yt= e02xc(j,kd)*dsin(pi-muki*pi/180.0d0)
     &      +e02yc(j,kd)*dcos(pi-muki*pi/180.0d0)

      write(20,8000) x0at,y0at,z0ac(j,kd),x0bt,y0bt,z0bc(j,kd)
     &              ,e01xt,e01yt,e01zc(j,kd),e02xt,e02yt,e02zc(j,kd)
 8000 format(1p12e15.5e3)
  490 continue

      do 492 j=1,numb
      do 492 j2=1,div
      xdistt= (xdistc(j,j2,kd)-xcamera)*dcos(pi-muki*pi/180.0d0)
     &       -(ydistc(j,j2,kd)-ycamera)*dsin(pi-muki*pi/180.0d0)
      ydistt= (xdistc(j,j2,kd)-xcamera)*dsin(pi-muki*pi/180.0d0)
     &       +(ydistc(j,j2,kd)-ycamera)*dcos(pi-muki*pi/180.0d0)

      fxdistt= fxdistc(j,j2,kd)*dcos(pi-muki*pi/180.0d0)
     &        -fydistc(j,j2,kd)*dsin(pi-muki*pi/180.0d0)
      fydistt= fxdistc(j,j2,kd)*dsin(pi-muki*pi/180.0d0)
     &        +fydistc(j,j2,kd)*dcos(pi-muki*pi/180.0d0)

      write(20,8001)  xdistt, ydistt, zdistc(j,j2,kd)
     &              ,fxdistt,fydistt,fzdistc(j,j2,kd)
 8001 format(1p6e15.5e3)
  492 continue

* writing cogcob.dat
      xgorgt= (xgorgd(kd)-xcamera)*dcos(pi-muki*pi/180.0d0)
     &       -(ygorgd(kd)-ycamera)*dsin(pi-muki*pi/180.0d0)
      ygorgt= (xgorgd(kd)-xcamera)*dsin(pi-muki*pi/180.0d0)
     &       +(ygorgd(kd)-ycamera)*dcos(pi-muki*pi/180.0d0)
      xcobt = (xcob(kd)-xcamera)*dcos(pi-muki*pi/180.0d0)
     &       -(ycob(kd)-ycamera)*dsin(pi-muki*pi/180.0d0)
      ycobt = (xcob(kd)-xcamera)*dsin(pi-muki*pi/180.0d0)
     &       +(ycob(kd)-ycamera)*dcos(pi-muki*pi/180.0d0)
      e1xorgt= (e1xorgd(kd))*dcos(pi-muki*pi/180.0d0)
     &        -(e1yorgd(kd))*dsin(pi-muki*pi/180.0d0)
      e1yorgt= (e1xorgd(kd))*dsin(pi-muki*pi/180.0d0)
     &        +(e1yorgd(kd))*dcos(pi-muki*pi/180.0d0)
      e2xorgt= (e2xorgd(kd))*dcos(pi-muki*pi/180.0d0)
     &        -(e2yorgd(kd))*dsin(pi-muki*pi/180.0d0)
      e2yorgt= (e2xorgd(kd))*dsin(pi-muki*pi/180.0d0)
     &        +(e2yorgd(kd))*dcos(pi-muki*pi/180.0d0)
      e3xorgt= (e3xorgd(kd))*dcos(pi-muki*pi/180.0d0)
     &        -(e3yorgd(kd))*dsin(pi-muki*pi/180.0d0)
      e3yorgt= (e3xorgd(kd))*dsin(pi-muki*pi/180.0d0)
     &        +(e3yorgd(kd))*dcos(pi-muki*pi/180.0d0)

      write(25,8003) xgorgt, ygorgt, zgorgd(kd), xcobt, ycobt, zcob(kd)
 8003 format(1p6e15.5e3)
      write(25,8004) e1xorgt,e1yorgt,e1zorgd(kd)
      write(25,8004) e2xorgt,e2yorgt,e2zorgd(kd)
      write(25,8004) e3xorgt,e3yorgt,e3zorgd(kd)
 8004 format(1p3e15.5e3)

      enddo
      end if

          
      write(*,*) '    direction              ',muki
      write(*,*) '    x                      ',xst
      write(*,*) '    y                      ',yst
      write(*,*) '    stroke length          ',stlen
      write(*,*) '    stroke length deviation'
     &                ,dabs((stlen-stlenold)/stlenold)

      open(unit=10,file='Output_data/swimming_speed.dat'
     &            ,status='old',access='append')
      write(10,*) loop, stlen
      close(10)

      open(unit=10,file='Output_data/stroke_length.dat'
     &            ,status='old',access='append')
      write(10,*) loop, stlen
      close(10)

      open(unit=10,file='Output_data/direction.dat'
     &            ,status='old',access='append')
      write(10,*) loop, muki*pi/180.0d0
      close(10)

      open(unit=10,file='Output_data/power_efficiency1.dat'
     &            ,status='old',access='append')
      write(10,*) loop, -draggli*stlen**3/poweraav
      close(10)

      open(unit=10,file='Output_data/power_efficiency2.dat'
     &            ,status='old',access='append')
      write(10,*) loop, powerpav/poweraav
      close(10)

      open(unit=10,file='Output_data/input_power.dat'
     &            ,status='old',access='append')
      write(10,*) loop, -poweraav
      close(10)

      open(unit=10,file='Output_data/output_power1.dat'
     &            ,status='old',access='append')
      write(10,*) loop, draggli*stlen**3
      close(10)

      open(unit=10,file='Output_data/output_power2.dat'
     &            ,status='old',access='append')
      write(10,*) loop, -powerpav
      close(10)

      xgloop = xg
      ygloop = yg
      zgloop = zg
      stlenold = stlen
      poweraav = 0.0d0
      powerpav = 0.0d0

      powtxav  = 0.0d0
      pownxav  = 0.0d0
      powixav  = 0.0d0
      powtyav  = 0.0d0
      pownyav  = 0.0d0
      powiyav  = 0.0d0
      powtzav  = 0.0d0
      pownzav  = 0.0d0
      powizav  = 0.0d0

 1000 continue

      open(unit=10,file='end_flag.dat',status='unknown')
      write(10,*) "1"
      close(10)

*     closing motion.dat
      close(20)
      close(25)

      stop
      end



