c--------------------------------------------------------------------+
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567--------------------------------------------------------------+
      include 'common_vib8.fil'
      real*8    px1(20,10000),px2(20,10000),
     &          py1(20,10000),py2(20,10000),
     &          ux(3,35,10000)
c
      character*14 afile,ifile,bfile,ffile,sfile
      character*14 ofile(100)
      character*14 vfile,nfile,infile
C 8/27/13
      common /inflow/Vin
C 8/16/13
      print *,' 8/16/13 ***** vfile format changed to formatted'
      print *,' 8/20/13 initsout,tsout to text out time series'
      print *,' 8/27/13 main,rdvelo: common/inflow/Vin'
      print *,'                set inflow vel to Vin'
      print *,' 8/29/13   text out vib in tsout'
      print *,' 8/30/13 uvwbdy: duc=0 for fixed body case'
      print *,' 8/31/13 uvwbdy: set vv(i,je,k)=vv(i,js,k) if iyper=0'
      print *,'                 else vv(i,je,k)=vv(i,js,k)=0'
      print *,'         rdvelo: remove pp(is,j,k),pp(is-1,j,k)=0'
      print *,'         hsmac2: remove pp(is,j,k)=0'
      print *,'         prsbdy: subtract pp(is,je/2,ke/2) from all pp'
      print *,' 9/ 2/13 uvwbdy: vv(i,je,k)=vv(i,js,k)=-yvel, iyper=1'
      print *,' 9/ 3/13 rduvwi: skip reading infile if fffff'
      print *,' 9/ 5/13 main:  call outvib every time'
      print *,'  ***************************************** '
c
c--------------------------------------------- 2012.05.15 rewrite
c      is=2
C 8/27/13
      print *,' enter inflow velocity, Vin'
      read (5,*) Vin
      is=3
      js=3
      ks=3
      t=0.0d0
c
      itm2=2
c
      call rdinit(afile,bfile,ifile,ffile,sfile,ofile,infile,nfile
     &           ,vfile,ifc,kout,nfl)
      call intsum
      call rdaxis(afile)
      call rdflag(bfile)
      call rduvwi(ifile,infile)
      if(ifile(1:5).eq.'fffff') call hsmac2(dtmin)
c
c
      print *,'dt=',dt
      print *,'T =',t
      ifcs=ifc
      icount=0
      nn=0
C 8/20/13 initialize time series output
      call initsout(t)
c
      ndivf1=ndivf
      nhcnt=0
      itm2c=0
      nn2=0 
c***********************************************************
 1000 nn=nn+1
c***********************************************************
      t=t+dt
      call prefgh
      call space3
      call euler2(ifile)
      call hsmac2(dtmin)
      call smagor
c************************************* add for cal. vibration
      call calfrm
      call calvib
c*************************************
c
      nn2=nn/2
      itm2c=itm2c+1
      print *,'itm2c nn2 nn =',itm2c,nn2,nn
c      if(itm2c.eq.itm2) then
c       call spcdat(px1,px2,py1,py2,ux,nn2)
c       itm2c=0
c      end if
c
      call change(dtmin)
      call calstt
c
      if(nn.eq.ifcs) then
       icount=icount+1
       if(icount.le.nfl) then
         call outdat(ofile,kout,icount,ifcs)
         call finalo(ffile)
C 9/ 5/13 output vib every time
         call outvib(vfile,nfile)
       end if
       ifcs=ifcs+ifc
      end if
C 8/20/13 time serie output every 10 steps
C 8/29/13 add xdis, ydis to output 7 quantities
      if(mod(nn,10).eq.0) then
        call tsout(nn,t,7,uu(is+1,je/2,ke/2),pfx,pfy,pft,
     &    xdis,ydis)
      endif
c
      print *,'nn = ',nn
c***********************************************************
      if(nn.lt.nt) go to 1000
c***********************************************************
c
c      call spcout(px1,px2,py1,py2,ux)
c
      call finalo(ffile)
      istt=nt
      call calfin(istt)
      call outstt(sfile)
c************************************* add for cal. vibration 
      call outvib(vfile,nfile)
c*************************************
c
      close(11)
      stop 
      end
c 8/20/13 ----------------------------------------- sub.initsout
      subroutine initsout(t)
c234567-------------------------------------------- sub.initsout
      double precision t
C  initialize time-series output
      if(t.eq.0.) open (25,file='DAT/breeze_vib8'//'.ts',
     &  status='UNKNOWN', form='formatted')
      if(t.gt.0.) then
        print *,' ts output append to exsisting file',
     &        ' or start from t=',t,'?(0/1)'
        read (5,*) itso
        if(itso.eq.0) open (25,file='DAT/breeze_vib8'//'.ts',
     &    status='OLD', form='formatted',access='APPEND')
        if(itso.eq.1) open (25,file='DAT/breeze_vib8'//'.ts',
     &    status='UNKNOWN', form='formatted')
      endif
      write (25,*) '#start at t=',t
C 8/29/13 add displacement 
      write (25,*) 
     & '#nn       t         Vin        pfx        pfy        pft'
     &,'      xdis      ydis' 
      print *, 'time series will be output in ',
     &    'DAT/breeze_vib8'//'.ts'
      close(25)
      return
      end
c 8/29/13 8/20/13  --------------------------------------- sub.tsout
      subroutine tsout(nn,t,nq,q1,q2,q3,q4,q5,q6)
c234567--------------------------------------------------- sub.tsout
      double precision t,q1,q2,q3,q4,q5,q6
C  output time series in unit 25
      open (25,file='DAT/breeze_vib8'//'.ts',
     &  status='OLD', form='formatted',access='APPEND')
      write(25,100) nn,t,q1,q2,q3,q4,q5,q6
  100 FORMAT(I6,F10.4,17E12.4)
      close(25)
      return
      end
c---------------------------------------------------------- sub.intsum
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.intsum 
      subroutine intsum
      include 'common_vib8.fil'
c
      do 10 k=1,kmax
      do 10 j=1,jmax
      do 10 i=1,imax
       aveu(i,j,k)=0.0d0
       avev(i,j,k)=0.0d0
       avew(i,j,k)=0.0d0
       avep(i,j,k)=0.0d0
       avet(i,j,k)=0.0d0
       rmsp(i,j,k)=0.0d0
       rs(i,j,k,1)=0.0d0
       rs(i,j,k,2)=0.0d0
       rs(i,j,k,3)=0.0d0
       rs(i,j,k,4)=0.0d0
       rs(i,j,k,5)=0.0d0
       rs(i,j,k,6)=0.0d0
       rs(i,j,k,7)=0.0d0
       rs(i,j,k,8)=0.0d0
       rs(i,j,k,9)=0.0d0
c
       rss(i,j,k,1)=0.0d0
       rss(i,j,k,2)=0.0d0
       rss(i,j,k,3)=0.0d0
   10 continue
c
      do 20 i=1,imax
       x(i)=0.0d0
       dx(i)=0.0d0
       dxh(i)=0.0d0
   20 continue
      do 30 j=1,jmax
       y(j)=0.0d0
       dy(j)=0.0d0
       dyh(j)=0.0d0
   30 continue
      do 40 k=1,kmax
       z(k)=0.0d0
       dz(k)=0.0d0
       dzh(k)=0.0d0
   40 continue
c
      do 50 k=1,kmax
      do 50 j=1,jmax
      do 50 i=1,imax
       uu(i,j,k)=0.0d0
       vv(i,j,k)=0.0d0
       ww(i,j,k)=0.0d0
       pp(i,j,k)=0.0d0
c
       u(i,j,k)=0.0d0
       v(i,j,k)=0.0d0
       w(i,j,k)=0.0d0
       p(i,j,k)=0.0d0
c
       fp(i,j,k)=0.0d0
       gp(i,j,k)=0.0d0
       hp(i,j,k)=0.0d0
c
       f(i,j,k)=0.0d0
       g(i,j,k)=0.0d0
       h(i,j,k)=0.0d0
c
       tnu(i,j,k)=0.0d0
       ibf(i,j,k)=0
   50 continue
c
      do 60 n=1,nmax
       frmt(1,n)=0.0d0
       frmt(2,n)=0.0d0
       frmt(3,n)=0.0d0
c
       resp(1,1,n)=0.0d0
       resp(1,2,n)=0.0d0
       resp(1,3,n)=0.0d0
c
       resp(2,1,n)=0.0d0
       resp(2,2,n)=0.0d0
       resp(2,3,n)=0.0d0
c
       resp(3,1,n)=0.0d0
       resp(3,2,n)=0.0d0
       resp(3,3,n)=0.0d0
c
   60 continue
c
      do 70 k=1,kmax
      do 70 j=1,jmax
      do 70 i=1,imax
       facc0(i,j,k)=0.0d0
       gacc0(i,j,k)=0.0d0
       facc1(i,j,k)=0.0d0
       gacc1(i,j,k)=0.0d0
   70 continue
      return
      end
c---------------------------------------------------------- sub.rdinit
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.rdinit
      subroutine rdinit
     &           (afile,bfile,ifile,ffile,sfile,ofile,infile,nfile,
     &            vfile,ifc,kout,nfl)
      include 'common_vib8.fil'
c
      character*14 afile,ifile,bfile,ffile,sfile
      character*14 ofile(100)
      character*14 vfile,nfile,infile
      
c
      open(10,file='breeze_vib8.inf',status='old',form='formatted')
       read(10,100)afile,bfile,ifile,ffile,sfile
       read(10,200)inum,jnum,knum
       read(10,300)iyper,izper
       read(10,400)re
       read(10,500)cs
       read(10,600)dt
       read(10,700)esp,omg
       read(10,700)ca1,ca2
       read(10,800)nt
       read(10,910)infile,nfile,vfile
       read(10,920)amas,aitm
       read(10,930)xkst,ykst,tkst
       read(10,940)xh0,yh0,th0
       read(10,950)beta,ndof
       read(10,960)ampxx,ampyy,amptt,tcntn
       read(10,800)ndivf

      close(10)
c
      open(10,file='breeze_vib8.fil',status='old',form='formatted')
       read(10,810)ifc,kout,nfl
       do 10 i=1,nfl
        read(10,820)ofile(i)
c         print *,ofile(i)
   10  continue
      close(10)
      print *,'ifc,kout,nfl = ',ifc,kout,nfl
c
  100 format(2x,a14,2x,a14,2x,a14,2x,a14,2x,a14)
  200 format(3i5)
  300 format(2i5)
  400 format(f10.0)
  500 format(f10.5)
  600 format(e13.5)
  700 format(2f10.5)
  800 format(i5)
  810 format(3i5)
  820 format(a14)
  910 format(2x,a14,2x,a14,2x,a14)
  920 format(2e13.5)
  930 format(3e13.5)
  940 format(3f10.5)
  950 format(f10.5,i5)
  960 format(4f10.5)
  
c
c--------------------------------------------- 2012.05.15 rewrite
c      ie=inum+1
      ie=inum+2
      je=jnum+2
      ke=knum+2
c
c************************************* add for cal. vibration
      pai=3.141592654d0
      xomg=dsqrt(xkst/amas)
      yomg=dsqrt(ykst/amas)
      tomg=dsqrt(tkst/aitm)
      xcvd=2.0d0*amas*xomg*xh0
      ycvd=2.0d0*amas*yomg*yh0
      tcvd=2.0d0*aitm*tomg*th0
c
      print *,'axis        file = ',afile
      print *,'flag        file = ',bfile
      print *,'input       file = ',ifile
      print *,'final       file = ',ffile
      print *,'statistics  file = ',sfile
      print *,'# of x,y,z  = ',ie-is+1,je-js+1,ke-ks+1
      print *,'Re Num.     = ',re
      print *,'Coef. Smg.  = ',cs
      print *,'Periodic in Y [yes=0,no=1  ]',iyper
      print *,'Periodic in Z [yes=0,no=1,2]',izper
      print *,'Dt      =',dt
      print *,'esp,omg =',esp,omg
      print *,'nt      =',nt
      print *,'ca1 ca2 =',ca1,ca2
      print *,'amas,aitm     =',amas,aitm
      print *,'xkst,ykst,tkst=',xkst,ykst,tkst
      print *,'xh0,yh0,th0   =',xh0,yh0,th0
      print *,'1/beta,ndof   =',1.0/beta,ndof
c
      ca1=1.0d0-(1.0d0-ca1)*2.0d0
      ca2=1.0d0-(1.0d0-ca2)*2.0d0
c
      print *,'ca1 ca2 =',ca1,ca2
c
      return 
      end
c---------------------------------------------------------- sub.rdaxis
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.rdaxis
      subroutine rdaxis(afile)
      include 'common_vib8.fil'
c
      character*14 afile
c
      open(10,file=afile,status='old',form='formatted')
       read(10,100)ic,jc,kc
        jd=0
c--------------------------------------------- 2012.05.15 rewrite
c        if(ie.ne.ic+1) jd=1
        if(ie.ne.ic+2) jd=1
        if(je.ne.jc+2) jd=2
        if(ke.ne.kc+2) jd=3
        if(jd.ne.0) then
         print *,'# of data are not same.'
         print *,'[1];x [2];y [3];z not same ',jd
         stop
        end if
c
c--------------------------------------------- 2012.05.15 rewrite
c       read(10,200)(x(i),i=is,ic+1)
c       read(10,200)(y(j),j=js,jc+2)
c       read(10,200)(z(k),k=ks,kc+2)
       read(10,200)(x(i),i=is,ie)
       read(10,200)(y(j),j=js,je)
       read(10,200)(z(k),k=ks,ke)
      close(10)
  100 format(3i5)
  200 format(5f10.5)
c
c---------------------------------- make dx,dxh
      do 10 i=is+1,ie
       dx(i)=x(i)-x(i-1)
   10 continue
      dx(is)=dx(is+1)
      dx(ie+1)=dx(ie)
c
      do 20 i=is,ie
       dxh(i)=0.5d0*(dx(i)+dx(i+1))
   20 continue
c
      x(is-1)=x(is)-dx(is)
      x(ie+1)=x(ie)+dx(ie+1)
c
c---------------------------------- make dy,dyh
      do 30 j=js+1,je
       dy(j)=y(j)-y(j-1)
   30 continue
c
      if(iyper.ne.0) then
       dy(js)=dy(js+1)
       dy(je+1)=dy(je)
        y(js-1)=y(js)-dy(js)
        y(je+1)=y(je)+dy(je+1)
      else
       dy(js)=dy(je)
       dy(js-1)=dy(je-1)
       dy(je+1)=dy(js+1)
       dy(je+2)=dy(js+2)
        y(js-1)=y(js)-dy(js)
        y(js-2)=y(js-1)-dy(js-1)
        y(je+1)=y(je)+dy(je+1)
        y(je+2)=y(je+1)+dy(je+2)
      end if
c
      do 40 j=js-1,je+1
       dyh(j)=0.5d0*(dy(j)+dy(j+1))
   40 continue
c
c---------------------------------- make dz,dzh
      do 50 k=ks+1,ke
       dz(k)=z(k)-z(k-1)
   50 continue
c
      if(izper.ne.0) then
       dz(ks)=dz(ks+1)
       dz(ke+1)=dz(ke)
        z(ks-1)=z(ks)-dz(ks)
        z(ke+1)=z(ke)+dz(ke+1)
      else
       dz(ks)=dz(ke)
       dz(ks-1)=dz(ke-1)
       dz(ke+1)=dz(ks+1)
       dz(ke+2)=dz(ks+2)
        z(ks-1)=z(ks)-dz(ks)
        z(ks-2)=z(ks-1)-dz(ks-1)
        z(ke+1)=z(ke)+dz(ke+1)
        z(ke+2)=z(ke+1)+dz(ke+2)
      end if
c
      do 60 k=ks-1,ke+1
       dzh(k)=0.5d0*(dz(k)+dz(k+1))
   60 continue
c
c
      dtmin=10000.0
      do 70 i=is,ie-1
      do 70 j=js,je-1
      do 70 k=ks,ke-1
       d1=re/2.0d0/(1.0d0/dx(i+1)**2+1.0d0/dy(j+1)**2+1.0d0/dz(k+1)**2) 
       dtmin=dmin1(dtmin,d1)
   70 continue
      delm=dtmin
c
      print *,'delm =',delm
c
c
      return
      end
c--------------------------------------------------------- sub. rdflag
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567--------------------------------------------------- sub. rdflag
      subroutine rdflag(bfile)
      include 'common_vib8.fil'
c
      character*14 bfile
c
      open(10,file=bfile,status='old',form='formatted')
       read(10,100)ic,jc,kc
        jd=0
c--------------------------------------------- 2012.05.15 rewrite
c        if(ie.ne.ic+1) jd=1
        if(ie.ne.ic+2) jd=1
        if(je.ne.jc+2) jd=2
        if(ke.ne.kc+2) jd=3
        if(jd.ne.0) then
         print *,'# of data are not same in rdflag.'
         print *,'[1];x [2];y [3];z not same ',jd
         stop
        end if
c
       read(10,200)(((ibf(i,j,k),i=is,ie-1),j=js,je-1),k=ks,ke-1)
      close(10)
  100 format(3i5)
  200 format(10i5)
c
      do 10 i=is,ie-1
      do 10 j=js,je-1
       if(ibf(i,j,ks).ne.0) then
        ibf(i,j,ks-1)=1
        ibf(i,j,ks-2)=1
       end if
       if(ibf(i,j,ke-1).ne.0) then
        ibf(i,j,ke)=1
        ibf(i,j,ke+1)=1
       end if
   10 continue
      do 20 i=is,ie-1
      do 20 k=ks,ke-1
       if(ibf(i,js,k).ne.0) then
        ibf(i,js-1,k)=1
        ibf(i,js-2,k)=1
       end if
       if(ibf(i,je-1,k).ne.0) then
        ibf(i,je,k)=1
        ibf(i,je+1,k)=1
       end if
   20 continue
c
      return
      end
c---------------------------------------------------------- sub.rduvwi
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.rduvwi
      subroutine rduvwi(ifile,infile)
      include 'common_vib8.fil'
C 8/27/13
      common /inflow/Vin
c
      character*14 ifile,infile
c
c
c------------------------------------
      if(ifile(1:5).ne.'fffff') then
c------------------------------------
      print *,'input file = ',ifile
      open(unit=11,file=ifile,status='old',form='unformatted')
       read(11)is2,ie2,js2,je2,ks2,ke2
       read(11)re2,cs2,dt2,t2,nt2,iype2,izpe2,ite2
       read(11)(((uu(i,j,k),i=is,ie+1),j=js,je-1),k=ks,ke-1)
       read(11)(((vv(i,j,k),i=is,ie),j=js,je),k=ks,ke-1)
       read(11)(((ww(i,j,k),i=is,ie),j=js,je-1),k=ks,ke)
       read(11)(((pp(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
c
       read(11)(((f(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       read(11)(((g(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       read(11)(((h(i,j,k),i=is,ie),j=js,je),k=ks,ke)
c
      close(11)
c
C 9/ 3/13 skip if infile='fffff'
      if(infile(1:5).ne.'fffff') then
c****************************************** add for cal. vibration
      open(unit=12,file=infile,status='old',form='unformatted')
       read(12)is3,ie3,js3,je3
       read(12)xacc,yacc,tacc
       read(12)xdis,ydis,tdis
       read(12)xvel,yvel,tvel
       read(12)(((facc1(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       read(12)(((gacc1(i,j,k),i=is,ie),j=js,je),k=ks,ke)
      close(12)
c****************************************** 
C 9/ 3/13
      endif
      iswch(1)=0
      iswch(2)=0
      iswch(3)=0
c
      print *,'t2 =',t2
      t=t2
c
c
c------------------------------------
      else
c------------------------------------
      do 10 k=ks,ke
      do 10 j=js,je
      do 10 i=is,ie
C 8/27/13
C       u(i,j,k)=1.0d0
       u(i,j,k)=Vin
       v(i,j,k)=0.0d0
       w(i,j,k)=0.0d0
       p(i,j,k)=0.0d0
   10 continue
c
      do 20 k=ks,ke
      do 20 j=js,je
C 8/27/13
C       u(is,j,k)=1.0d0
       u(is,j,k)=Vin
   20 continue
c
      ca01=(1.0d0+ca1)/2.0d0
      print *,'ca01=',ca01
      do 30 k=ks,ke
      do 30 j=js,je
      do 30 i=is,ie
       if(ibf(i,j,k).ne.0.and.ibf(i-1,j,k).eq.0) then
        u(i,j,k)=0.0d0
        u(i,j+1,k)=0.0d0
        u(i,j,k+1)=0.0d0
        u(i,j+1,k+1)=0.0d0
       end if
       if(ibf(i,j,k).ne.0.and.ibf(i+1,j,k).eq.0) then
        u(i+1,j,k)=0.0d0
        u(i+1,j+1,k)=0.0d0
        u(i+1,j,k+1)=0.0d0
        u(i+1,j+1,k+1)=0.0d0
       end if
c
       if(ibf(i,j,k).ne.0.and.ibf(i,j-1,k).eq.0) then
        u(i,j,k)=ca01
        u(i+1,j,k)=ca01
        u(i,j,k+1)=ca01
        u(i+1,j,k+1)=ca01
       end if
       if(ibf(i,j,k).ne.0.and.ibf(i,j+1,k).eq.0) then
        u(i,j+1,k)=ca01
        u(i+1,j+1,k)=ca01
        u(i,j+1,k+1)=ca01
        u(i+1,j+1,k+1)=ca01
       end if
c
       if(ibf(i,j,k).ne.0.and.ibf(i,j,k-1).eq.0) then
        u(i,j,k)=ca01
        u(i+1,j,k)=ca01
        u(i,j+1,k)=ca01
        u(i+1,j+1,k)=ca01
       end if
       if(ibf(i,j,k).ne.0.and.ibf(i,j,k+1).eq.0) then
        u(i,j,k+1)=ca01
        u(i+1,j,k+1)=ca01
        u(i,j+1,k+1)=ca01
        u(i+1,j+1,k+1)=ca01
       end if
   30 continue
c
c
      call uvwrvs
      call prsrvs
c
      do 70 k=ks,ke-1
      do 70 j=js,je-1
       uo(j,k,2)=3.0d0*uu(ie,j,k)-3.0d0*uu(ie-1,j,k)+uu(ie-2,j,k)
       uo(j,k,1)=uu(ie,j,k)
   70 continue
      do 80 k=ks,ke-1
      do 80 j=js,je
       vo(j,k,2)=3.0d0*vv(ie-1,j,k)-3.0d0*vv(ie-2,j,k)+vv(ie-3,j,k)
       vo(j,k,1)=vv(ie-1,j,k)
   80 continue
      do 90 k=ks,ke
      do 90 j=js,je-1
       wo(j,k,2)=3.0d0*ww(ie-1,j,k)-3.0d0*ww(ie-2,j,k)+ww(ie-3,j,k)
       wo(j,k,1)=ww(ie-1,j,k)
   90 continue
c
      iswch(1)=1
      iswch(2)=1
      iswch(3)=1
c
c--------------------------------
      end if
c--------------------------------
c
      do 40 j=js,je-1
      do 40 k=ks,ke-1
       ui(j,k)=uu(is,j,k)
   40 continue
c
      do 50 j=js,je
      do 50 k=ks,ke-1
       vi(j,k)=vv(is,j,k)
   50 continue
c
      do 60 j=js,je-1
      do 60 k=ks,ke
       wi(j,k)=ww(is,j,k)
   60 continue
c
      call uvwbdy
      call prsbdy
c
c
      return
      end
c----------------------------------------------------------- sub.rdvelo
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c234567----------------------------------------------------- sub.rdvelo
      subroutine rdvelo(ifile)
      include 'common_vib8.fil'
C 8/27/13
      common /inflow/Vin
      character*14 ifile
c
c
      do 10 k=ks,ke-1
      do 10 j=js-1,je
C 8/27/13
C       uu(is,j,k)=1.0d0
C       uu(is-1,j,k)=1.0d0
C       ui(j,k)=1.0d0
       uu(is,j,k)=Vin
       uu(is-1,j,k)=Vin
       ui(j,k)=Vin
   10 continue
c
      do 20 k=ks,ke-1
      do 20 j=js-1,je+1
       vv(is,j,k)=0.0d0 
       vv(is-1,j,k)=0.0d0
       vi(j,k)=0.0d0
   20 continue
c
      do 30 k=ks,ke
      do 30 j=js-1,je
       ww(is,j,k)=0.0d0
       ww(is-1,j,k)=0.0d0
       wi(j,k)=0.0d0
   30 continue
c
C 8/31/13 remove setting pp                                  sub.rdvelo
C      do 40 k=ks,ke-1
C      do 40 j=js-1,je
C       pp(is,j,k)=0.0d0
C       pp(is-1,j,k)=0.0d0
C   40 continue
c
c************************************* add for cal. vibration
      if((ndof.eq.2).or.(ndof.eq.4).or.(ndof.eq.6)) then
      do 11 k=ks,ke-1
      do 11 j=js-1,je
       uu(is,j,k)  =uu(is,j,k)  -xvel
c       uu(is+1,j,k)=uu(is+1,j,k)-xvel
       ui(j,k)     =ui(j,k)     -xvel
   11 continue
      end if
c
      if((ndof.eq.1).or.(ndof.eq.4).or.(ndof.eq.5).or.(ndof.eq.6)) then
      do 21 k=ks,ke-1
      do 21 j=js-1,je+1
       vv(is,j,k)=vv(is,j,k)-yvel 
       vi(j,k)=   vi(j,k)   -yvel
   21 continue
      end if
c
      if((ndof.eq.3).or.(ndof.eq.5).or.(ndof.eq.6)) then
      do 31 k=ks,ke-1
      do 31 j=js-1,je
      yyy=(y(j)+y(j+1))*0.5d0
      ftheta=datan2(yyy,x(is))
      fr=dsqrt(x(is)**2+yyy**2)
       uu(is,j,k)=ui(j,k)-tvel*fr*(-dsin(ftheta+tdis)) 
       ui(j,k)     =uu(is,j,k)
c
      yyy=(y(j)+y(j+1))*0.5d0
      ftheta1=datan2(yyy,x(is-1))
      fr=dsqrt(x(is-1)**2+yyy**2)
      uu(is-1,j,k)=uu(is-1,j,k)-tvel*fr*(-dsin(ftheta+tdis)) 
   31 continue
c
c
      do 32 k=ks,ke-1
      do 32 j=js-1,je+1
      xxx=(x(is)+x(is+1))*0.5d0
      gtheta=datan2(y(j),xxx)
      gr=dsqrt(xxx**2+y(j)**2)
       vv(is,j,k)=vi(j,k)-tvel*gr*( dcos(gtheta+tdis)) 
       vi(j,k)     =vv(is,j,k)
c
      xxx=(x(is-1)+x(is))*0.5d0
      gtheta=datan2(y(j),xxx)
      gr=dsqrt(xxx**2+y(j)**2)
       vv(is-1,j,k)=vv(is-1,j,k)-tvel*gr*( dcos(gtheta+tdis)) 
   32 continue
      end if
cc
c
c
c************************************* add for relative vel.
      if(ifile(1:5).eq.'fffff'.and.nn.le.2) then
c
      print *,'#############################  ',ifile,nn
      if((ndof.eq.2).or.(ndof.eq.4).or.(ndof.eq.6)) then
      do 41 k=ks,ke-1
      do 41 j=js,je-1
      do 41 i=is+2,ie+1
       uu(i,j,k)  =uu(i,j,k)  -(xvel-xvel0)
   41 continue
      end if
c
      if((ndof.eq.1).or.(ndof.eq.4).or.(ndof.eq.5).or.(ndof.eq.6)) then
      do 51 k=ks,ke-1
      do 51 j=js,je
      do 51 i=is+1,ie
       vv(i,j,k)=vv(i,j,k)-(yvel-yvel0) 
   51 continue
      end if
c
      do 61 k=ks,ke-1
      do 61 j=js,je
      do 61 i=is+1,ie+1
      yyy=(y(j)+y(j+1))*0.5d0
      ftheta=datan2(yyy,x(i))
      fr=dsqrt(x(i)**2+yyy**2)
      uu(i,j,k)=uu(i,j,k)-(tvel-tvel0)*fr*(-dsin(ftheta+tdis))
   61 continue
c
      do 62 k=ks,ke-1
      do 62 j=js,je
      do 62 i=is+1,ie
      xxx=(x(i)+x(i+1))*0.5d0
      gtheta=datan2(y(j),xxx)
      gr=dsqrt(xxx**2+y(j)**2)
      vv(i,j,k)=vv(i,j,k)-(tvel-tvel0)*gr*(dcos(gtheta+tdis))
   62 continue
c
      end if
c
      call uvwbdy
c
      write(1,111)nn,t,uu(is,80,10),uu(is+2,80,10),uu(is+5,80,10)
      write(2,111)nn,t,vv(is,80,10),vv(is+2,80,10),vv(is+5,80,10)
  111 format(i10,4f15.10)
c
c
      return
      end
c--------------------------------------------------------------- sub.hsmac2
c                                                                         |
c                                                                         |
c                                                                         |
c                                                                         |
c                                                                         |
c                                                                         |
c                                                                         |
c                                                                         |
c234567--------------------------------------------------------- sub.hsmac2
      subroutine hsmac2(dtmin)
      include 'common_vib8.fil'
c
      iter=0
      itmax=5000
c------------------------------------------------
      esp2=esp
      open(11,file='esp2.inf',status='unknown',form='formatted')
       read(11,112)icount2,icount3,sesp
      close(11)
  112 format(2i10,f10.5)
       print *,'esp esp2=',esp,esp2
       print *,'ic2 ic3 sesp',icount2,icount3,sesp
      nn2=nn-icount2
      if(nn2.le.1) nn2=1
      if(nn2.ge.icount3) nn2=icount3
      rnc=(esp/sesp)**(1.0/real(icount3-1))
      esp2=sesp*rnc**(nn2-1) 
      print *,'nn   =',nn
      print *,'nn2  =',nn2
      print *,'rnc  =',rnc
      print *,'sesp =',sesp
      print *,'esp  =',esp
      print *,'esp2 = ',esp2
c------------------------------------------------
c
 1000 iter=iter+1
      if(iter.gt.itmax) go to 2000
      dm=0.0d0
cccccccccccccccccccccccccc
      i2=0
      j2=0
      k2=0
cccccccccccccccccccccccccc
c
      do 10 k=ks,ke-1
      do 10 j=js,je-1
      do 10 i=is,ie-1
       if(ibf(i,j,k).ne.0) go to 3000
       d=(uu(i+1,j,k)-uu(i,j,k))/dx(i+1)
     &  +(vv(i,j+1,k)-vv(i,j,k))/dy(j+1)
     &  +(ww(i,j,k+1)-ww(i,j,k))/dz(k+1)
       d1=dabs(d)
cccccccccccccccccccccccccc
       if(d1.gt.dm) then
        i2=i
        j2=j
        k2=k
       end if
cccccccccccccccccccccccccc
       dm=dmax1(d1,dm)
 3000 dmy=0.0d0
   10 continue
c     print *,'iter dm esp =',iter,dm,esp
c------------------------------------------------
cc    if(dm.lt.esp) go to 2000
      if(dm.lt.esp2) go to 2000
c
      open(11,file='iter.dat',status='unknown',form='formatted')
       write(11,111)iter,dm,esp2,icount2
cccccccccccccccccccccccccc
       write(11,113)i2,j2,k2
cccccccccccccccccccccccccc
      close(11)
  111 format(i10,2f20.10,i10)
  113 format(3i5)
c------------------------------------------------
c
c
      do 20 i=is,ie-1
      do 20 j=js,je-1
      do 20 k=ks,ke-1
       d=(uu(i+1,j,k)-uu(i,j,k))/dx(i+1)
     &  +(vv(i,j+1,k)-vv(i,j,k))/dy(j+1)
     &  +(ww(i,j,k+1)-ww(i,j,k))/dz(k+1)
       if(i.ne.ie) dxyz=1.0d0/(dxh(i)*dxh(i+1))
     &                 +1.0d0/(dyh(j)*dyh(j+1))
     &                 +1.0d0/(dzh(k)*dzh(k+1))
c
       if(ibf(i,j,k).ne.0) d=0.0d0
       d2=dabs(d)
c
       dp=-omg*d/(2.0d0*dt*dxyz)
       if(d2.ge.esp2) pp(i,j,k)=pp(i,j,k)+dp
cc
       if(i.ne.ie-1) uu(i+1,j,k)=uu(i+1,j,k)+dt*dp/dxh(i+1)
       uu(i,j,k)=uu(i,j,k)-dt*dp/dxh(i)
       vv(i,j+1,k)=vv(i,j+1,k)+dt*dp/dyh(j+1)
       vv(i,j,k)=vv(i,j,k)-dt*dp/dyh(j)
       ww(i,j,k+1)=ww(i,j,k+1)+dt*dp/dzh(k+1)
       ww(i,j,k)=ww(i,j,k)-dt*dp/dzh(k)
cc
   20 continue
c
c
      do 30 k=ks,ke-1 
      do 30 j=js,je-1
       uu(is,j,k)=ui(j,k)
C 8/31/13 remove setting pp                            sub.hsmac2
C       pp(is,j,k)=0.0d0
   30 continue
c
      do 40 k=ks,ke-1 
      do 40 j=js,je
       vv(is,j,k)=vi(j,k)
   40 continue
c
      do 50 k=ks,ke
      do 50 j=js,je-1
       ww(is,j,k)=wi(j,k)
   50 continue
c
      do 60 k=ks,ke-1
      do 60 j=js,je-1
      do 60 i=is,ie-1
       if(ibf(i,j,k).ne.0) then
       uu(i,j,k)=0.0d0
       vv(i,j,k)=0.0d0
       ww(i,j,k)=0.0d0
       uu(i+1,j,k)=0.0d0
       vv(i,j+1,k)=0.0d0
       ww(i,j,k+1)=0.0d0
       pp(i,j,k)=0.0d0
       end if
   60 continue 
c
c
c
c      if(iyper.ne.0) then
c       do 70 k=ks,ke-1
c       do 70 i=is,ie-1
c        vv(i,js,k)=0.0d0
c        vv(i,je,k)=0.0d0
c   70  continue
c      end if
      if(izper.ne.0) then
       do 80 j=js,je-1
       do 80 i=is,ie-1
        ww(i,j,ks)=0.0d0
        ww(i,j,ke)=0.0d0
   80  continue
      end if
c
      go to 1000
c
 2000 print *,'Iteration =',iter
      print *,'esp , dm  =',esp,dm
c
c
      d4=1.0d0/delm
      dtmax=0.0d0
      do 90 i=is,ie-1
      do 90 j=js,je-1
      do 90 k=ks,ke-1
       d1=uu(i,j,k)/dx(i+1)
       d2=vv(i,j,k)/dy(j+1)
       d3=ww(i,j,k)/dz(k+1)
       dtmax=dmax1(dtmax,d1)
       dtmax=dmax1(dtmax,d2)
       dtmax=dmax1(dtmax,d3)
   90 continue
       dtmax=dmax1(dtmax,d4)
      dtmin=1.0d0/dtmax
      print *,'max dt  =',dtmin
c
c
      iswch(1)=0
      iswch(2)=1
      iswch(3)=1
      call uvwbdy
      call prsbdy
c
      return 
      end
c---------------------------------------------------------- sub.smagor
c             XNU by Smagorinsky model                               | 
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.smagor
      subroutine smagor
      include 'common_vib8.fil'
c
c-------------------------------------------- for around b-body
      do 10 k=ks,ke-1
      do 10 j=js,je-1
      do 10 i=is,ie-1
       if(ibf(i,j,k).ne.0) then
        uu(i,j,k)=0.0d0
        vv(i,j,k)=0.0d0
        ww(i,j,k)=0.0d0
        uu(i+1,j,k)=0.0d0
        vv(i,j+1,k)=0.0d0
        ww(i,j,k+1)=0.0d0
       end if
   10 continue
      do 15 k=ks,ke-1
      do 15 j=js,je-1
      do 15 i=is,ie-1
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+1,k).ne.0.and.ibf(i-1,j+1,k).ne.0) 
     &   uu(i,j+1,k)=ca1*uu(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &   uu(i,j-1,k)=ca1*uu(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i-1,j,k+1).ne.0)
     &   uu(i,j,k+1)=ca1*uu(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &   uu(i,j,k-1)=ca1*uu(i,j,k)
      end if 
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j-1,k).ne.0)
     &   vv(i+1,j,k)=ca1*vv(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &   vv(i-1,j,k)=ca1*vv(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i,j-1,k+1).ne.0)
     &   vv(i,j,k+1)=ca1*vv(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &   vv(i,j,k-1)=ca1*vv(i,j,k)
c
      end if 
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j,k-1).ne.0)
     &   ww(i+1,j,k)=ca1*ww(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &   ww(i-1,j,k)=ca1*ww(i,j,k)
      if(ibf(i,j+1,k).ne.0.and.ibf(i,j+1,k-1).ne.0)
     &   ww(i,j+1,k)=ca1*ww(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &   ww(i,j-1,k)=ca1*ww(i,j,k)
      end if
c
   15 continue
c
c-------------------------------------------- for around b-body
c
      do 100 i=is,ie-1
         delx=dx(i+1)
      do 100 j=js,je-1
         dely=dy(j+1)
      do 100 k=ks,ke-1
         delz=dz(k+1)
c
       if(ibf(i,j,k).ne.0) go to 1000
c
       del=(delx*dely*delz)**0.33333
c
       ujp=0.5d0*(dy(j+2)*(uu(i,j,k)+uu(i+1,j,k))
     &         +dy(j+1)*(uu(i,j+1,k)+uu(i+1,j+1,k)))/(dy(j+1)+dy(j+2))
       ujm=0.5d0*(dy(j)*(uu(i,j,k)+uu(i+1,j,k))
     &         +dy(j+1)*(uu(i,j-1,k)+uu(i+1,j-1,k)))/(dy(j)+dy(j+1))
       ukp=0.5d0*(dz(k+2)*(uu(i,j,k)+uu(i+1,j,k))
     &         +dz(k+1)*(uu(i,j,k+1)+uu(i+1,j,k+1)))/(dz(k+1)+dz(k+2))
       ukm=0.5d0*(dz(k)*(uu(i,j,k)+uu(i+1,j,k))
     &         +dz(k+1)*(uu(i,j,k-1)+uu(i+1,j,k-1)))/(dz(k)+dz(k+1))
c
       vip=0.5d0*(dx(i+2)*(vv(i,j,k)+vv(i,j+1,k))
     &         +dx(i+1)*(vv(i+1,j,k)+vv(i+1,j+1,k)))/(dx(i+1)+dx(i+2))
       vim=0.5d0*(dx(i)*(vv(i,j,k)+vv(i,j+1,k))
     &         +dx(i+1)*(vv(i-1,j,k)+vv(i-1,j+1,k)))/(dx(i)+dx(i+1))
       vkp=0.5d0*(dz(k+2)*(vv(i,j,k)+vv(i,j+1,k))
     &         +dz(k+1)*(vv(i,j,k+1)+vv(i,j+1,k+1)))/(dz(k+1)+dz(k+2))
       vkm=0.5d0*(dz(k)*(vv(i,j,k)+vv(i,j+1,k))
     &         +dz(k+1)*(vv(i,j,k-1)+vv(i,j+1,k-1)))/(dz(k)+dz(k+1))
c
       wip=0.5d0*(dx(i+2)*(ww(i,j,k)+ww(i,j,k+1))
     &         +dx(i+1)*(ww(i+1,j,k)+ww(i+1,j,k+1)))/(dx(i+1)+dx(i+2))
       wim=0.5d0*(dx(i)*(ww(i,j,k)+ww(i,j,k+1))
     &         +dx(i+1)*(ww(i-1,j,k)+ww(i-1,j,k+1)))/(dx(i)+dx(i+1))
       wjp=0.5d0*(dy(j+2)*(ww(i,j,k)+ww(i,j,k+1))
     &         +dy(j+1)*(ww(i,j+1,k)+ww(i,j+1,k+1)))/(dy(j+1)+dy(j+2))
       wjm=0.5d0*(dy(j)*(ww(i,j,k)+ww(i,j,k+1))
     &         +dy(j+1)*(ww(i,j-1,k)+ww(i,j-1,k+1)))/(dy(j)+dy(j+1))
c
c-------------------------------------------- for around b-body
      if(ca1.ne.-1.0) go to 999
      if(ibf(i+1,j,k).ne.0) then
       vip=0.0d0
       wip=0.0d0
      end if
      if(ibf(i-1,j,k).ne.0) then
       vim=0.0d0
       wim=0.0d0
      end if
c
      if(ibf(i,j+1,k).ne.0) then
       ujp=0.0d0
       wjp=0.0d0
      end if
      if(ibf(i,j-1,k).ne.0) then
       ujm=0.0d0
       wjm=0.0d0
      end if
c
      if(ibf(i,j,k+1).ne.0) then
       ukp=0.0d0
       vkp=0.0d0
      end if
      if(ibf(i,j,k-1).ne.0) then
       ukm=0.0d0
       vkm=0.0d0
      end if
  999 dmy=0.0
c
c-------------------------------------------- for around b-body
c
       d11=(uu(i+1,j,k)-uu(i,j,k))/dx(i+1)
       d22=(vv(i,j+1,k)-vv(i,j,k))/dy(j+1)
       d33=(ww(i,j,k+1)-ww(i,j,k))/dz(k+1)
       d12=(ujp-ujm)/dy(j+1)+(vip-vim)/dx(i+1)
       d13=(ukp-ukm)/dz(k+1)+(wip-wim)/dx(i+1)
       d23=(vkp-vkm)/dz(k+1)+(wjp-wjm)/dy(j+1)
c
        tnu(i,j,k)=(cs*del)**2
     &            *sqrt(2.0d0*d11**2+2.0d0*d22**2+2.0d0*d33**2
     &            +d12**2+d13**2+d23**2) 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      coef=2.0d0/3.0d0
      ck=0.094d0
      ck2=ck*ck
       tjp=(tnu(i,j,k)*dx(i)+tnu(i-1,j,k)*dx(i+1))/(dx(i)+dx(i+1))
       tjm=(tnu(i,j-1,k)*dx(i)+tnu(i-1,j-1,k)*dx(i+1))/(dx(i)+dx(i+1))
         tkp=(tjp*dy(j)+tjm*dy(j+1))/(dy(j)+dy(j+1))
c
       tjp2=(tnu(i,j,k-1)*dx(i)+tnu(i-1,j,k-1)*dx(i+1))/(dx(i)+dx(i+1))
       tjm2=(tnu(i,j-1,k-1)*dx(i)+tnu(i-1,j-1,k-1)*dx(i+1))
     &     /(dx(i)+dx(i+1))
         tkm=(tjp2*dy(j)+tjm2*dy(j+1))/(dy(j)+dy(j+1))
c
       tijk=(tkp*dz(k)+tkm*dz(k+1))/(dz(k)+dz(k+1))
c
       rss(i,j,k,1)=rss(i,j,k,1)-tijk*d12
       rss(i,j,k,2)=rss(i,j,k,2)-tijk*d23
       rss(i,j,k,3)=rss(i,j,k,3)-tijk*d13
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
 1000 dmy=0.0d0
c
  100 continue
c
c-------------------------------------------- for around b-body
      do 40 k=ks,ke-1
      do 40 j=js,je-1
      do 40 i=is,ie-1
       if(ibf(i,j,k).ne.0) then
        tnu(i,j,k)=0.0d0
       end if
   40 continue
c-------------------------------------------- for around b-body
c
c---------------------------------- inflow outflow
      do 110 j=js,je-1
      do 110 k=ks,ke-1
       tnu(is-1,j,k)=tnu(is,j,k)
c      tnu(ie,j,k)=tnu(ie-1,j,k)
       tnu(ie,j,k)=3.0d0*tnu(ie-1,j,k)-3.0d0*tnu(ie-2,j,k)
     &            +tnu(ie-3,j,k)
  110 continue
c
c---------------------------------- z-x plane (y-dir.)
      if(iyper.eq.0) then
       do 120 i=is,ie-1
       do 120 k=ks,ke-1
        tnu(i,js-1,k)=tnu(i,je-1,k)
        tnu(i,je,k)=tnu(i,js,k)
  120  continue
      else
       do 130 i=is,ie-1
       do 130 k=ks,ke-1
        tnu(i,js-1,k)=tnu(i,js,k)
        tnu(i,je,k)=tnu(i,je-1,k)
  130  continue
      end if
c
c---------------------------------- x-y plane (z-dir.)
      if(izper.eq.0) then
       do 140 i=is,ie-1
       do 140 j=js,je-1
        tnu(i,j,ks-1)=tnu(i,j,ke-1)
        tnu(i,j,ke)=tnu(i,j,ks)
  140  continue
      else
       do 150 i=is,ie-1
       do 150 j=js,je-1
        tnu(i,j,ks-1)=tnu(i,j,ks)
        tnu(i,j,ke)=tnu(i,j,ke-1)
  150  continue
      end if
c
c----------------------------------------- outside corner
c---------------------------------- y-z plane (x-dir.)
      if(iyper.eq.0) then
       if(izper.eq.0) then
        do 160 i=is,ie-1
         tnu(i,js-1,ks-1)=tnu(i,je-1,ke-1)
         tnu(i,je,ks-1)=tnu(i,js,ke-1)
         tnu(i,js-1,ke)=tnu(i,je-1,ks)
         tnu(i,je,ke)=tnu(i,js,ks)
  160   continue
       else
        do 170 i=is,ie-1
         tnu(i,js-1,ks-1)=tnu(i,je-1,ks)
         tnu(i,je,ks-1)=tnu(i,js,ks)
         tnu(i,js-1,ke)=tnu(i,je-1,ke-1)
         tnu(i,je,ke)=tnu(i,js,ke-1)
  170   continue
       end if
c
      else
c
       if(izper.eq.0) then
        do 180 i=is,ie-1
         tnu(i,js-1,ks-1)=tnu(i,js,ke-1)
         tnu(i,je,ks-1)=tnu(i,je-1,ke-1)
         tnu(i,js-1,ke)=tnu(i,js,ks)
         tnu(i,je,ke)=tnu(i,je-1,ks)
  180   continue
       else
        do 190 i=is,ie-1
         tnu(i,js-1,ks-1)=tnu(i,js,ks)
         tnu(i,je,ks-1)=tnu(i,je-1,ks)
         tnu(i,js-1,ke)=tnu(i,js,ke-1)
         tnu(i,je,ke)=tnu(i,je-1,ke-1)
  190   continue
       end if
c
      end if
c---------------------------------- z-x plane (y-dir.)
      if(izper.eq.0) then
       do 200 j=js,je-1
        tnu(is-1,j,ks-1)=tnu(is,j,ke-1)
        tnu(ie,j,ks-1)=tnu(ie-1,j,ke-1)
        tnu(is-1,j,ke)=tnu(is,j,ks)
        tnu(ie,j,ke)=tnu(ie-1,j,ks)
  200  continue
      else
       do 210 j=js,je-1
        tnu(is-1,j,ks-1)=tnu(is,j,ks)
        tnu(ie,j,ks-1)=tnu(ie-1,j,ks)
        tnu(is-1,j,ke)=tnu(is,j,ke-1)
        tnu(ie,j,ke)=tnu(ie-1,j,ke-1)
  210  continue
      end if
c---------------------------------- x-y plane (z-dir.)
      if(iyper.eq.0) then
       do 220 k=ks,ke-1
        tnu(is-1,js-1,k)=tnu(is,je-1,k)
        tnu(ie,js-1,k)=tnu(ie-1,je-1,k)
        tnu(is-1,je,k)=tnu(is,js,k)
        tnu(ie,je,k)=tnu(ie-1,js,k)
  220  continue
      else
       do 230 k=ks,ke-1
        tnu(is-1,js-1,k)=tnu(is,js,k)
        tnu(ie,js-1,k)=tnu(ie-1,js,k)
        tnu(is-1,je,k)=tnu(is,je-1,k)
        tnu(ie,je,k)=tnu(ie-1,je-1,k)
  230  continue
      end if
c
c
      return
      end
c---------------------------------------------------------- sub.space3
c                                                                    |
c        3rd order upwind diff for convective terms                  |
c        2nd order central diff for viscous terms                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.space3
      subroutine space3
      include 'common_vib8.fil'
c
c
c
c-----------------------------------------------------------+
c          calculation of F,  3rd order upwind              |
c-----------------------------------------------------------+
      do 100 i=is,ie-1
      do 100 j=js,je-1
      do 100 k=ks,ke-1
c
      icpy=0
      icpz=0
c------------------------------------------ for around b-body
      if(ibf(i,j,k).ne.0.and.ibf(i-1,j,k).ne.0) then
       f(i,j,k)=0.0d0
       go to 1001
      end if
c
c++++++++++++++++++++++++++++++++++++
c     if(ibf(i,j,k).ne.0.and.ibf(i-1,j,k).eq.0) then
c      if(ibf(i,j-1,k).ne.0) vv(i,j,k)  =-vv(i-1,j,k)
c      if(ibf(i,j+1,k).ne.0) vv(i,j+1,k)=-vv(i-1,j+1,k)
c      if(ibf(i,j,k-1).ne.0) ww(i,j,k)  =-ww(i-1,j,k)
c      if(ibf(i,j,k+1).ne.0) ww(i,j,k+1)=-ww(i-1,j,k+1)
c     end if
c
      if(ibf(i+1,j,k).ne.0) uu(i+2,j,k)=ca1*uu(i,j,k)
      if(ibf(i-2,j,k).ne.0) uu(i-2,j,k)=ca1*uu(i,j,k)
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+1,k).ne.0.and.ibf(i-1,j+1,k).ne.0)
     &    uu(i,j+1,k)=ca1*uu(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &    uu(i,j-1,k)=ca1*uu(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i-1,j,k+1).ne.0)
     &    uu(i,j,k+1)=ca1*uu(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &    uu(i,j,k-1)=ca1*uu(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+2,k).ne.0.and.ibf(i-1,j+2,k).ne.0)
     &    uu(i,j+2,k)=ca1*uu(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i-1,j-2,k).ne.0)
     &    uu(i,j-2,k)=ca1*uu(i,j-1,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i-1,j,k+2).ne.0)
     &    uu(i,j,k+2)=ca1*uu(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i-1,j,k-2).ne.0)
     &    uu(i,j,k-2)=ca1*uu(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j-1,k).ne.0)
     &    vv(i+1,j,k)=ca1*vv(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &    vv(i-1,j,k)=ca1*vv(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i,j-1,k+1).ne.0)
     &    vv(i,j,k+1)=ca1*vv(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &    vv(i,j,k-1)=ca1*vv(i,j,k)
      end if
c 
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j-1,k).ne.0)
     &    vv(i+2,j,k)=ca1*vv(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j-1,k).ne.0)
     &    vv(i-2,j,k)=ca1*vv(i-1,j,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i,j-1,k+2).ne.0)
     &    vv(i,j,k+2)=ca1*vv(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i,j-1,k-2).ne.0)
     &    vv(i,j,k-2)=ca1*vv(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j,k-1).ne.0)
     &    ww(i+1,j,k)=ca1*ww(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &    ww(i-1,j,k)=ca1*ww(i,j,k)
      if(ibf(i,j+1,k).ne.0.and.ibf(i,j+1,k-1).ne.0)
     &    ww(i,j+1,k)=ca1*ww(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &    ww(i,j-1,k)=ca1*ww(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j,k-1).ne.0)
     &    ww(i+2,j,k)=ca1*ww(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j,k-1).ne.0)
     &    ww(i-2,j,k)=ca1*ww(i-1,j,k)
      if(ibf(i,j+2,k).ne.0.and.ibf(i,j+2,k-1).ne.0)
     &    ww(i,j+2,k)=ca1*ww(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i,j-2,k-1).ne.0)
     &    ww(i,j-2,k)=ca1*ww(i,j-1,k)
      end if
c++++++++++++++++++++++++++++++++++++
c
c
      if(ibf(i,j+1,k).ne.0) icpy=1
      if(ibf(i,j-1,k).ne.0) icpy=2
      if(ibf(i,j,k+1).ne.0) icpz=1
      if(ibf(i,j,k-1).ne.0) icpz=2
c
c------------------------------------------ for around b-body
c
      uijk=uu(i,j,k)
      uipjk=uu(i+1,j,k)
      uimjk=uu(i-1,j,k)
      uijpk=uu(i,j+1,k)
      uijmk=uu(i,j-1,k)
      uijkp=uu(i,j,k+1)
      uijkm=uu(i,j,k-1)
c
c --------------------- Convective terms in F -----------------------
      vu=0.5d0*((dx(i)*(vv(i,j,k)+vv(i,j+1,k))
     &        +dx(i+1)*(vv(i-1,j,k)+vv(i-1,j+1,k)))/(dx(i)+dx(i+1)))
      wu=0.5d0*((dx(i)*(ww(i,j,k)+ww(i,j,k+1))
     &        +dx(i+1)*(ww(i-1,j,k)+ww(i-1,j,k+1)))/(dx(i)+dx(i+1)))
c 
c-------------------------------------------------  Udu/dx
      dudxp=(uipjk-uijk)/dx(i+1)
      dudxm=(uijk-uimjk)/dx(i)
c
c ---------------------- ududx, U>=0
      if(uijk.ge.0.) then
       if(i.eq.is) then
        ududx=uijk*dudxm
       else
        dx21=x(i-1)-x(i-2)
        dx24=x(i-1)-x(i+1)
        dx31=x(i)-x(i-2)
        dx32=x(i)-x(i-1)
        dx34=x(i)-x(i+1)
        dx41=x(i+1)-x(i-2)
         xa=-dx32*dx34/(dx21*dx31*dx41)
         xb=-dx31*dx34/(dx21*dx32*dx24)
         xc=(1.0d0/dx31+1.0d0/dx32+1.0d0/dx34)
         xd=dx31*dx32/(dx41*dx24*dx34)
          ududx=uijk*(xa*uu(i-2,j,k)+xb*uimjk+xc*uijk+xd*uipjk)
       endif
      endif
c
c ---------------------- ududx, U<0
      if(uijk.lt.0.) then
       if(i.eq.ie-1) then
        ududx=uijk*dudxp
       else
        dx24=x(i-1)-x(i+1)
        dx25=x(i-1)-x(i+2)
        dx32=x(i)-x(i-1)
        dx34=x(i)-x(i+1)
        dx35=x(i)-x(i+2)
        dx45=x(i+1)-x(i+2)
         xe=-dx34*dx35/(dx32*dx24*dx25)
         xf=(1.0d0/dx32+1.0d0/dx34+1.0d0/dx35)
         xg=dx32*dx35/(dx24*dx34*dx45)
         xh=-dx32*dx34/(dx25*dx35*dx45)
          ududx=uijk*(xe*uimjk+xf*uijk+xg*uipjk+xh*uu(i+2,j,k))
       endif
      endif
c
c-------------------------------------------------  Vdu/dy
      dudyp=(uijpk-uijk)/dyh(j+1)
      dudym=(uijk-uijmk)/dyh(j)
c
c ---------------------- vdudy, V>=0
      if(vu.ge.0.) then
       if((iyper.ne.0.and.j.eq.js).or.(icpy.eq.2)) then
        vdudy=vu*dudym
       else 
        dy21=0.5d0*(y(j)-y(j-2))
        dy24=0.5d0*(y(j-1)+y(j)-y(j+1)-y(j+2))
        dy31=0.5d0*(y(j)+y(j+1)-y(j-2)-y(j-1))
        dy32=0.5d0*(y(j+1)-y(j-1))
        dy34=0.5d0*(y(j)-y(j+2))
        dy41=0.5d0*(y(j+1)+y(j+2)-y(j-2)-y(j-1))
         ya=-dy32*dy34/(dy21*dy31*dy41)
         yb=-dy31*dy34/(dy21*dy32*dy24)
         yc=(1.0d0/dy31+1.0d0/dy32+1.0d0/dy34)
         yd=dy31*dy32/(dy41*dy24*dy34)
          vdudy=vu*(ya*uu(i,j-2,k)+yb*uijmk+yc*uijk+yd*uijpk)
       endif
      endif
c ---------------------- vdudy, V<0
      if(vu.lt.0.) then
       if((iyper.ne.0.and.j.eq.je-1).or.(icpy.eq.1)) then
        vdudy=vu*dudyp
       else
        dy24=0.5d0*(y(j-1)+y(j)-y(j+1)-y(j+2))
        dy25=0.5d0*(y(j-1)+y(j)-y(j+2)-y(j+3))
        dy32=0.5d0*(y(j+1)-y(j-1))
        dy34=0.5d0*(y(j)-y(j+2))
        dy35=0.5d0*(y(j)+y(j+1)-y(j+2)-y(j+3))
        dy45=0.5d0*(y(j+1)-y(j+3))
         ye=-dy34*dy35/(dy32*dy24*dy25)
         yf=(1.0d0/dy32+1.0d0/dy34+1.0d0/dy35)
         yg=dy32*dy35/(dy24*dy34*dy45)
         yh=-dy32*dy34/(dy25*dy35*dy45)
          vdudy=vu*(ye*uijmk+yf*uijk+yg*uijpk+yh*uu(i,j+2,k))
       endif
      endif
c
c-------------------------------------------------  Wdu/dy
      dudzp=(uijkp-uijk)/dzh(k+1)
      dudzm=(uijk-uijkm)/dzh(k)
c
c ---------------------- wdudz, W>=0
      if(wu.ge.0.) then
       if((izper.ne.0.and.k.eq.ks).or.(icpz.eq.2)) then
        wdudz=wu*dudzm
       else
        dz21=0.5d0*(z(k)-z(k-2))
        dz24=0.5d0*(z(k-1)+z(k)-z(k+1)-z(k+2))
        dz31=0.5d0*(z(k)+z(k+1)-z(k-2)-z(k-1))
        dz32=0.5d0*(z(k+1)-z(k-1))
        dz34=0.5d0*(z(k)-z(k+2))
        dz41=0.5d0*(z(k+1)+z(k+2)-z(k-2)-z(k-1))
         za=-dz32*dz34/(dz21*dz31*dz41)
         zb=-dz31*dz34/(dz21*dz32*dz24)
         zc=(1.0d0/dz31+1.0d0/dz32+1.0d0/dz34)
         zd=dz31*dz32/(dz41*dz24*dz34)
          wdudz=wu*(za*uu(i,j,k-2)+zb*uijkm+zc*uijk+zd*uijkp)
       endif
      endif
c ---------------------- wdudz, W<0
      if(wu.lt.0.) then
       if((izper.ne.0.and.k.eq.ke-1).or.(icpz.eq.1)) then    
        wdudz=wu*dudzp
       else
        dz24=0.5d0*(z(k-1)+z(k)-z(k+1)-z(k+2))
        dz25=0.5d0*(z(k-1)+z(k)-z(k+2)-z(k+3))
        dz32=0.5d0*(z(k+1)-z(k-1))
        dz34=0.5d0*(z(k)-z(k+2))
        dz35=0.5d0*(z(k)+z(k+1)-z(k+2)-z(k+3))
        dz45=0.5d0*(z(k+1)-z(k+3))
         ze=-dz34*dz35/(dz32*dz24*dz25)
         zf=(1.0d0/dz32+1.0d0/dz34+1.0d0/dz35)
         zg=dz32*dz35/(dz24*dz34*dz45)
         zh=-dz32*dz34/(dz25*dz35*dz45)
          wdudz=wu*(ze*uijkm+zf*uijk+zg*uijkp+zh*uu(i,j,k+2))
       endif
      endif
c
c
c --------------    Viscous and turbulence terms in F 
      dudxx=(dudxp-dudxm)/dxh(i)
      dudyy=(dudyp-dudym)/((dyh(j)+dyh(j+1))/2.0d0)
      dudzz=(dudzp-dudzm)/((dzh(k)+dzh(k+1))/2.0d0)
c------------------ Sij for U 
      dudx=0.5d0*(dudxm+dudxp)
      dudy=0.5d0*(dudyp+dudym)
      dudz=0.5d0*(dudzp+dudzm)
c
      vipjk=(vv(i  ,j,k)+vv(i  ,j+1,k))/2.0d0
      vimjk=(vv(i-1,j,k)+vv(i-1,j+1,k))/2.0d0
      dvdx=(vipjk-vimjk)/dxh(i)
c
      wipjk=(ww(i  ,j,k)+ww(i  ,j,k+1))/2.0d0
      wimjk=(ww(i-1,j,k)+ww(i-1,j,k+1))/2.0d0
      dwdx=(wipjk-wimjk)/dxh(i)
c
c*******************************************
      if(ibf(i,j,k).ne.0.and.ibf(i-1,j,k).eq.0)
     &   tnu(i,j,k)=tnu(i-1,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0)
     &   tnu(i-1,j,k)=tnu(i,j,k)
c*******************************************
c
      tnuijk=(tnu(i-1,j,k)*dx(i+1)+tnu(i,j,k)*dx(i))/(dxh(i)*2.0d0)
      dndx=(tnu(i,j,k)-tnu(i-1,j,k))/dxh(i)
c
      if(ibf(i,j,k).eq.0.and.ibf(i,j+1,k).ne.0) tnu(i,j+1,k)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i,j-1,k).ne.0) tnu(i,j-1,k)=tnu(i,j,k)
      if(ibf(i-1,j,k).eq.0.and.ibf(i-1,j+1,k).ne.0)
     &   tnu(i-1,j+1,k)=tnu(i-1,j,k)
      if(ibf(i-1,j,k).eq.0.and.ibf(i-1,j-1,k).ne.0)
     &   tnu(i-1,j-1,k)=tnu(i-1,j,k)
c
      tnijpk=(tnu(i-1,j+1,k)*dx(i+1)+tnu(i,j+1,k)*dx(i))
     &      /(dxh(i)*2.0d0)
      tnijmk=(tnu(i-1,j-1,k)*dx(i+1)+tnu(i,j-1,k)*dx(i))
     &      /(dxh(i)*2.0d0)
      if(ibf(i-1,j+1,k).ne.0.and.ibf(i,j+1,k).eq.0)
     &  tnijpk=tnu(i,j+1,k)
      if(ibf(i-1,j+1,k).eq.0.and.ibf(i,j+1,k).ne.0)
     &  tnijpk=tnu(i-1,j+1,k)
      if(ibf(i-1,j-1,k).ne.0.and.ibf(i,j-1,k).eq.0)
     &  tnijmk=tnu(i,j-1,k)
      if(ibf(i-1,j-1,k).eq.0.and.ibf(i,j-1,k).ne.0)
     &  tnijmk=tnu(i-1,j-1,k)
c
      if(ibf(i,j+1,k).ne.0.and.ibf(i-1,j+1,k).ne.0) then
       if(ibf(i-1,j,k).eq.0.and.ibf(i,j,k).ne.0) tnijpk=tnu(i-1,j,k)
       if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0) tnijpk=tnu(i,j,k)
      end if
      if(ibf(i,j-1,k).ne.0.and.ibf(i-1,j-1,k).ne.0) then
       if(ibf(i-1,j,k).eq.0.and.ibf(i,j,k).ne.0) tnijmk=tnu(i-1,j,k)
       if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0) tnijmk=tnu(i,j,k)
      end if
      dndy=(tnijpk-tnijmk)/(dy(j+1)+(dy(j)+dy(j+2))/2.0d0)
c
      if(ibf(i,j,k).eq.0.and.ibf(i,j,k+1).ne.0) tnu(i,j,k+1)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i,j,k-1).ne.0) tnu(i,j,k-1)=tnu(i,j,k)
      if(ibf(i-1,j,k).eq.0.and.ibf(i-1,j,k+1).ne.0)
     &   tnu(i-1,j,k+1)=tnu(i-1,j,k)
      if(ibf(i-1,j,k).eq.0.and.ibf(i-1,j,k-1).ne.0)
     &   tnu(i-1,j,k-1)=tnu(i-1,j,k)
c
      tnijkp=(tnu(i-1,j,k+1)*dx(i+1)+tnu(i,j,k+1)*dx(i))
     &      /(dxh(i)*2.0d0)
      tnijkm=(tnu(i-1,j,k-1)*dx(i+1)+tnu(i,j,k-1)*dx(i))
     &      /(dxh(i)*2.0d0)
      if(ibf(i-1,j,k+1).ne.0.and.ibf(i,j,k+1).eq.0)
     &  tnijkp=tnu(i,j,k+1)
      if(ibf(i-1,j,k+1).eq.0.and.ibf(i,j,k+1).ne.0)
     &  tnijkp=tnu(i-1,j,k+1)
      if(ibf(i-1,j,k-1).ne.0.and.ibf(i,j,k-1).eq.0)
     &  tnijkm=tnu(i,j,k-1)
      if(ibf(i-1,j,k-1).eq.0.and.ibf(i,j,k-1).ne.0)
     &  tnijkm=tnu(i-1,j,k-1)
c
      if(ibf(i,j,k+1).ne.0.and.ibf(i-1,j,k+1).ne.0) then
       if(ibf(i-1,j,k).eq.0.and.ibf(i,j,k).ne.0) tnijkp=tnu(i-1,j,k)
       if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0) tnijkp=tnu(i,j,k)
      end if
      if(ibf(i,j,k-1).ne.0.and.ibf(i-1,j,k-1).ne.0) then
       if(ibf(i-1,j,k).eq.0.and.ibf(i,j,k).ne.0) tnijkm=tnu(i-1,j,k)
       if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0) tnijkm=tnu(i,j,k)
      end if
      dndz=(tnijkp-tnijkm)/(dz(k+1)+(dz(k)+dz(k+2))/2.0d0)
c
      f(i,j,k)=-(ududx+vdudy+wdudz)
     &        +(1.0d0/re+tnuijk)*(dudxx+dudyy+dudzz)
     &        +2.0d0*dndx*dudx+dndy*(dudy+dvdx)+dndz*(dudz+dwdx)
c
c
c
 1001 dmy=0.0d0
c
  100 continue
c
c
c-----------------------------------------------------------+
c          calculation of G, 3rd order upwind               |
c-----------------------------------------------------------+
      do 200 i=is,ie-1
      do 200 j=js,je
      do 200 k=ks,ke-1
c
      jcpy=0
      jcp=0
      jcpz=0
      jcpx=0
c------------------------------------------ for around b-body
      if(ibf(i,j,k).ne.0.and.ibf(i,j-1,k).ne.0) then
       g(i,j,k)=0.0d0
       go to 2001
      end if
c
c++++++++++++++++++++++++++++++++++++
      if(ibf(i,j+1,k).ne.0) vv(i,j+2,k)=ca1*vv(i,j,k)
      if(ibf(i,j-2,k).ne.0) vv(i,j-2,k)=ca1*vv(i,j,k)
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+1,k).ne.0.and.ibf(i-1,j+1,k).ne.0)
     &   uu(i,j+1,k)=ca1*uu(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &   uu(i,j-1,k)=ca1*uu(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i-1,j,k+1).ne.0)
     &   uu(i,j,k+1)=ca1*uu(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i-1,j,k-1).ne.0) 
     &   uu(i,j,k-1)=ca1*uu(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+2,k).ne.0.and.ibf(i-1,j+2,k).ne.0)
     &   uu(i,j+2,k)=ca1*uu(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i-1,j-2,k).ne.0)
     &   uu(i,j-2,k)=ca1*uu(i,j-1,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i-1,j,k+2).ne.0)
     &   uu(i,j,k+2)=ca1*uu(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i-1,j,k-2).ne.0)
     &   uu(i,j,k-2)=ca1*uu(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j-1,k).ne.0)
     &   vv(i+1,j,k)=ca1*vv(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &   vv(i-1,j,k)=ca1*vv(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i,j-1,k+1).ne.0)
     &   vv(i,j,k+1)=ca1*vv(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &   vv(i,j,k-1)=ca1*vv(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j-1,k).ne.0)
     &   vv(i+2,j,k)=ca1*vv(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j-1,k).ne.0)
     &   vv(i-2,j,k)=ca1*vv(i-1,j,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i,j-1,k+2).ne.0)
     &   vv(i,j,k+2)=ca1*vv(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i,j-1,k-2).ne.0)
     &   vv(i,j,k-2)=ca1*vv(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j,k-1).ne.0)
     &   ww(i+1,j,k)=ca1*ww(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &   ww(i-1,j,k)=ca1*ww(i,j,k)
      if(ibf(i,j+1,k).ne.0.and.ibf(i,j+1,k-1).ne.0)
     &   ww(i,j+1,k)=ca1*ww(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &   ww(i,j-1,k)=ca1*ww(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j,k-1).ne.0)
     &   ww(i+2,j,k)=ca1*ww(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j,k-1).ne.0)
     &   ww(i-2,j,k)=ca1*ww(i-1,j,k)
      if(ibf(i,j+2,k).ne.0.and.ibf(i,j+2,k-1).ne.0)
     &   ww(i,j+2,k)=ca1*ww(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i,j-2,k-1).ne.0)
     &   ww(i,j-2,k)=ca1*ww(i,j-1,k)
      end if
c++++++++++++++++++++++++++++++++++++
c
      if(ibf(i,j,k+1).ne.0) jcpz=1
      if(ibf(i,j,k-1).ne.0) jcpz=2
      if(ibf(i+1,j,k).ne.0) jcpx=1
      if(ibf(i-1,j,k).ne.0) jcpx=2
c
c------------------------------------------ for around b-body
c
      vijk=vv(i,j,k)
      vipjk=vv(i+1,j,k)
      vimjk=vv(i-1,j,k)
      vijpk=vv(i,j+1,k)
      vijmk=vv(i,j-1,k)
      vijkp=vv(i,j,k+1)
      vijkm=vv(i,j,k-1)
c
c --------------------- Convective terms in G -----------------------
      uv=0.5d0*((dy(j)*(uu(i,j,k)+uu(i+1,j,k))
     &        +dy(j+1)*(uu(i,j-1,k)+uu(i+1,j-1,k)))/(dy(j)+dy(j+1)))
      wv=0.5d0*((dy(j)*(ww(i,j,k)+ww(i,j,k+1))
     &        +dy(j+1)*(ww(i,j-1,k)+ww(i,j-1,k+1)))/(dy(j)+dy(j+1)))
c
c-------------------------------------------------  Udv/dx
      dvdxp=(vipjk-vijk)/dxh(i+1)
      dvdxm=(vijk-vimjk)/dxh(i)
c
c ---------------------- Udvdx, U>=0
      if(uv.ge.0.0) then
       if((i.eq.is).or.(jcpx.eq.2)) then
        udvdx=uv*dvdxm
       else
        dx21=0.5d0*(x(i)-x(i-2))
        dx24=0.5d0*(x(i-1)+x(i)-x(i+1)-x(i+2))
        dx31=0.5d0*(x(i)+x(i+1)-x(i-2)-x(i-1))
        dx32=0.5d0*(x(i+1)-x(i-1))
        dx34=0.5d0*(x(i)-x(i+2))
        dx41=0.5d0*(x(i+1)+x(i+2)-x(i-2)-x(i-1))
         xa=-dx32*dx34/(dx21*dx31*dx41)
         xb=-dx31*dx34/(dx21*dx32*dx24)
         xc=(1.0d0/dx31+1.0d0/dx32+1.0d0/dx34)
         xd=dx31*dx32/(dx41*dx24*dx34)
          udvdx=uv*(xa*vv(i-2,j,k)+xb*vimjk+xc*vijk+xd*vipjk)
       endif
      endif
c
c ---------------------- Udvdx, U<0
      if(uv.lt.0.0) then
       if((i.eq.ie-1).or.(jcpx.eq.1)) then
        udvdx=uv*dvdxp
       else
        dx24=0.5d0*(x(i-1)+x(i)-x(i+1)-x(i+2))
        dx25=0.5d0*(x(i-1)+x(i)-x(i+2)-x(i+3))
        dx32=0.5d0*(x(i+1)-x(i-1))
        dx34=0.5d0*(x(i)-x(i+2))
        dx35=0.5d0*(x(i)+x(i+1)-x(i+2)-x(i+3))
        dx45=0.5d0*(x(i+1)-x(i+3))
         xe=-dx34*dx35/(dx32*dx24*dx25)
         xf=(1.0d0/dx32+1.0d0/dx34+1.0d0/dx35)
         xg=dx32*dx35/(dx24*dx34*dx45)
         xh=-dx32*dx34/(dx25*dx35*dx45)
          udvdx=uv*(xe*vimjk+xf*vijk+xg*vipjk+xh*vv(i+2,j,k))
       endif
      endif
c
c-------------------------------------------------  Vdv/dy
      dvdyp=(vijpk-vijk)/dy(j+1)
      dvdym=(vijk-vijmk)/dy(j)
c
c ---------------------- Vdvdy, V>=0
      if(vijk.ge.0.) then
       if(iyper.ne.0.and.j.eq.js) then
        vdvdy=vijk*dvdym
       else
        dy21=y(j-1)-y(j-2)
        dy24=y(j-1)-y(j+1)
        dy31=y(j)-y(j-2)
        dy32=y(j)-y(j-1)
        dy34=y(j)-y(j+1)
        dy41=y(j+1)-y(j-2)
         ya=-dy32*dy34/(dy21*dy31*dy41)
         yb=-dy31*dy34/(dy21*dy32*dy24)
         yc=(1.0d0/dy31+1.0d0/dy32+1.0d0/dy34)
         yd=dy31*dy32/(dy41*dy24*dy34)
          vdvdy=vijk*(ya*vv(i,j-2,k)+yb*vijmk+yc*vijk+yd*vijpk)
       endif
      endif
c
c ---------------------- Vdvdy, V<0
      if(vijk.lt.0.) then
       if(iyper.ne.0.and.j.eq.je) then
        vdvdy=vijk*dvdyp
       else
        dy24=y(j-1)-y(j+1)
        dy25=y(j-1)-y(j+2)
        dy32=y(j)-y(j-1)
        dy34=y(j)-y(j+1)
        dy35=y(j)-y(j+2)
        dy45=y(j+1)-y(j+2)
         ye=-dy34*dy35/(dy32*dy24*dy25)
         yf=(1.0d0/dy32+1.0d0/dy34+1.0d0/dy35)
         yg=dy32*dy35/(dy24*dy34*dy45)
         yh=-dy32*dy34/(dy25*dy35*dy45)
          vdvdy=vijk*(ye*vijmk+yf*vijk+yg*vijpk+yh*vv(i,j+2,k))
        endif
      endif
c
c-------------------------------------------------  Wdv/dz
      dvdzp=(vijkp-vijk)/dzh(k+1)
      dvdzm=(vijk-vijkm)/dzh(k)
c
c ---------------------- wdvdz, W>=0
      if(wv.ge.0.0) then
       if((izper.ne.0.and.k.eq.ks).or.(jcpz.eq.2)) then
        wdvdz=wv*dvdzm
       else
        dz21=0.5d0*(z(k)-z(k-2))
        dz24=0.5d0*(z(k-1)+z(k)-z(k+1)-z(k+2))
        dz31=0.5d0*(z(k)+z(k+1)-z(k-2)-z(k-1))
        dz32=0.5d0*(z(k+1)-z(k-1))
        dz34=0.5d0*(z(k)-z(k+2))
        dz41=0.5d0*(z(k+1)+z(k+2)-z(k-2)-z(k-1))
         za=-dz32*dz34/(dz21*dz31*dz41)
         zb=-dz31*dz34/(dz21*dz32*dz24)
         zc=(1.0d0/dz31+1.0d0/dz32+1.0d0/dz34)
         zd=dz31*dz32/(dz41*dz24*dz34)
          wdvdz=wv*(za*vv(i,j,k-2)+zb*vijkm+zc*vijk+zd*vijkp)
       endif 
      endif
c
c ---------------------- wdvdz, W<0
      if(wv.lt.0.0) then
       if((izper.ne.0.and.k.eq.ke-1).or.(jcpz.eq.1)) then
        wdvdz=wv*dvdzp
       else
        dz24=0.5d0*(z(k-1)+z(k)-z(k+1)-z(k+2))
        dz25=0.5d0*(z(k-1)+z(k)-z(k+2)-z(k+3))
        dz32=0.5d0*(z(k+1)-z(k-1))
        dz34=0.5d0*(z(k)-z(k+2))
        dz35=0.5d0*(z(k)+z(k+1)-z(k+2)-z(k+3))
        dz45=0.5d0*(z(k+1)-z(k+3))
         ze=-dz34*dz35/(dz32*dz24*dz25)
         zf=(1.0d0/dz32+1.0d0/dz34+1.0d0/dz35)
         zg=dz32*dz35/(dz24*dz34*dz45)
         zh=-dz32*dz34/(dz25*dz35*dz45)
          wdvdz=wv*(ze*vijkm+zf*vijk+zg*vijkp+zh*vv(i,j,k+2))
       endif
      endif
c
c
c --------------    Viscous and turbulence terms in G
      dvdxx=(dvdxp-dvdxm)/((dxh(i)+dxh(i+1))/2.0d0)
      dvdyy=(dvdyp-dvdym)/dyh(j)
      dvdzz=(dvdzp-dvdzm)/((dzh(k)+dzh(k+1))/2.0d0)
c------------------ Sij for U 
      dvdx=0.5d0*(dvdxp+dvdxm)
      dvdy=0.5d0*(dvdyp+dvdym)
      dvdz=0.5d0*(dvdzp+dvdzm)
c
      uijpk=(uu(i,j  ,k)+uu(i+1,j  ,k))/2.0d0
      uijmk=(uu(i,j-1,k)+uu(i+1,j-1,k))/2.0d0
      dudy=(uijpk-uijmk)/dyh(j)
c
      wijpk=(ww(i,j  ,k)+ww(i,j  ,k+1))/2.0d0
      wijmk=(ww(i,j-1,k)+ww(i,j-1,k+1))/2.0d0
      dwdy=(wijpk-wijmk)/dyh(j)
c
c
      if(ibf(i,j,k).ne.0.and.ibf(i,j-1,k).eq.0)
     &   tnu(i,j,k)=tnu(i,j-1,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0)
     &   tnu(i,j-1,k)=tnu(i,j,k)
c
      tnuijk=(tnu(i,j-1,k)*dy(j+1)+tnu(i,j,k)*dy(j))/(dyh(j)*2.0d0)
c
      if(ibf(i,j,k).eq.0.and.ibf(i+1,j,k).ne.0) tnu(i+1,j,k)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i-1,j,k).ne.0) tnu(i-1,j,k)=tnu(i,j,k)
      if(ibf(i,j-1,k).eq.0.and.ibf(i+1,j-1,k).ne.0)
     &   tnu(i+1,j-1,k)=tnu(i,j-1,k)
      if(ibf(i,j-1,k).eq.0.and.ibf(i-1,j-1,k).ne.0)
     &   tnu(i-1,j-1,k)=tnu(i,j-1,k)
c
      tnipjk=(tnu(i+1,j-1,k)*dy(j+1)+tnu(i+1,j,k)*dy(j))
     &      /(dyh(j)*2.0d0)
      tnimjk=(tnu(i-1,j-1,k)*dy(j+1)+tnu(i-1,j,k)*dy(j))
     &      /(dyh(j)*2.0d0)
      if(ibf(i+1,j-1,k).ne.0.and.ibf(i+1,j,k).eq.0)
     &  tnipjk=tnu(i+1,j,k)
      if(ibf(i+1,j-1,k).eq.0.and.ibf(i+1,j,k).ne.0)
     &  tnipjk=tnu(i+1,j-1,k)
      if(ibf(i-1,j-1,k).ne.0.and.ibf(i-1,j,k).eq.0)
     &  tnimjk=tnu(i-1,j,k)
      if(ibf(i-1,j-1,k).eq.0.and.ibf(i-1,j,k).ne.0)
     &  tnimjk=tnu(i-1,j-1,k)
c
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j-1,k).ne.0) then
       if(ibf(i,j-1,k).eq.0.and.ibf(i,j,k).ne.0) tnipjk=tnu(i,j-1,k)
       if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0) tnipjk=tnu(i,j,k)
      end if
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j-1,k).ne.0) then
       if(ibf(i,j-1,k).eq.0.and.ibf(i,j,k).ne.0) tnimjk=tnu(i,j-1,k)
       if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0) tnimjk=tnu(i,j,k)
      end if
      dndx=(tnipjk-tnimjk)/(dx(i+1)+(dx(i)+dx(i+2))/2.0d0)
c
      dndy=(tnu(i,j,k)-tnu(i,j-1,k))/dyh(j)
c
      if(ibf(i,j,k).eq.0.and.ibf(i,j,k+1).ne.0) tnu(i,j,k+1)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i,j,k-1).ne.0) tnu(i,j,k-1)=tnu(i,j,k)
      if(ibf(i,j-1,k).eq.0.and.ibf(i,j-1,k+1).ne.0)
     &   tnu(i,j-1,k+1)=tnu(i,j-1,k)
      if(ibf(i,j-1,k).eq.0.and.ibf(i,j-1,k-1).ne.0)
     &   tnu(i,j-1,k-1)=tnu(i,j-1,k)
c
      tnijkp=(tnu(i,j-1,k+1)*dy(j+1)+tnu(i,j,k+1)*dy(j))
     &      /(dyh(j)*2.0d0)
      tnijkm=(tnu(i,j-1,k-1)*dy(j+1)+tnu(i,j,k-1)*dy(j))
     &      /(dyh(j)*2.0d0)
      if(ibf(i,j-1,k+1).ne.0.and.ibf(i,j,k+1).eq.0)
     &  tnijkp=tnu(i,j,k+1)
      if(ibf(i,j-1,k+1).eq.0.and.ibf(i,j,k+1).ne.0)
     &  tnijkp=tnu(i,j-1,k+1)
      if(ibf(i,j-1,k-1).ne.0.and.ibf(i,j,k-1).eq.0)
     &  tnijkm=tnu(i,j,k-1)
      if(ibf(i,j-1,k-1).eq.0.and.ibf(i,j,k-1).ne.0)
     &  tnijkm=tnu(i,j-1,k-1)
c
      if(ibf(i,j,k+1).ne.0.and.ibf(i,j-1,k+1).ne.0) then
       if(ibf(i,j-1,k).eq.0.and.ibf(i,j,k).ne.0) tnijkp=tnu(i,j-1,k)
       if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0) tnijkp=tnu(i,j,k)
      end if
      if(ibf(i,j,k-1).ne.0.and.ibf(i,j-1,k-1).ne.0) then
       if(ibf(i,j-1,k).eq.0.and.ibf(i,j,k).ne.0) tnijkm=tnu(i,j-1,k)
       if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0) tnijkm=tnu(i,j,k)
      end if
c
      dndz=(tnijkp-tnijkm)/(dz(k+1)+(dz(k)+dz(k+2))/2.0d0)
c
      g(i,j,k)=-(udvdx+vdvdy+wdvdz)
     &        +(1.0d0/re+tnuijk)*(dvdxx+dvdyy+dvdzz)
     &        +dndx*(dvdx+dudy)+2.0d0*dndy*dvdy+dndz*(dvdz+dwdy)
c
c
c
 2001 dmy=0.0d0
c
  200   continue
c
c
c-----------------------------------------------------------+
c          calculation of H,  3rd order upwind              |
c-----------------------------------------------------------+
      do 300 i=is,ie-1
      do 300 j=js,je-1
      do 300 k=ks,ke
c 
      kcpz=0
      kcp=0
      kcpx=0
      kcpy=0
c------------------------------------------ for around b-body
      if(ibf(i,j,k).ne.0.and.ibf(i,j,k-1).ne.0) then
       h(i,j,k)=0.0d0
       go to 3001
      end if
c
c++++++++++++++++++++++++++++++++++++
      if(ibf(i,j,k+1).ne.0) ww(i,j,k+2)=ca1*ww(i,j,k)
      if(ibf(i,j,k-2).ne.0) ww(i,j,k-2)=ca1*ww(i,j,k)
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+1,k).ne.0.and.ibf(i-1,j+1,k).ne.0)
     &    uu(i,j+1,k)=ca1*uu(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &    uu(i,j-1,k)=ca1*uu(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i-1,j,k+1).ne.0)
     &    uu(i,j,k+1)=ca1*uu(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &    uu(i,j,k-1)=ca1*uu(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i,j+2,k).ne.0.and.ibf(i-1,j+2,k).ne.0)
     &    uu(i,j+2,k)=ca1*uu(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i-1,j-2,k).ne.0)
     &    uu(i,j-2,k)=ca1*uu(i,j-1,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i-1,j,k+2).ne.0)
     &    uu(i,j,k+2)=ca1*uu(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i-1,j,k-2).ne.0)
     &    uu(i,j,k-2)=ca1*uu(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j-1,k).ne.0)
     &    vv(i+1,j,k)=ca1*vv(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j-1,k).ne.0)
     &    vv(i-1,j,k)=ca1*vv(i,j,k)
      if(ibf(i,j,k+1).ne.0.and.ibf(i,j-1,k+1).ne.0)
     &    vv(i,j,k+1)=ca1*vv(i,j,k)
      if(ibf(i,j,k-1).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &    vv(i,j,k-1)=ca1*vv(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j-1,k).ne.0)
     &    vv(i+2,j,k)=ca1*vv(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j-1,k).ne.0)
     &    vv(i-2,j,k)=ca1*vv(i-1,j,k)
      if(ibf(i,j,k+2).ne.0.and.ibf(i,j-1,k+2).ne.0)
     &    vv(i,j,k+2)=ca1*vv(i,j,k+1)
      if(ibf(i,j,k-2).ne.0.and.ibf(i,j-1,k-2).ne.0)
     &    vv(i,j,k-2)=ca1*vv(i,j,k-1)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j,k-1).ne.0)
     &    ww(i+1,j,k)=ca1*ww(i,j,k)
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j,k-1).ne.0)
     &    ww(i-1,j,k)=ca1*ww(i,j,k)
      if(ibf(i,j+1,k).ne.0.and.ibf(i,j+1,k-1).ne.0)
     &    ww(i,j+1,k)=ca1*ww(i,j,k)
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j-1,k-1).ne.0)
     &    ww(i,j-1,k)=ca1*ww(i,j,k)
      end if
c
      if(ibf(i,j,k).eq.0) then
      if(ibf(i+2,j,k).ne.0.and.ibf(i+2,j,k-1).ne.0)
     &    ww(i+2,j,k)=ca1*ww(i+1,j,k)
      if(ibf(i-2,j,k).ne.0.and.ibf(i-2,j,k-1).ne.0)
     &    ww(i-2,j,k)=ca1*ww(i-1,j,k)
      if(ibf(i,j+2,k).ne.0.and.ibf(i,j+2,k-1).ne.0)
     &    ww(i,j+2,k)=ca1*ww(i,j+1,k)
      if(ibf(i,j-2,k).ne.0.and.ibf(i,j-2,k-1).ne.0)
     &    ww(i,j-2,k)=ca1*ww(i,j-1,k)
      end if
c++++++++++++++++++++++++++++++++++++
c
      if(ibf(i+1,j,k).ne.0) kcpx=1
      if(ibf(i-1,j,k).ne.0) kcpx=2 
      if(ibf(i,j+1,k).ne.0) kcpy=1
      if(ibf(i,j-1,k).ne.0) kcpy=2
c
c------------------------------------------ for around b-body
c
      wijk=ww(i,j,k)
      wipjk=ww(i+1,j,k)
      wimjk=ww(i-1,j,k)
      wijpk=ww(i,j+1,k)
      wijmk=ww(i,j-1,k)
      wijkp=ww(i,j,k+1)
      wijkm=ww(i,j,k-1)
c
c --------------------- Convective terms in H -----------------------
      uw=0.5d0*((dz(k)*(uu(i,j,k)+uu(i+1,j,k))
     &        +dz(k+1)*(uu(i,j,k-1)+uu(i+1,j,k-1)))/(dz(k)+dz(k+1)))      
      vw=0.5d0*((dz(k)*(vv(i,j,k)+vv(i,j+1,k))
     &        +dz(k+1)*(vv(i,j,k-1)+vv(i,j+1,k-1)))/(dz(k)+dz(k+1)))
c
c-------------------------------------------------  Udw/dx
      dwdxp=(wipjk-wijk)/dxh(i+1)
      dwdxm=(wijk-wimjk)/dxh(i)
c
c ---------------------- udwdx, U>=0
      if(uw.ge.0.0) then
       if((i.eq.is).or.(kcpx.eq.2)) then
        udwdx=uw*dwdxm
       else
        dx21=0.5d0*(x(i)-x(i-2))
        dx24=0.5d0*(x(i-1)+x(i)-x(i+1)-x(i+2))
        dx31=0.5d0*(x(i)+x(i+1)-x(i-2)-x(i-1))
        dx32=0.5d0*(x(i+1)-x(i-1))
        dx34=0.5d0*(x(i)-x(i+2))
        dx41=0.5d0*(x(i+1)+x(i+2)-x(i-2)-x(i-1))
         xa=-dx32*dx34/(dx21*dx31*dx41)
         xb=-dx31*dx34/(dx21*dx32*dx24)
         xc=(1.0d0/dx31+1.0d0/dx32+1.0d0/dx34)
         xd=dx31*dx32/(dx41*dx24*dx34)
          udwdx=uw*(xa*ww(i-2,j,k)+xb*wimjk+xc*wijk+xd*wipjk)
       endif
      endif
c
c ---------------------- udwdx, U<0
      if(uw.lt.0.0) then
       if((i.eq.ie-1).or.(kcpx.eq.1)) then
        udwdx=uw*dwdxp
       else
        dx24=0.5d0*(x(i-1)+x(i)-x(i+1)-x(i+2))
        dx25=0.5d0*(x(i-1)+x(i)-x(i+2)-x(i+3))
        dx32=0.5d0*(x(i+1)-x(i-1))
        dx34=0.5d0*(x(i)-x(i+2))
        dx35=0.5d0*(x(i)+x(i+1)-x(i+2)-x(i+3))
        dx45=0.5d0*(x(i+1)-x(i+3))
         xe=-dx34*dx35/(dx32*dx24*dx25)
         xf=(1.0d0/dx32+1.0d0/dx34+1.0d0/dx35)
         xg=dx32*dx35/(dx24*dx34*dx45)
         xh=-dx32*dx34/(dx25*dx35*dx45)
         udwdx=uw*(xe*wimjk+xf*wijk+xg*wipjk+xh*ww(i+2,j,k))
       endif
      endif
c
c-------------------------------------------------  Vdw/dy
      dwdyp=(wijpk-wijk)/dyh(j+1)
      dwdym=(wijk-wijmk)/dyh(j)
c
c ---------------------- vdwdy, V>=0
      if(vw.ge.0.0) then
       if((iyper.ne.0.and.j.eq.js).or.(kcpy.eq.2)) then
        vdwdy=vw*dwdym
       else 
        dy21=0.5d0*(y(j)-y(j-2))
        dy24=0.5d0*(y(j-1)+y(j)-y(j+1)-y(j+2))
        dy31=0.5d0*(y(j)+y(j+1)-y(j-2)-y(j-1))
        dy32=0.5d0*(y(j+1)-y(j-1))
        dy34=0.5d0*(y(j)-y(j+2))
        dy41=0.5d0*(y(j+1)+y(j+2)-y(j-2)-y(j-1))
         ya=-dy32*dy34/(dy21*dy31*dy41)
         yb=-dy31*dy34/(dy21*dy32*dy24)
         yc=(1.0d0/dy31+1.0d0/dy32+1.0d0/dy34)
         yd=dy31*dy32/(dy41*dy24*dy34)
          vdwdy=vw*(ya*ww(i,j-2,k)+yb*wijmk+yc*wijk+yd*wijpk)
        endif
      endif
c
c ---------------------- vdwdy, V>=0
      if(vw.lt.0.0) then
       if((iyper.ne.0.and.j.eq.je-1).or.(kcpy.eq.1)) then
        vdwdy=vw*dwdyp
       else
        dy24=0.5d0*(y(j-1)+y(j)-y(j+1)-y(j+2))
        dy25=0.5d0*(y(j-1)+y(j)-y(j+2)-y(j+3))
        dy32=0.5d0*(y(j+1)-y(j-1))
        dy34=0.5d0*(y(j)-y(j+2))
        dy35=0.5d0*(y(j)+y(j+1)-y(j+2)-y(j+3))
        dy45=0.5d0*(y(j+1)-y(j+3))
         ye=-dy34*dy35/(dy32*dy24*dy25)
         yf=(1.0d0/dy32+1.0d0/dy34+1.0d0/dy35)
         yg=dy32*dy35/(dy24*dy34*dy45)
         yh=-dy32*dy34/(dy25*dy35*dy45)
          vdwdy=vw*(ye*wijmk+yf*wijk+yg*wijpk+yh*ww(i,j+2,k))
       endif
      endif
c
c-------------------------------------------------  Wdw/dz
      dwdzp=(wijkp-wijk)/dz(k+1)
      dwdzm=(wijk-wijkm)/dz(k)
c
c ---------------------- wdwdz, W>=0
      if(wijk.ge.0.) then
       if(izper.ne.0.and.k.eq.ks) then 
        wdwdz=wijk*dwdzm
       else
        dz21=z(k-1)-z(k-2)
        dz24=z(k-1)-z(k+1)
        dz31=z(k)-z(k-2)
        dz32=z(k)-z(k-1)
        dz34=z(k)-z(k+1)
        dz41=z(k+1)-z(k-2)
         za=-dz32*dz34/(dz21*dz31*dz41)
         zb=-dz31*dz34/(dz21*dz32*dz24)
         zc=(1.0d0/dz31+1.0d0/dz32+1.0d0/dz34)
         zd=dz31*dz32/(dz41*dz24*dz34)
          wdwdz=wijk*(za*ww(i,j,k-2)+zb*wijkm+zc*wijk+zd*wijkp)
       endif
      endif
c
c ---------------------- wdwdz, W<0
      if(wijk.lt.0.) then
       if(izper.ne.0.and.k.eq.ke) then
        wdwdz=wijk*dwdzp
       else
        dz24=z(k-1)-z(k+1)
        dz25=z(k-1)-z(k+2)
        dz32=z(k)-z(k-1)
        dz34=z(k)-z(k+1)
        dz35=z(k)-z(k+2)
        dz45=z(k+1)-z(k+2)
         ze=-dz34*dz35/(dz32*dz24*dz25)
         zf=(1.0d0/dz32+1.0d0/dz34+1.0d0/dz35)
         zg=dz32*dz35/(dz24*dz34*dz45)
         zh=-dz32*dz34/(dz25*dz35*dz45)
          wdwdz=wijk*(ze*wijkm+zf*wijk+zg*wijkp+zh*ww(i,j,k+2))
        endif
      endif
c
c
c --------------    Viscous and turbulence terms in H
      dwdxx=(dwdxp-dwdxm)/((dxh(i)+dxh(i+1))/2.0d0)
      dwdyy=(dwdyp-dwdym)/((dyh(j)+dyh(j+1))/2.0d0)
      dwdzz=(dwdzp-dwdzm)/dzh(k)
c------------------ Sij for W 
      dwdx=0.5d0*(dwdxp+dwdxm)
      dwdy=0.5d0*(dwdyp+dwdym)
      dwdz=0.5d0*(dwdzp+dwdzm)
c
      uijkp=(uu(i,j,k  )+uu(i+1,j,k  ))/2.0d0
      uijkm=(uu(i,j,k-1)+uu(i+1,j,k-1))/2.0d0
      dudz=(uijkp-uijkm)/dzh(k)
c
      vijkp=(vv(i,j,k  )+vv(i,j+1,k  ))/2.0d0
      vijkm=(vv(i,j,k-1)+vv(i,j+1,k-1))/2.0d0
      dvdz=(vijkp-vijkm)/dzh(k)
c
c
      if(ibf(i,j,k).ne.0.and.ibf(i,j,k-1).eq.0)
     &   tnu(i,j,k)=tnu(i,j,k-1)
      if(ibf(i,j,k).eq.0.and.ibf(i,j,k-1).ne.0)
     &   tnu(i,j,k-1)=tnu(i,j,k)
c
      tnuijk=(tnu(i,j,k-1)*dz(k+1)+tnu(i,j,k)*dz(k))/(dzh(k)*2.0d0)
c
      if(ibf(i,j,k).eq.0.and.ibf(i+1,j,k).ne.0) tnu(i+1,j,k)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i-1,j,k).ne.0) tnu(i-1,j,k)=tnu(i,j,k)
      if(ibf(i,j,k-1).eq.0.and.ibf(i+1,j,k-1).ne.0)
     &   tnu(i+1,j,k-1)=tnu(i,j,k-1)
      if(ibf(i,j,k-1).eq.0.and.ibf(i-1,j,k-1).ne.0)
     &   tnu(i-1,j,k-1)=tnu(i,j,k-1)
c
      tnipjk=(tnu(i+1,j,k-1)*dz(k+1)+tnu(i+1,j,k)*dz(k))
     &      /(dzh(k)*2.0d0)
      tnimjk=(tnu(i-1,j,k-1)*dz(k+1)+tnu(i-1,j,k)*dz(k))
     &      /(dzh(k)*2.0d0)
      if(ibf(i+1,j,k-1).ne.0.and.ibf(i+1,j,k).eq.0)
     &  tnipjk=tnu(i+1,j,k)
      if(ibf(i+1,j,k-1).eq.0.and.ibf(i+1,j,k).ne.0)
     &  tnipjk=tnu(i+1,j,k-1)
      if(ibf(i-1,j,k-1).ne.0.and.ibf(i-1,j,k).eq.0)
     &  tnimjk=tnu(i-1,j,k)
      if(ibf(i-1,j,k-1).eq.0.and.ibf(i-1,j,k).ne.0)
     &  tnimjk=tnu(i-1,j,k-1)
c
      if(ibf(i+1,j,k).ne.0.and.ibf(i+1,j,k-1).ne.0) then
       if(ibf(i,j,k-1).eq.0.and.ibf(i,j,k).ne.0) tnipjk=tnu(i,j,k-1)
       if(ibf(i,j,k-1).ne.0.and.ibf(i,j,k).eq.0) tnipjk=tnu(i,j,k)
      end if
      if(ibf(i-1,j,k).ne.0.and.ibf(i-1,j,k-1).ne.0) then
       if(ibf(i,j,k-1).eq.0.and.ibf(i,j,k).ne.0) tnimjk=tnu(i,j,k-1)
       if(ibf(i,j,k-1).ne.0.and.ibf(i,j,k).eq.0) tnimjk=tnu(i,j,k)
      end if
c
      dndx=(tnipjk-tnimjk)/(dx(i+1)+(dx(i+2)+dx(i))/2.0d0)
c
      if(ibf(i,j,k).eq.0.and.ibf(i,j+1,k).ne.0) tnu(i,j+1,k)=tnu(i,j,k)
      if(ibf(i,j,k).eq.0.and.ibf(i,j-1,k).ne.0) tnu(i,j-1,k)=tnu(i,j,k)
      if(ibf(i,j,k-1).eq.0.and.ibf(i,j+1,k-1).ne.0)
     &   tnu(i,j+1,k-1)=tnu(i,j,k-1)
      if(ibf(i,j,k-1).eq.0.and.ibf(i,j-1,k-1).ne.0)
     &   tnu(i,j-1,k-1)=tnu(i,j,k-1)
c
      tnijpk=(tnu(i,j+1,k-1)*dz(k+1)+tnu(i,j+1,k)*dz(k))
     &      /(dzh(k)*2.0d0)
      tnijmk=(tnu(i,j-1,k-1)*dz(k+1)+tnu(i,j-1,k)*dz(k))
     &      /(dzh(k)*2.0d0)
      if(ibf(i,j+1,k-1).ne.0.and.ibf(i,j+1,k).eq.0)
     &  tnijpk=tnu(i,j+1,k)
      if(ibf(i,j+1,k-1).eq.0.and.ibf(i,j+1,k).ne.0)
     &  tnijpk=tnu(i,j+1,k-1)
      if(ibf(i,j-1,k-1).ne.0.and.ibf(i,j-1,k).eq.0)
     &  tnijmk=tnu(i,j-1,k)
      if(ibf(i,j-1,k-1).eq.0.and.ibf(i,j-1,k).ne.0)
     &  tnijmk=tnu(i,j-1,k-1)
c
      if(ibf(i,j+1,k).ne.0.and.ibf(i,j+1,k-1).ne.0) then
       if(ibf(i,j,k-1).eq.0.and.ibf(i,j,k).ne.0) tnijpk=tnu(i,j,k-1)
       if(ibf(i,j,k-1).ne.0.and.ibf(i,j,k).eq.0) tnijpk=tnu(i,j,k)
      end if
      if(ibf(i,j-1,k).ne.0.and.ibf(i,j-1,k-1).ne.0) then
       if(ibf(i,j,k-1).eq.0.and.ibf(i,j,k).ne.0) tnijmk=tnu(i,j,k-1)
       if(ibf(i,j,k-1).ne.0.and.ibf(i,j,k).eq.0) tnijmk=tnu(i,j,k)
      end if
      dndy=(tnijpk-tnijmk)/(dy(j+1)+(dy(j+2)+dy(j))/2.0d0)
c
      dndz=(tnu(i,j,k)-tnu(i,j,k-1))/dzh(k)
c
      h(i,j,k)=-(udwdx+vdwdy+wdwdz)
     &        +(1.0d0/re+tnuijk)*(dwdxx+dwdyy+dwdzz)
     &        +dndx*(dwdx+dudz)+dndy*(dwdy+dvdz)+2.0d0*dndz*dwdz
c
c
c
 3001 dmy=0.0d0
c
  300 continue
c             
c
c
      return
      end
c----------------------------------------------------------- sub.euler2
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c234567----------------------------------------------------- sub.euler2
      subroutine euler2(ifile)
      include 'common_vib8.fil'
      character*14 ifile
c
      resu=0.0d0
      resv=0.0d0
      resw=0.0d0
c
c-----------------------------------------    RESU   
      do 100 k=ks,ke-1
      do 100 j=js,je-1
      do 100 i=is+1,ie-1
c********************************************* mod. for cal. vibration
       ures=dt*(0.5d0*(3.0d0*(f(i,j,k) -xacc -facc1(i,j,k))
     &                      -(fp(i,j,k)-xacc0-facc0(i,j,k)))
     &     -(pp(i,j,k)-pp(i-1,j,k))/dxh(i))
       uu(i,j,k)=uu(i,j,k)+ures
         resu=resu+ures*ures
c
  100 continue
c
c-----------------------------------------    RESV   
      do 200 k=ks,ke-1
      do 200 j=js,je
      do 200 i=is+1,ie-1
c********************************************* mod. for cal. vibration
       vres=dt*(0.5d0*(3.0d0*(g(i,j,k) -yacc -gacc1(i,j,k))
     &                      -(gp(i,j,k)-yacc0-gacc0(i,j,k)))
     &     -(pp(i,j,k)-pp(i,j-1,k))/dyh(j))
       vv(i,j,k)=vv(i,j,k)+vres
          resv=resv+vres*vres
  200 continue
c
c      if(iyper.ne.0) then
c       do 210 k=ks,ke-1
c       do 210 i=is,ie-1
c        vv(i,js,k)=0.0d0
c        vv(i,je,k)=0.0d0
c  210  continue
c      end if 
c
c-----------------------------------------    RESW   
      do 300 k=ks,ke
      do 300 j=js,je-1
      do 300 i=is+1,ie-1
       wres=dt*(0.5d0*(3.0d0*h(i,j,k)-hp(i,j,k))
     &     -(pp(i,j,k)-pp(i,j,k-1))/dzh(k))
       ww(i,j,k)=ww(i,j,k)+wres
          resw=resw+wres*wres
  300 continue
c
      if(izper.ne.0) then
       do 310 i=is,ie-1
       do 310 j=js,je-1
        ww(i,j,ks)=0.0d0
        ww(i,j,ke)=0.0d0
  310 continue
      end if 
c
c----------------------------------------- for around b-body
      do 10 k=ks,ke-1
      do 10 j=js,je-1
      do 10 i=is,ie-1
       if(ibf(i,j,k).ne.0) then
        uu(i,j,k)=0.0d0
        vv(i,j,k)=0.0d0
        ww(i,j,k)=0.0d0
        uu(i+1,j,k)=0.0d0
        vv(i,j+1,k)=0.0d0
        ww(i,j,k+1)=0.0d0
       end if
   10 continue
c
c
      resu=dsqrt(resu/dble((ie-1-is)*(je-js)*(ke-ks)))
      resv=dsqrt(resv/dble((ie-1-is)*(je-js+1)*(ke-ks)))
      resw=dsqrt(resw/dble((ie-1-is)*(je-js)*(ke-ks+1)))
c
      call rdvelo(ifile)
      iswch(1)=1
      iswch(2)=1
      iswch(3)=1
      call uvwbdy
c
      return
      end
c---------------------------------------------------------- sub.uvwrvs
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.uvwrvs
      subroutine uvwrvs
      include 'common_vib8.fil'
c
c---------------     u(i,j,k) -> uu(i,j,k)
      do 10 i=is,ie
       do 10 j=js,je-1
        do 10 k=ks,ke-1
         uu(i,j,k)=0.25d0*(u(i,j,k)+u(i,j+1,k)
     &            +u(i,j+1,k+1)+u(i,j,k+1))
   10 continue
c
c---------------     v(i,j,k) -> vv(i,j,k) 
      do 20 i=is,ie-1
       do 20 j=js,je
        do 20 k=ks,ke-1
         vv(i,j,k)=0.25d0*(v(i,j,k)+v(i+1,j,k)
     &            +v(i+1,j,k+1)+v(i,j,k+1))
   20 continue
c
c---------------     w(i,j,k) -> ww(i,j,k) 
      do 30 i=is,ie-1
       do 30 j=js,je-1
        do 30 k=ks,ke
         ww(i,j,k)=0.25d0*(w(i,j,k)+w(i+1,j,k)
     &            +w(i+1,j+1,k)+w(i,j+1,k))
   30 continue
c
c
      return
      end
c---------------------------------------------------------- sub.prsrvs
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.prsrvs
      subroutine prsrvs
      include 'common_vib8.fil'
c
c---------------     p(i,j,k) -> pp(i,j,k)   
      do 10 i=is,ie-1
      do 10 j=js,je-1
      do 10 k=ks,ke-1
       pp(i,j,k)=(p(i,j,k)+p(i+1,j,k)
     &          +p(i+1,j,k+1)+p(i,j,k+1)+p(i,j+1,k)
     &          +p(i+1,j+1,k)+p(i+1,j+1,k+1)+p(i,j+1,k+1))*0.125d0
   10 continue
c
      call prsbdy
c
c
      return
      end
c----------------------------------------------------------- sub.prefgh
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c234567----------------------------------------------------- sub.prefgh
      subroutine prefgh
      include 'common_vib8.fil'
c
      do 10 k=ks,ke
      do 10 j=js,je
      do 10 i=is,ie
       fp(i,j,k)=f(i,j,k)
       gp(i,j,k)=g(i,j,k)
       hp(i,j,k)=h(i,j,k)
   10 continue
c
      do 20 k=ks,ke-1
      do 20 j=js,je-1
       uo(j,k,2)=uu(ie,j,k)
       uo(j,k,1)=uu(ie-1,j,k)
c**********************************************************      11/19/98
       if(uo(j,k,2).lt.0.0) uo(j,k,2)=0.0
       if(uo(j,k,1).lt.0.0) uo(j,k,1)=0.0
c************************************************************************
   20 continue
      do 30 k=ks,ke-1
      do 30 j=js,je
       vo(j,k,2)=vv(ie,j,k)
       vo(j,k,1)=vv(ie-1,j,k)
   30 continue
      do 40 k=ks,ke
      do 40 j=js,je-1
       wo(j,k,2)=ww(ie,j,k)
       wo(j,k,1)=ww(ie-1,j,k)
   40 continue
c
c
      return
      end
c---------------------------------------------------------- sub.chguvw
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.chguvw
      subroutine chguvw
      include 'common_vib8.fil'
c
      do 10 k=ks,ke
      do 10 j=js,je
      do 10 i=is,ie
       ukp=(dy(j)*uu(i,j,k)+dy(j+1)*uu(i,j-1,k))/(2.0d0*dyh(j))
       ukm=(dy(j)*uu(i,j,k-1)+dy(j+1)*uu(i,j-1,k-1))/(2.0d0*dyh(j))
        u(i,j,k)=(dz(k)*ukp+dz(k+1)*ukm)/(2.0d0*dzh(k))
c
       vip=(dz(k)*vv(i,j,k)+dz(k+1)*vv(i,j,k-1))/(2.0d0*dzh(k))
       vim=(dz(k)*vv(i-1,j,k)+dz(k+1)*vv(i-1,j,k-1))/(2.0d0*dzh(k))
        v(i,j,k)=(dx(i)*vip+dx(i+1)*vim)/(2.0d0*dxh(i))
c
       wjp=(dx(i)*ww(i,j,k)+dx(i+1)*ww(i-1,j,k))/(2.0d0*dxh(i))
       wjm=(dx(i)*ww(i,j-1,k)+dx(i+1)*ww(i-1,j-1,k))/(2.0d0*dxh(i))
        w(i,j,k)=(dy(j)*wjp+dy(j+1)*wjm)/(2.0d0*dyh(j))
   10 continue
c
c
      return
      end
c----------------------------------------------------------- sub.uvwbd2
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c                                                                     |
c234567----------------------------------------------------- sub.uvwbd2
      subroutine uvwbdy
      include 'common_vib8.fil'
c
      uinit=1.0d0
c
      do 11 k=ks,ke-1
      do 11 j=js,je-1
      do 11 i=is,ie-1
       if(ibf(i,j,k).ne.0.and.ibf(i-1,j,k).eq.0) then
        vv(i,j,k)=ca1*vv(i-1,j,k)
        ww(i,j,k)=ca1*ww(i-1,j,k)
        vv(i,j+1,k)=ca1*vv(i-1,j+1,k)
        ww(i,j,k+1)=ca1*ww(i-1,j,k+1)
       end if
       if(ibf(i,j,k).eq.0.and.ibf(i-1,j,k).ne.0) then
        vv(i-1,j,k)=ca1*vv(i,j,k)
        ww(i-1,j,k)=ca1*ww(i,j,k)
        vv(i-1,j+1,k)=ca1*vv(i,j+1,k)
        ww(i-1,j,k+1)=ca1*ww(i,j,k+1)
       end if
c
       if(ibf(i,j,k).ne.0.and.ibf(i,j-1,k).eq.0) then
        uu(i,j,k)=ca1*uu(i,j-1,k)
        ww(i,j,k)=ca1*ww(i,j-1,k)
        uu(i+1,j,k)=ca1*uu(i+1,j-1,k)
        ww(i,j,k+1)=ca1*ww(i,j-1,k+1)
       end if
       if(ibf(i,j,k).eq.0.and.ibf(i,j-1,k).ne.0) then
        uu(i,j-1,k)=ca1*uu(i,j,k)
        ww(i,j-1,k)=ca1*ww(i,j,k)
        uu(i+1,j-1,k)=ca1*uu(i+1,j,k)
        ww(i,j-1,k+1)=ca1*ww(i,j,k+1)
       end if
c
       if(ibf(i,j,k).ne.0.and.ibf(i,j,k-1).eq.0) then
        uu(i,j,k)=ca1*uu(i,j,k-1)
        vv(i,j,k)=ca1*vv(i,j,k-1)
        uu(i+1,j,k)=ca1*uu(i+1,j,k-1)
        vv(i,j+1,k)=ca1*vv(i,j+1,k-1)
       end if
       if(ibf(i,j,k).eq.0.and.ibf(i,j,k-1).ne.0) then
        uu(i,j,k-1)=ca1*uu(i,j,k)
        vv(i,j,k-1)=ca1*vv(i,j,k)
        uu(i+1,j,k-1)=ca1*uu(i+1,j,k)
        vv(i,j+1,k-1)=ca1*vv(i,j+1,k)
       end if

   11 continue
c
      do 12 k=ks,ke-1
      do 12 j=js,je-1
      do 12 i=is,ie-1
       if(ibf(i-1,j,k).eq.0.and.ibf(i,j,k).ne.0) uu(i,j,k)=0.0d0
       if(ibf(i,j-1,k).eq.0.and.ibf(i,j,k).ne.0) vv(i,j,k)=0.0d0
       if(ibf(i,j,k-1).eq.0.and.ibf(i,j,k).ne.0) ww(i,j,k)=0.0d0
       if(ibf(i-1,j,k).ne.0.and.ibf(i,j,k).eq.0) uu(i,j,k)=0.0d0
       if(ibf(i,j-1,k).ne.0.and.ibf(i,j,k).eq.0) vv(i,j,k)=0.0d0
       if(ibf(i,j,k-1).ne.0.and.ibf(i,j,k).eq.0) ww(i,j,k)=0.0d0
   12 continue
c
c
c------------------------------------------ inflow , outflow
      rdx=1.0d0/dx(ie)
c
c      do 40 k=ks,ke-1
c      do 40 j=js,je-1
c        uu(is-1,j,k)=uu(is,j,k)
c   40 continue
c
      qall=0.0d0
      do 410 k=ks,ke-1
      do 410 j=js,je-1
       qall=qall+uu(is,j,k)*dy(j+1)*dz(k+1)
       qall=qall-uu(ie,j,k)*dy(j+1)*dz(k+1)
  410 continue
      do 420 i=is,ie-1
      do 420 k=ks,ke-1
       qall=qall+vv(i,js,k)*dz(k+1)*dx(i+1)
       qall=qall-vv(i,je,k)*dz(k+1)*dx(i+1)
  420 continue
      do 430 i=is,ie-1
      do 430 j=js,je-1
       qall=qall+ww(i,j,ks)*dx(i+1)*dy(j+1)
       qall=qall-ww(i,j,ke)*dx(i+1)*dy(j+1)
  430 continue
      sa=(z(ke)-z(ks))*(y(je)-y(js))
     &  +2.0d0*(x(ie)-x(is))*(y(je)-y(js))
      duc=qall/sa
C 8/30/13 set duc=0 for fixed body case
      duc=0.d0
      print *,'Qall ,cancel U=',qall,duc
c
          
      if(iswch(1).eq.1) then
      do 41 k=ks,ke-1
      do 41 j=js,je-1
c        uu(ie,j,k)=uo(j,k,2)-ui(j,k)*dt*rdx*(uo(j,k,2)-uo(j,k,1))
       yyy=(y(j)+y(j+1))*0.5d0
       fth=datan2(yyy,x(ie))
       fr=dsqrt(yyy**2+x(ie)**2)
       uinit2=uinit-tvel*fr*(-dsin(fth+tdis))
c        uu(ie,j,k)=uo(j,k,2)-uinit2*dt*rdx*(uo(j,k,2)-uo(j,k,1))
c
      uu(ie,j,k)=uu(ie-1,j,k)+duc
   41 continue
      end if
c
c
c      do 50 k=ks,ke-1
c      do 50 j=js,je
c        vv(is-1,j,k)=vv(is,j,k)
c   50 continue
      if(iswch(2).eq.1) then
      do 51 k=ks,ke-1
      do 51 j=js,je
c        vv(ie,j,k)=vo(j,k,2)-ui(j,k)*dt*rdx*(vo(j,k,2)-vo(j,k,1))
       xxx=(x(ie)+x(ie+1))*0.5d0
       fth=datan2(y(j),xxx)
       fr=dsqrt(y(j)**2+xxx**2)
       uinit2=uinit-tvel*fr*(-dsin(fth+tdis))
c        vv(ie,j,k)=vo(j,k,2)-uinit2*dt*rdx*(vo(j,k,2)-vo(j,k,1))
c
      vv(ie,j,k)=vv(ie-1,j,k)
   51 continue
      end if
c
c      do 60 k=ks,ke
c      do 60 j=js,je-1
c        ww(is-1,j,k)=ww(is,j,k)
c   60 continue
      if(iswch(3).eq.1) then
      do 61 k=ks,ke
      do 61 j=js,je-1
c        ww(ie,j,k)=wo(j,k,2)-ui(j,k)*dt*rdx*(wo(j,k,2)-wo(j,k,1))
       xxx=(x(ie)+x(ie+1))*0.5d0
       fth=datan2(y(j),xxx)
       fr=dsqrt(y(j)**2+xxx**2)
       uinit2=uinit-tvel*fr*(-dsin(fth+tdis))
c        ww(ie,j,k)=wo(j,k,2)-uinit2*dt*rdx*(wo(j,k,2)-wo(j,k,1))
c
      ww(ie,j,k)=ww(ie-1,j,k)
   61 continue
      end if
c
c
c------------------------------------------ z-x boundary (y dir.)
c----------------------------------+  
      if(iyper.eq.0) then
       do 70 k=ks,ke-1
       do 70 i=is,ie
         uu(i,js-1,k)=uu(i,je-1,k)
         uu(i,js-2,k)=uu(i,je-2,k)
         uu(i,je,k)=uu(i,js,k)
         uu(i,je+1,k)=uu(i,js+1,k)
   70  continue
       do 80 k=ks,ke-1
       do 80 i=is,ie-1
         vv(i,js-1,k)=vv(i,je-1,k)
         vv(i,js-2,k)=vv(i,je-2,k)
         vv(i,je+1,k)=vv(i,js+1,k)
         vv(i,je+2,k)=vv(i,js+2,k)
C 8/31/13  periodic BC for vv(i,je,k)                         sub.uvwbdy
         vv(i,je,k)=vv(i,js,k)
   80  continue
       do 90 k=ks,ke
       do 90 i=is,ie-1
         ww(i,js-1,k)=ww(i,je-1,k)
         ww(i,js-2,k)=ww(i,je-2,k)
         ww(i,je,k)=ww(i,js,k)
         ww(i,je+1,k)=ww(i,js+1,k)
   90  continue
c
c-----------+  
      else
c-----------+  
c
      do 100 k=ks,ke-1
      do 100 i=is+1,ie
        uu(i,js-1,k)=uu(i,js,k)
        uu(i,je,k)=uu(i,je-1,k)
  100 continue
      do 110 k=ks,ke-1
      do 110 i=is+1,ie-1
cc        vv(i,js-1,k)=-vv(i,js+1,k)
cc        vv(i,je+1,k)=-vv(i,je-1,k)
cc        vv(i,js,k)=0.0d0
cc        vv(i,je,k)=0.0d0
        vv(i,js,k)=vv(i,js,k)-duc
        vv(i,je,k)=vv(i,je,k)+duc
        vv(i,js-1,k)=vv(i,js,k)
        vv(i,je+1,k)=vv(i,je,k)
C 8/31/13 slip cond. for fixed body                      sub.uvwbdy
C 9/ 2/13 slip cond. for transating body                 sub.uvwbdy
C        vv(i,js,k)=0.0d0
C        vv(i,je,k)=0.0d0
        vv(i,js,k)=-yvel
        vv(i,je,k)=-yvel
  110 continue
      do 120 k=ks,ke
      do 120 i=is+1,ie-1
        ww(i,js-1,k)=ww(i,js,k)
        ww(i,je,k)=ww(i,je-1,k)
  120 continue
c
       end if
c----------------------------------+  
c
c
c------------------------------------------ x-y boundary (z dir.)
c----------------------------------+  
      if(izper.eq.0) then
       do 130 j=js,je-1
       do 130 i=is,ie
         uu(i,j,ks-1)=uu(i,j,ke-1)
         uu(i,j,ks-2)=uu(i,j,ke-2)
         uu(i,j,ke)=uu(i,j,ks)
         uu(i,j,ke+1)=uu(i,j,ks+1)
  130  continue
       do 140 j=js,je
       do 140 i=is,ie-1
         vv(i,j,ks-1)=vv(i,j,ke-1)
         vv(i,j,ks-2)=vv(i,j,ke-2)
         vv(i,j,ke)=vv(i,j,ks)
         vv(i,j,ke+1)=vv(i,j,ks+1)
  140  continue
       do 150 j=js,je-1
       do 150 i=is,ie-1
         ww(i,j,ks-1)=ww(i,j,ke-1)
         ww(i,j,ks-2)=ww(i,j,ke-2)
         ww(i,j,ke+1)=ww(i,j,ks+1)
         ww(i,j,ke+2)=ww(i,j,ks+2)
c--------------------------------------------add. 120927
         ww(i,j,ke)=ww(i,j,ks)
  150  continue
c
c-----------+  
      else
c-----------+  
c
      do 160 j=js,je-1
      do 160 i=is,ie
        uu(i,j,ks-1)=uu(i,j,ks)
        uu(i,j,ke)=uu(i,j,ke-1)
  160 continue
      do 170 j=js,je
      do 170 i=is,ie-1
        vv(i,j,ks-1)=vv(i,j,ks)
        vv(i,j,ke)=vv(i,j,ke-1)
  170 continue
      do 180 j=js,je-1
      do 180 i=is,ie-1
        ww(i,j,ks-1)=-ww(i,j,ks+1)
        ww(i,j,ke+1)=-ww(i,j,ke-1)
        ww(i,j,ks)=0.0d0
        ww(i,j,ke)=0.0d0
  180 continue
c
      end if
c----------------------------------+  
c----------------------------------+  
      if(izper.eq.2) then
       do 190 j=js,je-1
       do 190 i=is,ie
         uu(i,j,ks-1)=ca2*uu(i,j,ks)
  190 continue
       do 200 j=js,je
       do 200 i=is,ie-1
         vv(i,j,ks-1)=ca2*vv(i,j,ks)
  200 continue
      end if
c----------------------------------+  
c
c
c----------------------------------------- corner y-z plane (x dir.)
c----------------------------------+ 1 top 
      if(iyper.eq.0) then
c----------------------------------+ 2 top 
       if(izper.eq.0) then
        do 210 i=is,ie
         uu(i,js-1,ks-1)=uu(i,je-1,ke-1)
         uu(i,je,ks-1)=uu(i,js,ke-1)
         uu(i,js-1,ke)=uu(i,je-1,ks)
         uu(i,je,ke)=uu(i,js,ks)
  210   continue
       else
        do 220 i=is,ie
         uu(i,js-1,ks-1)=uu(i,je-1,ks)
         uu(i,je,ks-1)=uu(i,js,ks)
         uu(i,js-1,ke)=uu(i,je-1,ke-1)
         uu(i,je,ke)=uu(i,js,ke-1)
  220   continue
       end if
c----------------------------------+ 2 btm
      if(izper.eq.2) then
       do 230 i=is,ie 
        uu(i,js-1,ks-1)=-uu(i,je-1,ks)
        uu(i,je,ks-1)=-uu(i,js,ks)
  230  continue
      end if
c
      else
c----------------------------------+ 2 top 
      if(izper.eq.0) then
       do 240 i=is,ie
         uu(i,js-1,ks-1)=uu(i,js,ke-1)
         uu(i,je,ks-1)=uu(i,je-1,ke-1)
         uu(i,js-1,ke)=uu(i,js,ks)
         uu(i,je,ke)=uu(i,je-1,ks)
  240  continue
      else
       do 250 i=is,ie
         uu(i,js-1,ks-1)=uu(i,js,ks)
         uu(i,je,ks-1)=uu(i,je-1,ks)
         uu(i,js-1,ke)=uu(i,js,ke-1)
         uu(i,je,ke)=uu(i,je-1,ke-1)
  250  continue
      end if
c----------------------------------+ 2 btm 
      if(izper.eq.2) then
       do 260 i=is,ie
        uu(i,js-1,ks-1)=ca2*uu(i,js,ks)
        uu(i,je,ks-1)=ca2*uu(i,je-1,ks)
  260  continue
      end if
c
      end if       
c----------------------------------+ 1 btm 
c
c
c----------------------------------------- corner z-x plane (y dir.)
      if(izper.eq.0) then
       do 270 j=js,je
        vv(is-1,j,ks-1)=vv(is,j,ke-1)
        vv(ie,j,ks-1)=vv(ie-1,j,ke-1)
        vv(is-1,j,ke)=vv(is,j,ks)
        vv(ie,j,ke)=vv(ie-1,j,ks)
  270  continue
      else
       do 280 j=js,je
        vv(is-1,j,ks-1)=vv(is,j,ks)
        vv(ie,j,ks-1)=vv(ie-1,j,ks)
        vv(is-1,j,ke)=vv(is,j,ke-1)
        vv(ie,j,ke)=vv(ie-1,j,ke-1)
  280  continue
       if(izper.eq.2) then
        do 290 j=js,je
         vv(is-1,j,ks-1)=ca2*vv(is,j,ks)
         vv(ie,j,ks-1)=ca2*vv(ie-1,j,ks)
  290   continue
       end if
      end if
c
c
c----------------------------------------- corner x-y plane (z dir.)
      if(iyper.eq.0) then
       do 300 k=ks,ke
        ww(is-1,js-1,k)=ww(is,je-1,k)
        ww(ie,js-1,k)=ww(ie-1,je-1,k)
        ww(is-1,je,k)=ww(is,js,k)
        ww(ie,je,k)=ww(ie-1,js,k)
  300  continue
      else
       do 310 k=ks,ke
        ww(is-1,js-1,k)=ww(is,js,k)
        ww(ie,js-1,k)=ww(ie-1,js,k)
        ww(is-1,je,k)=ww(is,je-1,k)
        ww(ie,je,k)=ww(ie-1,je-1,k)
  310  continue
      end if
c
c
      return
      end
c---------------------------------------------------------- sub.prsbdy
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.prsbdy
      subroutine prsbdy
      include 'common_vib8.fil'
c
      do 11 k=ks,ke-1
      do 11 j=js,je-1
      do 11 i=is,ie-1
       if(ibf(i,j,k).ne.0) pp(i,j,k)=0.0d0
   11 continue
c
c------------------------------------------ inflow 
      do 10 k=ks,ke-1
      do 10 j=js,je-1
        pp(is-1,j,k)=pp(is,j,k)
c        pp(ie,j,k)=3.0d0*pp(ie-1,j,k)-3.0d0*pp(ie-2,j,k)+pp(ie-3,j,k)
        pp(ie-1,j,k)=pp(ie-1,j,k)-duc*dx(ie)/dt
        pp(ie,j,k)=pp(ie-1,j,k)
   10 continue
c
c
c------------------------------------------ z-x boundary (y dir.)
      if(iyper.eq.0) then
       do 20 k=ks,ke-1
       do 20 i=is,ie-1
         pp(i,js-1,k)=pp(i,je-1,k)
         pp(i,je,k)=pp(i,js,k)
   20  continue
      else
      do 30 k=ks,ke-1
      do 30 i=is,ie-1
c        pp(i,js-1,k)=pp(i,js,k)
c        pp(i,je,k)=pp(i,je-1,k)
c       if(vv(i,js,k).ge.0.0d0) pp(i,js-1,k)=0.0d0
c       if(vv(i,je,k).le.0.0d0) pp(i,je  ,k)=0.0d0
        pp(i,js,k)=pp(i,js,k)-duc*dy(js+1)*rdt
        pp(i,je-1,k)=pp(i,je-1,k)+duc*dy(je)*rdt
        pp(i,js-1,k)=pp(i,js,k)
        pp(i,je,k)=pp(i,je-1,k)
   30 continue
      end if
c
c
c------------------------------------------ x-y boundary (z dir.)
      if(izper.eq.0) then
       do 40 j=js,je-1
       do 40 i=is,ie-1
         pp(i,j,ks-1)=pp(i,j,ke-1)
         pp(i,j,ke)=pp(i,j,ks)
   40  continue
      else
      do 50 j=js,je-1
      do 50 i=is,ie-1
        pp(i,j,ks-1)=pp(i,j,ks)
        pp(i,j,ke)=pp(i,j,ke-1)
   50 continue
      end if
c
c
c----------------------------------------- corner y-z plane (x dir.)
c----------------------------------+ 1 top 
      if(iyper.eq.0) then
c----------------------------------+ 2 top 
       if(izper.eq.0) then
        do 60 i=is,ie-1
         pp(i,js-1,ks-1)=pp(i,je-1,ke-1)
         pp(i,je,ks-1)=pp(i,js,ke-1)
         pp(i,js-1,ke)=pp(i,je-1,ks)
         pp(i,je,ke)=pp(i,js,ks)
   60   continue
       else
        do 70 i=is,ie-1
         pp(i,js-1,ks-1)=pp(i,je-1,ks)
         pp(i,je,ks-1)=pp(i,js,ks)
         pp(i,js-1,ke)=pp(i,je-1,ke-1)
         pp(i,je,ke)=pp(i,js,ke-1)
   70   continue
       end if
c----------------------------------+ 2 btm
      else
c----------------------------------+ 2 top 
      if(izper.eq.0) then
       do 80 i=is,ie-1
         pp(i,js-1,ks-1)=pp(i,js,ke-1)
         pp(i,je,ks-1)=pp(i,je-1,ke-1)
         pp(i,js-1,ke)=pp(i,js,ks)
         pp(i,je,ke)=pp(i,je-1,ks)
  80  continue
      else
       do 90 i=is,ie-1
         pp(i,js-1,ks-1)=pp(i,js,ks)
         pp(i,je,ks-1)=pp(i,je-1,ks)
         pp(i,js-1,ke)=pp(i,js,ke-1)
         pp(i,je,ke)=pp(i,je-1,ke-1)
  90  continue
      end if
c----------------------------------+ 2 btm 
      end if       
c----------------------------------+ 1 btm 
c
c
c----------------------------------------- corner z-x plane (y dir.)
      if(izper.eq.0) then
       do 100 j=js,je-1
        pp(is-1,j,ks-1)=pp(is,j,ke-1)
        pp(ie,j,ks-1)=pp(ie-1,j,ke-1)
        pp(is-1,j,ke)=pp(is,j,ks)
        pp(ie,j,ke)=pp(ie-1,j,ks)
  100  continue
      else
       do 110 j=js,je-1
        pp(is-1,j,ks-1)=pp(is,j,ks)
        pp(ie,j,ks-1)=pp(ie-1,j,ks)
        pp(is-1,j,ke)=pp(is,j,ke-1)
        pp(ie,j,ke)=pp(ie-1,j,ke-1)
  110  continue
      end if
c
c
c----------------------------------------- corner x-y plane (z dir.)
      if(iyper.eq.0) then
       do 120 k=ks,ke-1
        pp(is-1,js-1,k)=pp(is,je-1,k)
        pp(ie,js-1,k)=pp(ie-1,je-1,k)
        pp(is-1,je,k)=pp(is,js,k)
        pp(ie,je,k)=pp(ie-1,js,k)
  120  continue
      else
       do 130 k=ks,ke-1
        pp(is-1,js-1,k)=pp(is,js,k)
        pp(ie,js-1,k)=pp(ie-1,js,k)
        pp(is-1,je,k)=pp(is,je-1,k)
        pp(ie,je,k)=pp(ie-1,je-1,k)
  130  continue
      end if
c
c-------------------------------------------- corner
      pp(is-1,js-1,ks-1)=pp(is,js-1,ks-1)
      pp(is-1,je,ks-1)=pp(is,je,ks-1)
      pp(is-1,js-1,ke)=pp(is,js-1,ke)
      pp(is-1,je,ke)=pp(is,je,ke)
c
      pp(ie,js-1,ks-1)=pp(ie-1,js-1,ks-1)
      pp(ie,je,ks-1)=pp(ie-1,je,ks-1)
      pp(ie,js-1,ke)=pp(ie-1,js-1,ke)
      pp(ie,je,ke)=pp(ie-1,je,ke)
c
c
C 8/31/13 subtract ref pres                         sub.prsbdy
      pref=pp(is,je/2,ke/2)
      do 200 i=is-1,ie+1
      do 200 j=js-1,je+1
      do 200 k=ks-1,ke+1
        pp(i,j,k)=pp(i,j,k)-pref
  200 continue
      return
      end
c---------------------------------------------------------- sub.outdat
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.outdat
      subroutine outdat(ofile,kout,icount,ifcs)
      include 'common_vib8.fil'
c
      dimension p2(350,200,200)
      character*14 ofile(100)
c
      coef=2.0d0/3.0d0
      ck=0.094d0
      ck2=ck*ck
c
      do 10 k=ks,ke-1
      do 10 j=js,je-1
      do 10 i=is,ie
       del=(dx(i+1)*dy(j+1)*dz(k+1))**0.33333
       rks=tnu(i,j,k)*tnu(i,j,k)/(ck2*del*del)
       p2(i,j,k)=pp(i,j,k)-coef*rks
   10 continue
c
c
      iii=999
      print *,'********************* output DATA *********'
      print *,' file =',ofile(icount)
      print *,'********************* output DATA *********'
      open(12,file=ofile(icount),status='new',form='unformatted')
       write(12)is,ie,js,je,ks,ke
       write(12)re,cs,dt,t,nt,iyper,izper,iter
       write(12)(((uu(i,j,k),i=is,ie+1),j=js,je-1),k=ks,ke-1)
       write(12)(((vv(i,j,k),i=is,ie),j=js,je),k=ks,ke-1)
       write(12)(((ww(i,j,k),i=is,ie),j=js,je-1),k=ks,ke)
       write(12)(((p2(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
c
       write(12)(((f(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(12)(((g(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(12)(((h(i,j,k),i=is,ie),j=js,je),k=ks,ke)
c
       write(12)(((tnu(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
      close(12)      
c
  100 format(2i5,2f10.5)
  200 format(6i5)
  300 format(i5)
  400 format(5f12.7)
c
      return
      end
c---------------------------------------------------------- sub.finalo
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.finalo
      subroutine finalo(ffile)
      include 'common_vib8.fil'
c
c
      character*14 ffile
c
c
c
      open(12,file=ffile,status='unknown',form='unformatted')
       write(12)is,ie,js,je,ks,ke
       write(12)re,cs,dt,t,nt,iyper,izper,iter
       write(12)(((uu(i,j,k),i=is,ie+1),j=js,je-1),k=ks,ke-1)
       write(12)(((vv(i,j,k),i=is,ie),j=js,je),k=ks,ke-1)
       write(12)(((ww(i,j,k),i=is,ie),j=js,je-1),k=ks,ke)
       write(12)(((pp(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
c
       write(12)(((f(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(12)(((g(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(12)(((h(i,j,k),i=is,ie),j=js,je),k=ks,ke)
c
       write(12)(((tnu(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
      close(12)
c
      return
      end
c---------------------------------------------------------- sub.change
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.change
      subroutine change(dtmin)
      include 'common_vib8.fil'
c
      open(13,file='change.inf',status='unknown',form='formatted')
       read(13,100)nt
cc     read(13,101)dt,esp
      close(13)
c
      open(13,file='change.out',status='unknown',form='formatted')
       write(13,100)nt
       write(13,200)nn,t
       write(13,100)iter
       write(13,300)dt,dtmin,esp
       write(13,400)
       write(13,400)
       write(13,400)
       write(13,401)
       write(13,402)
       write(13,403)
       write(13,404)
       write(13,405)
c
      close(13)
  100 format(i10)
cc101 format(2f12.7)
  200 format(i10,f12.7)
  300 format(3f12.7)
  400  format(1h ,'  ')
  401  format(1h ,'c---------------------------- information')
  402  format(1h ,'nt           ;cycle number  ')
  403  format(1h ,'nn,t         ;present number , present time')
  404  format(1h ,'iter         ;iteration number  ')
  405  format(1h ,'dt,dtmin,esp ;dt , possible dt ,esp ')
c
c
      return
      end
c---------------------------------------------------------- sub.calstt
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.calstt
      subroutine calstt
      include 'common_vib8.fil'
c
      coef=2.0d0/3.0d0
      ck=0.094d0
      ck2=ck*ck
c
      call chguvw
c
      do 10 k=ks,ke
      do 10 j=js,je
      do 10 i=is,ie
       aveu(i,j,k)=aveu(i,j,k)+u(i,j,k)
       avev(i,j,k)=avev(i,j,k)+v(i,j,k)
       avew(i,j,k)=avew(i,j,k)+w(i,j,k)
       avet(i,j,k)=avet(i,j,k)+tnu(i,j,k)
   10 continue
c
      do 20 k=ks,ke-1
      do 20 j=js,je-1
      do 20 i=is,ie
       del=(dx(i+1)*dy(j+1)*dz(k+1))**0.33333
       rks=tnu(i,j,k)*tnu(i,j,k)/(ck2*del*del)
       avep(i,j,k)=avep(i,j,k)+(pp(i,j,k)-coef*rks)
       rmsp(i,j,k)=rmsp(i,j,k)
     &            +(pp(i,j,k)-coef*rks)*(pp(i,j,k)-coef*rks)
   20 continue
c
      do 25 k=ks,ke
      do 25 j=js,je
      do 25 i=is,ie
       del=(dxh(i)*dyh(j)*dzh(k))**0.33333
c
       tjp=(tnu(i,j,k)*dx(i)+tnu(i-1,j,k)*dx(i+1))/(dx(i)+dx(i+1))
       tjm=(tnu(i,j-1,k)*dx(i)+tnu(i-1,j-1,k)*dx(i+1))/(dx(i)+dx(i+1))
         tkp=(tjp*dy(j)+tjm*dy(j+1))/(dy(j)+dy(j+1))
c
       tjp2=(tnu(i,j,k-1)*dx(i)+tnu(i-1,j,k-1)*dx(i+1))/(dx(i)+dx(i+1))
       tjm2=(tnu(i,j-1,k-1)*dx(i)+tnu(i-1,j-1,k-1)*dx(i+1))
     &     /(dx(i)+dx(i+1))
         tkm=(tjp2*dy(j)+tjm2*dy(j+1))/(dy(j)+dy(j+1))
c
       tijk=(tkp*dz(k)+tkm*dz(k+1))/(dz(k)+dz(k+1))
c
       rks=tijk*tijk/(ck2*del*del)
c
       d11=(u(i+1,j,k)-u(i-1,j,k))/(dx(i+1)+dx(i))
       d22=(v(i,j+1,k)-v(i,j-1,k))/(dy(j+1)+dy(j))
       d33=(w(i,j,k+1)-w(i,j,k-1))/(dz(k+1)+dz(k))
c
       rs(i,j,k,7)=rs(i,j,k,7)+coef*rks-2.0d0*tijk*d11
       rs(i,j,k,8)=rs(i,j,k,8)+coef*rks-2.0d0*tijk*d22
       rs(i,j,k,9)=rs(i,j,k,9)+coef*rks-2.0d0*tijk*d33
c
   25 continue
c
      do 30 k=ks,ke
      do 30 j=js,je
      do 30 i=is,ie
       rs(i,j,k,1)=rs(i,j,k,1)+u(i,j,k)*u(i,j,k) 
       rs(i,j,k,2)=rs(i,j,k,2)+v(i,j,k)*v(i,j,k) 
       rs(i,j,k,3)=rs(i,j,k,3)+w(i,j,k)*w(i,j,k) 
       rs(i,j,k,4)=rs(i,j,k,4)+u(i,j,k)*v(i,j,k) 
       rs(i,j,k,5)=rs(i,j,k,5)+v(i,j,k)*w(i,j,k) 
       rs(i,j,k,6)=rs(i,j,k,6)+w(i,j,k)*u(i,j,k) 
   30 continue
c
c
      return
      end
c---------------------------------------------------------- sub.calfin
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.calfin
      subroutine calfin(istt)
      include 'common_vib8.fil'
c
c
       rstt=1.0d0/dble(istt)
       print *,'rstt =',rstt
c
      do 10 k=ks,ke
      do 10 j=js,je
      do 10 i=is,ie
       aveu(i,j,k)=aveu(i,j,k)*rstt
       avev(i,j,k)=avev(i,j,k)*rstt
       avew(i,j,k)=avew(i,j,k)*rstt
       avet(i,j,k)=avet(i,j,k)*rstt
   10 continue
      do 20 k=ks,ke
      do 20 j=js,je
      do 20 i=is,ie
       rs(i,j,k,1)=rs(i,j,k,1)*rstt-aveu(i,j,k)*aveu(i,j,k)
       rs(i,j,k,2)=rs(i,j,k,2)*rstt-avev(i,j,k)*avev(i,j,k)
       rs(i,j,k,3)=rs(i,j,k,3)*rstt-avew(i,j,k)*avew(i,j,k)
       rs(i,j,k,4)=rs(i,j,k,4)*rstt-aveu(i,j,k)*avev(i,j,k)
       rs(i,j,k,5)=rs(i,j,k,5)*rstt-avev(i,j,k)*avew(i,j,k)
       rs(i,j,k,6)=rs(i,j,k,6)*rstt-avew(i,j,k)*aveu(i,j,k)
c
       rs(i,j,k,7)=rs(i,j,k,7)*rstt
       rs(i,j,k,8)=rs(i,j,k,8)*rstt
       rs(i,j,k,9)=rs(i,j,k,9)*rstt
c
       rss(i,j,k,1)=rs(i,j,k,1)*rstt
       rss(i,j,k,2)=rs(i,j,k,2)*rstt
       rss(i,j,k,3)=rs(i,j,k,3)*rstt
   20 continue
c
      do 30 k=ks,ke-1
      do 30 j=js,je-1
      do 30 i=is,ie
       avep(i,j,k)=avep(i,j,k)*rstt
   30 continue
      do 40 k=ks,ke-1
      do 40 j=js,je-1
      do 40 i=is,ie
       rmsp(i,j,k)=rmsp(i,j,k)*rstt-avep(i,j,k)*avep(i,j,k)
   40 continue
c
c
      return
      end
c---------------------------------------------------------- sub.outstt
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.outstt
      subroutine outstt(sfile)
      include 'common_vib8.fil'
c
      character*14 sfile
c
c
      open(13,file=sfile,status='new',form='unformatted')
       write(13)is,ie,js,je,ks,ke
       write(13)re,cs,dt,t,nt,iyper,izper,iter
       write(13)(((aveu(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((avev(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((avew(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,1),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,2),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,3),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,4),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,5),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,6),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,7),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,8),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rs(i,j,k,9),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((avep(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
       write(13)(((rmsp(i,j,k),i=is,ie),j=js,je-1),k=ks,ke-1)
       write(13)(((avet(i,j,k),i=is,ie),j=js,je),k=ks,ke)
c
       write(13)(((rss(i,j,k,1),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rss(i,j,k,2),i=is,ie),j=js,je),k=ks,ke)
       write(13)(((rss(i,j,k,3),i=is,ie),j=js,je),k=ks,ke)
      close(13)
c
c
      return
      end
c---------------------------------------------------------- sub.calfrm
c                                                                    |
c        pfx  : force in x direction                                 |
c        pfy  : force in y direction                                 |
c        pitch: torsional moment                                     |
c        theta: degree of object                                     |
c                                                                    |
c        ndof: 0 not cal.                                            |
c              1 cal. x force                                        |
c              2 cal. y force                                        |
c              3 cal. torsional moment                               |
c              4 cal. x,y forces                                     |
c              5 cal. y force, torsional moment                      |
c              6 cal. 3 degree forces and moment(all)                |
c234567---------------------------------------------------- sub.calfrm 
      subroutine calfrm
      include 'common_vib8.fil'
c
       pfx=0.0d0
       pfy=0.0d0
       pitch=0.0d0
       theta=tdis
c      if(ndof.eq.0) go to 1000
c
c---------------------------------------------------cal.x,y forces
c      if(ndof.ne.3) then
        pfup=0.0d0
        pfdn=0.0d0
        pfrt=0.0d0
        pflt=0.0d0
       do 100 k=ks,ke-1
       do 100 j=js,je-1
       do 100 i=is,ie-1
       if((ibf(i,j,k).eq.0).and.(ibf(i,j-1,k).eq.1)) then
         pfup=pfup+pp(i,j,k)*dx(i+1)*dz(k+1)
        else
	 pfup=pfup
       end if
       if((ibf(i,j,k).eq.0).and.(ibf(i,j+1,k).eq.1)) then
	 pfdn=pfdn+pp(i,j,k)*dx(i+1)*dz(k+1)
	else
	 pfdn=pfdn
       end if
       if((ibf(i,j,k).eq.0).and.(ibf(i+1,j,k).eq.1)) then
         pflt=pflt+pp(i,j,k)*dy(j+1)*dz(k+1)
        else
         pflt=pflt
       end if
       if((ibf(i,j,k).eq.0).and.(ibf(i-1,j,k).eq.1)) then
         pfrt=pfrt+pp(i,j,k)*dy(j+1)*dz(k+1)
       else
         pfrt=pfrt
       end if
c
c**************************************************************
       pfx=pflt-pfrt
       pfy=pfdn-pfup
       pfd=pfx*dcos(theta)-pfy*dsin(theta)
       pfl=pfx*dsin(theta)+pfy*dcos(theta)
c**************************************************************
c
  100 continue
c      end if
c      if((ndof.eq.1).or.(ndof.eq.5)) pfx=0.0d0
c      if (ndof.eq.2)                 pfy=0.0d0
c
c
c---------------------------------------------------cal. moment
c      if((ndof.eq.3).or.(ndof.eq.5).or.(ndof.eq.6)) then
        pitup=0.0d0
        pitdn=0.0d0
       do 200 k=ks,ke-1
       do 200 j=js,je-1
       do 200 i=is,ie-1
       if((ibf(i,j,k).eq.0).and.(ibf(i,j-1,k).eq.1)) then
         pitup=pitup+0.5d0*pp(i,j,k)*(x(i)+x(i+1))*dx(i+1)*dz(k+1)
        else if((ibf(i,j,k).eq.0).and.(ibf(i+1,j,k).eq.1)) then
         pitup=pitup+0.5d0*pp(i,j,k)*(y(j)+y(j+1))*dy(j+1)*dz(k+1)
        else
         pitup=pitup
       end if
       if((ibf(i,j,k).eq.0).and.(ibf(i,j+1,k).eq.1)) then
         pitdn=pitdn+0.5d0*pp(i,j,k)*(x(i)+x(i+1))*dx(i+1)*dz(k+1)
        else if((ibf(i,j,k).eq.0).and.(ibf(i-1,j,k).eq.1)) then
         pitdn=pitdn+0.5d0*pp(i,j,k)*(y(j)+y(j+1))*dy(j+1)*dz(k+1)
        else
         pitdn=pitdn
       end if
       pitch=pitup-pitdn
  200 continue
c
      pft=pitch
c
c      end if
c
c 1000 continue
c
c      frmt(1,nn)=pfx
c      frmt(2,nn)=pfy
c      frmt(3,nn)=pitch
c
      write(*,111)nn,pfx,pfy,pft
111   format('###   in calfrm nn,pfx,pfy,pft =',i5,3f10.5)  
c
      return
      end
c---------------------------------------------------------- sub.calvib
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.calvib
      subroutine calvib
      include 'common_vib8.fil'
c
      pai=3.1415926d0
      xomg=dsqrt(xkst/amas)
      yomg=dsqrt(ykst/amas)
      tomg=dsqrt(tkst/aitm)
c
      if(ndof.eq.0) go to 1000
c
      if((ndof.eq.2).or.(ndof.eq.4).or.(ndof.eq.6)) then
      xdis0=xdis
      xvel0=xvel
      xacc0=xacc
c
      xup=pfx-xcvd*(xvel0+0.5d0*xacc0*dt)
     &   -xkst*(xdis0+xvel0*dt+(0.5d0-beta)*xacc0*dt*dt)
      xdn=amas+0.5d0*xcvd*dt+beta*xkst*dt*dt
      xacc=xup/xdn
      xvel=xvel0+0.5d0*(xacc+xacc0)*dt
      xdis=xdis0+xvel0*dt+(0.5d0-beta)*xacc0*dt*dt+beta*xacc*dt*dt
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc add. forced disp.
      if(ampxx.ne.0.0) then
      if(t.lt.tcntn) then
       xdis=ampxx*dsin(xomg*t)
       xvel=ampxx*xomg*dcos(xomg*t)
       xacc=-ampxx*xomg*xomg*dsin(xomg*t)
      end if
      end if 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      else
      xacc=0.0d0
      xvel=0.0d0
      xdis=0.0d0
      endif
c
      if((ndof.eq.1).or.(ndof.eq.4).or.(ndof.eq.5).or.(ndof.eq.6)) then
      ydis0=ydis
      yvel0=yvel
      yacc0=yacc
c
      yup=pfy-ycvd*(yvel0+0.5d0*yacc0*dt)
     &   -ykst*(ydis0+yvel0*dt+(0.5d0-beta)*yacc0*dt*dt)
      ydn=amas+0.5d0*ycvd*dt+beta*ykst*dt*dt
      yacc=yup/ydn
      yvel=yvel0+0.5d0*(yacc+yacc0)*dt
      ydis=ydis0+yvel0*dt+(0.5d0-beta)*yacc0*dt*dt+beta*yacc*dt*dt
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc add. forced disp.
      if(ampyy.ne.0.0) then
      if(t.lt.tcntn) then
       ydis=ampyy*dsin(yomg*t)
       yvel=ampyy*yomg*dcos(yomg*t)
       yacc=-ampyy*yomg*yomg*dsin(yomg*t)
      end if
      end if 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      else
      yacc=0.0d0
      yvel=0.0d0
      ydis=0.0d0
      endif
c
c
      if((ndof.eq.3).or.(ndof.eq.5).or.(ndof.eq.6)) then
      tdis0=tdis
      tvel0=tvel
      tacc0=tacc
c
      tup=pft-tcvd*(tvel0+0.5d0*tacc0*dt)
     &   -tkst*(tdis0+tvel0*dt+(0.5d0-beta)*tacc0*dt*dt)
      tdn=aitm+0.5d0*tcvd*dt+beta*tkst*dt*dt
      tacc=tup/tdn
      tvel=tvel0+0.5d0*(tacc+tacc0)*dt
      tdis=tdis0+tvel0*dt+(0.5d0-beta)*tacc0*dt*dt+beta*tacc*dt*dt
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc add. forced disp.
      if(amptt.ne.0.0) then
      if(t.lt.tcntn) then
       tdis=amptt*dsin(tomg*t)
       tvel=amptt*tomg*dcos(tomg*t)
       tacc=-amptt*tomg*tomg*dsin(tomg*t)
      end if
      end if 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      else
      tacc=0.0d0
      tvel=0.0d0
      tdis=0.0d0
      endif
c
c
 1000 continue
c
      write(*,111)nn,pfx,pfy,pft
111   format('###   in calvib nn,pfx,pfy,pft =',i5,3f10.5)  
c
      if(nn.eq.ndivf1) then
       ndivf1=ndivf1+ndivf
       nhcnt=nhcnt+1
       frmt(1,nhcnt)=pfx
       frmt(2,nhcnt)=pfy
       frmt(3,nhcnt)=pft
c
      resp(1,1,nhcnt)=xacc
      resp(1,2,nhcnt)=xvel
      resp(1,3,nhcnt)=xdis
c
      resp(2,1,nhcnt)=yacc
      resp(2,2,nhcnt)=yvel
      resp(2,3,nhcnt)=ydis
c
      resp(3,1,nhcnt)=tacc
      resp(3,2,nhcnt)=tvel
      resp(3,3,nhcnt)=tdis
      endif
c
      if(nn.eq.1) then
       xvels=xvel
       yvels=yvel
       tvels=tvel
       xdiss=xdis
       ydiss=ydis
       tdiss=tdis
      end if
c
      do 400 k=ks,ke-1
      do 400 j=js,je-1
      do 400 i=is,ie
      vu=0.5d0*((dx(i)*(vv(i,j,k)+vv(i,j+1,k))
     &        +dx(i+1)*(vv(i-1,j,k)+vv(i-1,j+1,k)))/(dx(i)+dx(i+1)))
      if(i.eq.is) vu=0.5d0*(vv(i  ,j,k)+vv(i  ,j+1,k))
      if(i.eq.ie) vu=0.5d0*(vv(i-1,j,k)+vv(i-1,j+1,k))
c
       facc0(i,j,k)=facc1(i,j,k)
       yyy=(y(j)+y(j+1))*0.5d0
       facc1(i,j,k)=-2.0d0*tvel*vu-tacc*yyy-tvel*tvel*x(i)
c
c********************************** non I
ccc       facc1(i,j,k)=0.0d0
c********************************** only inertial force
ccc       facc1(i,j,k)=-tacc*yyy
c********************************** only coriolis force
ccc       facc1(i,j,k)=-2.0d0*tvel*vu
c********************************** only centrifugal force
ccc       facc1(i,j,k)=-tvel*tvel*x(i)
c
  400 continue
c
      do 450 k=ks,ke-1
      do 450 j=js,je
      do 450 i=is,ie-1
      uv=0.5d0*((dy(j)*(uu(i,j,k)+uu(i+1,j,k))
     &        +dy(j+1)*(uu(i,j-1,k)+uu(i+1,j-1,k)))/(dy(j)+dy(j+1)))
      if(j.eq.js) uv=0.5d0*(uu(i,j  ,k)+uu(i+1,j  ,k))
      if(j.eq.je) uv=0.5d0*(uu(i,j-1,k)+uu(i+1,j-1,k))
c
       gacc0(i,j,k)=gacc1(i,j,k)
       xxx=(x(i)+x(i+1))*0.5d0
       gtheta=datan2(y(j),xxx)
       gacc1(i,j,k)=2.0d0*tvel*uv+tacc*xxx-tvel*tvel*y(j)
c
c********************************** non I
ccc       gacc1(i,j,k)=0.0d0
c********************************** only inertial force
ccc       gacc1(i,j,k)=tacc*xxx
c********************************** only coriolis force
ccc       gacc1(i,j,k)=2.0d0*tvel*uv
c********************************** only centrifugal force
ccc       gacc1(i,j,k)=-tvel*tvel*y(j)
c
  450 continue
c 
c
      return
      end
c---------------------------------------------------------- sub.outvib
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c                                                                    |
c234567---------------------------------------------------- sub.outvib
      subroutine outvib(vfile,nfile)
      include 'common_vib8.fil'
c
      character*14 vfile,nfile
c
      print *,'vfile=',vfile
      print *,'nfile=',nfile
      dt2=dt*real(ndivf)
C 8/20/13 back to original 8/16/13 text out
      open(unit=11,file=vfile,status='unknown',form='unformatted')
      write(11)nhcnt,dt2,beta
      write(11)xdiss,ydiss,tdiss
      write(11)xvels,yvels,tvels
      write(11)((frmt(i,n),i=1,3),n=1,nhcnt)
      write(11)((resp(1,i,n),i=1,3),n=1,nhcnt)
      write(11)((resp(2,i,n),i=1,3),n=1,nhcnt)
      write(11)((resp(3,i,n),i=1,3),n=1,nhcnt)
C      open(unit=11,file=vfile,status='unknown',form='formatted')
C      write(11,'(i10,2e12.4)')nhcnt,dt2,beta
C      write(11,'(5e12.5)')xdiss,ydiss,tdiss
C      write(11,'(5e12.5)')xvels,yvels,tvels
C      write(11,'(5e12.5)')((frmt(i,n),i=1,3),n=1,nhcnt)
C      write(11,'(5e12.5)')((resp(1,i,n),i=1,3),n=1,nhcnt)
C      write(11,'(5e12.5)')((resp(2,i,n),i=1,3),n=1,nhcnt)
C      write(11,'(5e12.5)')((resp(3,i,n),i=1,3),n=1,nhcnt)
      close(11)
c
      open(unit=12,file=nfile,status='unknown',form='unformatted')
       write(12)is,ie,js,ks,ke
       write(12)xacc,yacc,tacc
       write(12)xdis,ydis,tdis
       write(12)xvel,yvel,tvel
       write(12)(((facc1(i,j,k),i=is,ie),j=js,je),k=ks,ke)
       write(12)(((gacc1(i,j,k),i=is,ie),j=js,je),k=ks,ke)
      close(12)
c
c
      return
      end
