c     *************************************************************
c     ** three-dimensional incompressible viscous flow solved by **
c     **                                                         **
c     ** finite difference method for rectilinear coordinate     **
c     **                                                         **
c     ** 2nd order central scheme,  mac staggered grid           **
c     **                                                         **
c     ** fractional step method by kim and moin                  **
c     **                                                         **
c     **                                                         **
c     **                                                         **
c     **  unit 1:grid data reading                               **
c     **  unit 2:pre-result reading                              **
c     **  unit 3:result output                                   **
c     **                                                         **
c     **  1999/07/24  version 3.0    interpolation scheme        **
c     **                                                         **
c     **                                                         **
c     **                                       taro Nihonn       **
c     *************************************************************
c
c     ******************
c     ** main routine **
c     ******************
c
c
         call init
c
         call grid
c
         call solve
c
         call output
c
c
      stop
      end
c
c     ********************************************
c     ** solve the flow by time marching method **
c     ********************************************
c
      subroutine solve
c
      implicit double precision(a-h,o-z)
      common /tmcnd/ dt,nstep
      common /past/ time
c
c
         call boundp
c
         call boundv
c
c     ***** implicit time marching loop *****
c
      do 500 istep=1,nstep
c
         time=time+dt
c
      if(mod(istep,10).eq.0) then
      write(*,600) time
  600 format(' time=',f10.4)
      end if
c
         call convx
c
         call convy
c
         call convz
c
         call convec
c
         call visc
c
         call explicit
c
         call boundv
c
         call poiss(istep)
c
         call boundp
c
         call euler
c
       if(mod(istep,10).eq.0) then
         call judge(errs)
        if(errs.le.1.d-6) goto 10
       end if
c
c
  500 continue
c
      write(*,*)' not converged!!'
c
   10 continue
c
c
      return
      end
c
c     ***************************
c     ** setting constant data **
c     ***************************
c
      subroutine init
c
      implicit double precision(a-h,o-z)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /conv1/ eps1,iter1
      common /tmcnd/ dt,nstep
      common /reyno/ re,rei
      common /cont/ icon
      common /choice/ iselc
c
c     ***** iteration parameter *****
c
         iter1=50
c
c     **** convergence parameter ****
c
         eps1=1.0d-3
c
c     *** to be continue yes=1 no=0 ***
c
         icon=0
c
c     **** select convective term ****
c     **** 0=2nd  1=consistent ****
c
         iselc=1
c
c     ****** time step *****
c
         dt=0.002d0
c
c     ***** reynolds number *****
c
         re =1000.d0
         rei=1.d0/re
c
c     ***** time iteration *****
c
         nstep=10000
c
c
      return
      end
c
c     **************************
c     ** subroutine grid data **
c     **************************
c
      subroutine grid
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /oldvel/ vo(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
      common /cont/ icon
      common /past/ time
c
c
        nx=21
        ny=21
        nz=21
c
      dx=1.d0/dble(nx-1)
      dy=1.d0/dble(ny-1)
      dz=1.d0/dble(nz-1)
c
      dxi=1.d0/dx
      dyi=1.d0/dy
      dzi=1.d0/dz
      dx2i=1.d0/dx**2
      dy2i=1.d0/dy**2
      dz2i=1.d0/dz**2
c
c     ***** number of divide *****
c
         nxm1=nx-1
         nym1=ny-1
         nzm1=nz-1
         nxp1=nx+1
         nyp1=ny+1
         nzp1=nz+1
c
c     ****  impulsive start  ****
c
      do 110 k=1,nz
      do 110 j=1,ny
      do 110 i=1,nx
c
         p(i,j,k)   =0.d0
         v(i,j,k,1) =0.d0
         v(i,j,k,2) =0.d0
  110    v(i,j,k,3) =0.d0
c
      if(icon.eq.1) then
c
c     *** pre-results reading ***
c
      open(2,file='rd1.dat')
c
      read(2,510) nx,ny,nz,time,dmy
      read(2,520)((((v(i,j,k,l),l=1,3),p(i,j,k),i=1,nx)
     &           ,j=1,ny),k=1,nz)
c
  510 format(3i10,f10.5,f10.0)
  520 format(3f10.6,e15.7)
c
      end if
c
c
      return
      end
c
c     ****************************
c     ** output result to file  **
c     ****************************
c
      subroutine output
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /tmcnd/ dt,nstep
      common /past/ time
      common /reyno/ re,rei
c
c
c     *** output 3d data ***
c
      open(3,file='rd1.dat')
c
      write(3,600) nx,ny,nz,time,re
      write(3,610)((((v(i,j,k,l),l=1,3),p(i,j,k),i=1,nx)
     &            ,j=1,ny),k=1,nz)
c
      close(3)
c
c     *** output 3d data interpolation for regular position ***
c
      open(20,file='rr1.dat')
c
      write(20,600) nx,ny,nz,time,re
c
      do 100 k=1,nz
      do 100 j=1,ny
      do 100 i=1,nx
c
      vre1=0.25d0*(v(i,j,k  ,1)+v(i,j+1,k  ,1)
     &            +v(i,j,k+1,1)+v(i,j+1,k+1,1))
c
      vre2=0.25d0*(v(i  ,j,k,2)+v(i  ,j,k+1,2)
     &            +v(i+1,j,k,2)+v(i+1,j,k+1,2))
c
      vre3=0.25d0*(v(i  ,j,k,3)+v(i  ,j+1,k,3)
     &            +v(i+1,j,k,3)+v(i+1,j+1,k,3))
c
      pre=0.125d0*
     &    (p(i  ,j,k)+p(i  ,j,k+1)+p(i  ,j+1,k)+p(i  ,j+1,k+1)
     &    +p(i+1,j,k)+p(i+1,j,k+1)+p(i+1,j+1,k)+p(i+1,j+1,k+1))
c
      write(20,610) vre1,vre2,vre3,pre
c
  100 continue
c
      close(20)
c
  600 format(3i10,f10.5,f10.0)
  610 format(3f10.6,e15.7)
c
c
      return
      end
c
c     *******************************
c     ** judgement of steady state **
c     *******************************
c
      subroutine judge(errs)
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /oldvel/ vo(mx,my,mz,3)
      dimension errv(3)
c
c     **** judge for steady state ****
c
         errs=0.d0
         errv(1)=0.d0
         errv(2)=0.d0
         errv(3)=0.d0
c
c     *** maximum error ***
c
      do 100 l=1,3
      do 100 k=2,nzm1
      do 100 j=2,nym1
      do 100 i=2,nxm1
c
         ersq=(v(i,j,k,1)-vo(i,j,k,1))**2
     &       +(v(i,j,k,2)-vo(i,j,k,2))**2
     &       +(v(i,j,k,3)-vo(i,j,k,3))**2
c
         errs=errs+ersq
c
         if(dabs(v(i,j,k,l)).gt.errv(l)) then 
           errv(l)=dabs(v(i,j,k,l))
         end if
c
       if(errv(l).gt.5.d0) then
         write(*,*)'***************************'
         write(*,*)'*** faulty terminated ! ***'
         write(*,*)'***************************'
         goto 99
       end if
c
  100 continue
c
c
         errs=dsqrt(errs/dfloat((nx-2)*(ny-2)*(nz-2)))
c
      write(*,*)'     '
      write(*,*)'**********************************************'
      write(*,660) errs
      do 120 l=1,3
  120 write(*,670) l,errv(l)
      write(*,*)'**********************************************'
      write(*,*)'     '
c
  660 format('***  rms check=',e10.4)
  670 format('***  max. velocity check=',i3,2x,f8.4)
c
c
      return
c
  99   call output
c
c
       stop
       end
c
c     *****************************************
c     ** setting velocity boundary condition **
c     *****************************************
c
      subroutine boundv
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
c
c     **** no-slip condition ****
c
      do 100 j=1,ny
      do 100 i=1,nx
c
         v(i,j,nz+1,1)=-v(i,j,nz,1)
         v(i,j,nz+1,2)=-v(i,j,nz,2)
         v(i,j,nz,3)=0.d0
c
         v(i,j,1 ,1)=-v(i,j,2 ,1)
         v(i,j,1 ,2)=-v(i,j,2 ,2)
         v(i,j,1 ,3)=0.d0
c
  100 continue
c
c      **** top wall condition ****
c
      do 120 k=1,nz
      do 120 i=1,nx
c
         v(i,1 ,k,1)=-v(i,2 ,k,1)
         v(i,1 ,k,2)=0.d0
         v(i,1 ,k,3)=-v(i,2 ,k,3)
c
         v(i,ny+1,k,1)=2.d0-v(i,ny,k,1)
         v(i,ny  ,k,2)=0.d0
         v(i,ny+1,k,3)=-v(i,ny,k,3)
c
  120 continue
c
c     **** wall condition ****
c
      do 130 k=1,nz
      do 130 j=1,ny
c
         v(1 ,j,k,1)=0.d0
         v(1 ,j,k,2)=-v(2 ,j,k,2)
         v(1 ,j,k,3)=-v(2 ,j,k,3)
c
         v(nx,j,k,1)=0.d0
         v(nx+1,j,k,2)=-v(nx,j,k,2)
         v(nx+1,j,k,3)=-v(nx,j,k,3)
c
  130 continue
c
c
      return
      end
c
c     *****************************************
c     ** setting pressure boundary condition **
c     *****************************************
c
      subroutine boundp
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /diffp/ dp,dpdx
c
c     **** on the wall condition ****
c
      do 100 j=1,ny+1
      do 100 i=1,nx+1
c
          p(i,j,1   )=(4.d0*p(i,j,2   )-p(i,j,3   ))/3.d0
          p(i,j,nzp1)=(4.d0*p(i,j,nz  )-p(i,j,nzm1))/3.d0
c
  100 continue
c
c     **** on the wall condition ****
c
      do 110 k=1,nz+1
      do 110 i=1,nx+1
c
          p(i,1   ,k)=(4.d0*p(i,2   ,k)-p(i,3   ,k))/3.d0
          p(i,nyp1,k)=(4.d0*p(i,ny  ,k)-p(i,nym1,k))/3.d0
c
  110 continue
c
c     **** on the wall condition ****
c
      do 120 k=1,nz+1
      do 120 j=1,ny+1
c
          p(1   ,j,k)=(4.d0*p(2   ,j,k)-p(3   ,j,k))/3.d0
          p(nxp1,j,k)=(4.d0*p(nx  ,j,k)-p(nxm1,j,k))/3.d0
c
  120 continue
c
c
      return
      end
c
c     *****************************
c     ** compute convective term **
c     *****************************
c
      subroutine convx
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /conv/ cx(mx,my,mz,3),cy(mx,my,mz,3),cz(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
c
c     *** compute convective & viscous term ***
c
      do 100 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 100 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 100 i=2,nxp1
c
         ip1=i+1
         im1=i-1
c
      cx(i,j,k,1)=0.5d0*(v(i,j,k,1)+v(im1,j,k,1))
     &                 *(v(i,j,k,1)-v(im1,j,k,1))*dxi
c
  100 continue
c
      do 110 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      cz(i,j,k,1)=0.5d0*(v(ip1,j,k,3)+v(i,j,k,3))
     &                 *(v(i,j,kp1,1)-v(i,j,k,1))*dzi
c
  110 continue
c
c
      do 120 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 120 j=2,nyp1
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      cy(i,j,k,1)=0.5d0*(v(ip1,j,k,2)+v(i,j,k,2))
     &                 *(v(i,jp1,k,1)-v(i,j,k,1))*dyi
c
  120 continue
c
c
      return
      end
c
c     *****************************
c     ** compute convective term **
c     *****************************
c
      subroutine convy
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /conv/ cx(mx,my,mz,3),cy(mx,my,mz,3),cz(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
c
c     *** compute convective & viscous term ***
c
      do 100 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 100 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 100 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      cx(i,j,k,2)=0.5d0*(v(i,jp1,k,1)+v(i,j,k,1))
     &                 *(v(ip1,j,k,2)-v(i,j,k,2))*dxi
c
  100 continue
c
c
      do 110 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      cz(i,j,k,2)=0.5d0*(v(i,jp1,k,3)+v(i,j,k,3))
     &                 *(v(i,j,kp1,2)-v(i,j,k,2))*dzi
c
c
  110 continue
c
      do 120 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 120 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      cy(i,j,k,2)=0.5d0*(v(i,j,k,2)+v(i,jm1,k,2))
     &                 *(v(i,j,k,2)-v(i,jm1,k,2))*dyi
c
  120 continue
c
c
      return
      end
c
c     *****************************
c     ** compute convective term **
c     *****************************
c
      subroutine convz
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /conv/ cx(mx,my,mz,3),cy(mx,my,mz,3),cz(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
c
c     *** compute convective & viscous term ***
c
      do 100 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 100 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 100 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      cx(i,j,k,3)=0.5d0*(v(i,j,kp1,1)+v(i,j,k,1))
     &                 *(v(ip1,j,k,3)-v(i,j,k,3))*dxi
c
  100 continue
c
c
      do 110 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      cz(i,j,k,3)=0.5d0*(v(i,j,k,3)+v(i,j,km1,3))
     &                 *(v(i,j,k,3)-v(i,j,km1,3))*dzi
c
  110 continue
c
      do 120 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 120 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      cy(i,j,k,3)=0.5d0*(v(i,j,kp1,2)+v(i,j,k,2))
     &                 *(v(i,jp1,k,3)-v(i,j,k,3))*dyi
c
  120 continue
c
c
      return
      end
c
c     *****************************
c     ** compute convective term **
c     *****************************
c
      subroutine convec
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /aright/ h1(mx,my,mz,3),h2(mx,my,mz,3)
      common /conv/ cx(mx,my,mz,3),cy(mx,my,mz,3),cz(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
      common /choice/ iselc
c
c     *** compute convective term ***
c
      do 100 l=1,3
      do 100 k=2,nz
      do 100 j=2,ny
      do 100 i=2,nx
c
  100    h2(i,j,k,l)=h1(i,j,k,l)
c
c
      if(iselc.eq.0) then
c
      do 120 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 120 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c     *****  convective velocity  *****
c
      cux=v(i,j,k,1)
      cvx=0.25d0*(v(i,j,k,2)+v(i,jm1,k,2)+v(ip1,j,k,2)+v(ip1,jm1,k,2))
      cwx=0.25d0*(v(i,j,k,3)+v(i,j,km1,3)+v(ip1,j,k,3)+v(ip1,j,km1,3))
c
      cuy=0.25d0*(v(i,j,k,1)+v(i,jp1,k,1)+v(im1,j,k,1)+v(im1,jp1,k,1))
      cvy=v(i,j,k,2)
      cwy=0.25d0*(v(i,j,k,3)+v(i,jp1,k,3)+v(i,j,km1,3)+v(i,jp1,km1,3))
c
      cuz=0.25d0*(v(i,j,k,1)+v(i,j,kp1,1)+v(im1,j,k,1)+v(im1,j,kp1,1))
      cvz=0.25d0*(v(i,j,k,2)+v(i,j,kp1,2)+v(i,jm1,k,2)+v(i,jm1,kp1,2))
      cwz=v(i,j,k,3)
c
c     *****  convective term  *****
c
      u1=.5d0*(v(ip1,j,k,1)-v(im1,j,k,1))*dxi
      u3=.5d0*(v(i,jp1,k,1)-v(i,jm1,k,1))*dyi
      u5=.5d0*(v(i,j,kp1,1)-v(i,j,km1,1))*dzi
c
      v1=.5d0*(v(ip1,j,k,2)-v(im1,j,k,2))*dxi
      v3=.5d0*(v(i,jp1,k,2)-v(i,jm1,k,2))*dyi
      v5=.5d0*(v(i,j,kp1,2)-v(i,j,km1,2))*dzi
c
      w1=.5d0*(v(ip1,j,k,3)-v(im1,j,k,3))*dxi
      w3=.5d0*(v(i,jp1,k,3)-v(i,jm1,k,3))*dyi
      w5=.5d0*(v(i,j,kp1,3)-v(i,j,km1,3))*dzi
c
      conv(i,j,k,1)=cux*u1+cvx*u3+cwx*u5
      conv(i,j,k,2)=cuy*v1+cvy*v3+cwy*v5
      conv(i,j,k,3)=cuz*w1+cvz*w3+cwz*w5
c
  120 continue
c
      else if(iselc.eq.1) then
c
      do 130 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 130 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 130 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c
      conv(i,j,k,1)=(cx(i,j,k,1)+cx(i+1,j,k,1)+cy(i,j,k,1)+cy(i,j-1,k,1)
     &              +cz(i,j,k,1)+cz(i,j,k-1,1))*.5d0
c
  130 continue
c
c
      do 140 k=2,nzm1
c
         kp1=k+1
         km1=k-1
c
      do 140 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 140 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      conv(i,j,k,3)=(cx(i,j,k,3)+cx(i-1,j,k,3)+cy(i,j,k,3)+cy(i,j-1,k,3)
     &              +cz(i,j,k,3)+cz(i,j,k+1,3))*.5d0
c
  140 continue
c
      do 150 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 150 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 150 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      conv(i,j,k,2)=(cx(i,j,k,2)+cx(i-1,j,k,2)+cy(i,j,k,2)+cy(i,j+1,k,2)
     &              +cz(i,j,k,2)+cz(i,j,k-1,2))*.5d0
c
  150 continue
c
      end if
c
c
      return
      end
c
c     **************************
c     ** compute viscous term **
c     **************************
c
      subroutine visc
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /reyno/ re,rei
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
c
c     *** compute viscous term ***
c
      do 100 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 100 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 100 i=2,nxm1
c
         ip1=i+1
         im1=i-1
c
c     *****  viscous term  *****
c
      vx1=(v(ip1,j,k,1)-2.d0*v(i,j,k,1)+v(im1,j,k,1))*dx2i
     &   +(v(i,jp1,k,1)-2.d0*v(i,j,k,1)+v(i,jm1,k,1))*dy2i
     &   +(v(i,j,kp1,1)-2.d0*v(i,j,k,1)+v(i,j,km1,1))*dz2i
c
        vis(i,j,k,1)=vx1*rei
c
  100 continue
c
      do 110 k=2,nzm1
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c     *****  viscous term  *****
c
      vz1=(v(ip1,j,k,3)-2.d0*v(i,j,k,3)+v(im1,j,k,3))*dx2i
     &   +(v(i,jp1,k,3)-2.d0*v(i,j,k,3)+v(i,jm1,k,3))*dy2i
     &   +(v(i,j,kp1,3)-2.d0*v(i,j,k,3)+v(i,j,km1,3))*dz2i
c
        vis(i,j,k,3)=vz1*rei
c
  110 continue
c
      do 120 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 120 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 120 i=2,nx
c
         ip1=i+1
         im1=i-1
c
c     *****  viscous term  *****
c
      vy1=(v(ip1,j,k,2)-2.d0*v(i,j,k,2)+v(im1,j,k,2))*dx2i
     &   +(v(i,jp1,k,2)-2.d0*v(i,j,k,2)+v(i,jm1,k,2))*dy2i
     &   +(v(i,j,kp1,2)-2.d0*v(i,j,k,2)+v(i,j,km1,2))*dz2i
c
        vis(i,j,k,2)=vy1*rei
c
  120 continue
c
c
      return
      end
c
c     *******************************
c     ** explicit time integration **
c     *******************************
c
      subroutine explicit
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /tmcnd/ dt,nstep
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /metric/ xi(mx),et(my),ze(mz)
      common /oldvel/ vo(mx,my,mz,3)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
c
c     *** predictor step by euler forward method ***
c
c     *** velocity interchange ***
c
      do 100 k=2,nz
      do 100 j=2,ny
      do 100 i=2,nx
c
         vo(i,j,k,1)=v(i,j,k,1)
         vo(i,j,k,2)=v(i,j,k,2)
  100    vo(i,j,k,3)=v(i,j,k,3)
c
c     *** tentative velocity ***
c
      do 110 k=2,nz
      do 110 j=2,ny
      do 110 i=2,nxm1
c
      v(i,j,k,1)=vo(i,j,k,1)+dt*(-conv(i,j,k,1)+vis(i,j,k,1))
c
  110 continue
c
      do 120 k=2,nzm1
      do 120 j=2,ny
      do 120 i=2,nx
c
      v(i,j,k,3)=vo(i,j,k,3)+dt*(-conv(i,j,k,3)+vis(i,j,k,3))
c
  120 continue
c
      do 130 k=2,nz
      do 130 j=2,nym1
      do 130 i=2,nx
c
      v(i,j,k,2)=vo(i,j,k,2)+dt*(-conv(i,j,k,2)+vis(i,j,k,2))
c
  130 continue
c
c
      return
      end
c
c     **********************************
c     ** poisson equation of pressure **
c     **********************************
c
      subroutine poiss(is)
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /conv1/ eps1,iter1
      common /tmcnd/ dt,nstep
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
      dimension rhs(mx,my,mz)
c
c     *** right hand side of poisson eqation ***
c
      do 100 k=2,nz
c
         km1=k-1
c
      do 100 j=2,ny
c
         jm1=j-1
c
      do 100 i=2,nx
c
         im1=i-1
c
      ux=(v(i,j,k,1)-v(im1,j,k,1))*dxi
      vy=(v(i,j,k,2)-v(i,jm1,k,2))*dyi
      wz=(v(i,j,k,3)-v(i,j,km1,3))*dzi
c
        div=ux+vy+wz
c
        rhs(i,j,k)=-div/dt
c
  100 continue
c
c     ***** sor iteration loop *****
c
      do 500 loop=1,iter1
c
         err=0.d0
c
         call boundp
c
c     ***** cal. internal point *****
c
      do 130 k=2,nz
c
         kp1=k+1
         km1=k-1
c
      do 130 j=2,ny
c
         jp1=j+1
         jm1=j-1
c
      do 130 i=2,nx
c
         ip1=i+1
         im1=i-1
c
      rp=(p(ip1,j,k)+p(im1,j,k))*dx2i
     &  +(p(i,jp1,k)+p(i,jm1,k))*dy2i
     &  +(p(i,j,kp1)+p(i,j,km1))*dz2i+rhs(i,j,k)
c
c     **** cal. residual ****
c
        dp=.5d0*rp/(dx2i+dy2i+dz2i)-p(i,j,k)
c
        p(i,j,k)=p(i,j,k)+1.8d0*dp
c
        err=err+dp*dp
c
  130 continue
c
        err=dsqrt(err/dfloat((nxm1-1)*(nym1-1)*(nzm1-1)))
        if(err.lt.eps1) goto 90
c
  500 continue
c
   90 continue
c
       if(mod(is,10).eq.0) then
      write(*,600) loop,err
  600 format(' poisson iteration & error=',i5,2x,e12.5)
       end if
c
c
      return
      end
c
c     ***************************************
c     ** corrector stepof time integration **
c     ***************************************
c
      subroutine euler
c
      implicit double precision(a-h,o-z)
      parameter (mx=51,my=51,mz=51)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /tmcnd/ dt,nstep
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(0:mx,0:my,0:mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
      common /space/ dxi,dyi,dzi,dx2i,dy2i,dz2i
c
c     *** time integration by euler forward method ***
c
      do 100 k=1,nz
      do 100 j=1,ny
      do 100 i=1,nxm1
c
c     *****  pressure term  *****
c
      px=(p(i+1,j,k)-p(i,j,k))*dxi
c
c     *** new velocity ***
c
      v(i,j,k,1)=v(i,j,k,1)-dt*px
c
  100 continue
c
c
      do 110 k=1,nzm1
      do 110 j=1,ny
      do 110 i=1,nx
c
c     *****  pressure term  *****
c
      pz=(p(i,j,k+1)-p(i,j,k))*dzi
c
c     *** new velocity ***
c
      v(i,j,k,3)=v(i,j,k,3)-dt*pz
c
  110 continue
c
      do 120 k=1,nz
      do 120 j=1,nym1
      do 120 i=1,nx
c
c     *****  pressure term  *****
c
      py=(p(i,j+1,k)-p(i,j,k))*dyi
c
c     *** new velocity ***
c
      v(i,j,k,2)=v(i,j,k,2)-dt*py
c
  120 continue
c
c
      return
      end