c     *************************************************************
c     ** three-dimensional incompressible viscous flow solved by **
c     **                                                         **
c     ** finite difference method for rectilinear coordinate     **
c     **                                                         **
c     ** 2nd order central scheme                                **
c     **                                                         **
c     ** regular grid                                            **
c     **                                                         **
c     ** euler forward scheme                                    **
c     **                                                         **
c     **  unit 1:grid data reading                               **
c     **  unit 2:pre-result reading                              **
c     **  unit 3:result output                                   **
c     **                                                         **
c     **                                                         **
c     **   さぁ，頑張って計算しよう！                            **
c     **                                                         **
c     **                                                         **
c     **  2005/06/12  version 1.0                                **
c     **                                         name            **
c     **                                                         **
c     *************************************************************
c
c     ******************
c     ** main routine **
c     ******************
c
c
         call init
c
         call grid
c
         call vtrans
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
      write(*,600) time
  600 format(' time=',f10.4)
c
         call poiss
c
         call boundp
c
         call convec
c
         call visc
c
         call euler
c
         call boundv
c
         call judge(errs)
c
        if(errs.le.1.d-6) goto 10
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/ omega1,eps1,iter1
      common /tmcnd/ dt,nstep
      common /reyno/ re,rei
      common /cont/ icon
c
c     ***** iteration parameter *****
c
         iter1=100
c
c     **** convergence parameter ****
c
         eps1=1.0d-3
c
c     **** acceleration parameter ****
c
         omega1=1.8d0
c
c     *** to be continue yes=1 no=0 ***
c
         icon=0
c
c     ****** time step *****
c
         dt=0.001d0
c
c     ***** reynolds number *****
c
         re =400.d0
         rei=1.d0/re
c
c     ***** time iteration *****
c
         nstep=5000
c
c
      return
      end
c
c     **************************
c     ** subroutine grid data **
c     **************************
c
      subroutine grid
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
      common /oldvel/ vo(mx,my,mz,3)
      common /coord/ x(mx),y(my),z(mz)
      common /cont/ icon
      common /past/ time
c
c     *** grid generation ***
c
      nx=31
      ny=nx
      nz=nx
c
      alpha=0.85d0
      at=dlog((1.d0+alpha)/(1.d0-alpha))/2.d0
c
      do i=1,nx
c
        eta=at*(-1.d0+2.d0*dble(i-1)/dble(nx-1))
        x(i)=(dtanh(eta)/alpha+1.d0)-1.d0
c
      end do
c
      do j=1,ny
        y(j)=x(j)
      end do
c
      do k=1,nz
        z(k)=x(k)
      end do
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
         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='rd2.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.1)
  520 format(3f10.5,e15.7)
c
      end if
c
c
      return
      end
c
c     ****************************
c     ** output result to file  **
c     ****************************
c
      subroutine output
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
      common /tmcnd/ dt,nstep
      common /coord/ x(mx),y(my),z(mz)
      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
c
      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
c
  600 format(3i10,f10.5,f10.1)
  610 format(3f10.5,e15.7)
c
      close(3)
c
c
      return
      end
c
c     *******************************
c     ** judgement of steady state **
c     *******************************
c
      subroutine judge(errs)
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
      common /oldvel/ vo(mx,my,mz,3)
c
c     **** judge for steady state ****
c
         errs=0.d0
         errv=0.d0
c
c     *** maximum error ***
c
      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,1)).gt.errv) then 
           errv=dabs(v(i,j,k,1))
           idi=i
           idj=j
           idk=k
           idv=1
         end if
         if(dabs(v(i,j,k,2)).gt.errv) then
           errv=dabs(v(i,j,k,2))
           idi=i
           idj=j
           idk=k
           idv=2
         end if
         if(dabs(v(i,j,k,3)).gt.errv) then
           errv=dabs(v(i,j,k,3))
           idi=i
           idj=j
           idk=k
           idv=3
         end if
c
  100 continue
c
c
       if(errv.gt.3.d0) then
         write(*,*)'***************************'
         write(*,*)'*** faulty terminated ! ***'
         write(*,*)'***************************'
         goto 99
       end if
c
         errs=dsqrt(errs/dfloat((nx-2)*(ny-2)*(nz-2)))
c
      write(*,*)'     '
      write(*,*)'**********************************************'
      write(*,660) errs
      write(*,670) idv,errv
      write(*,680) idi,idj,idk
      write(*,*)'**********************************************'
      write(*,*)'     '
c
  660 format('***  steady state  check=',e10.4)
  670 format('***  max. velocity check  ','v',i1,'=',f8.4)
  680 format('***  max. velocity point=',3i5)
c
c
      return
c
  99   call output
c
c
       stop
       end
c
c     ************************************
c     ** independent variable translate **
c     ************************************
c
      subroutine vtrans
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /coord/ x(mx),y(my),z(mz)
      common /coeff1/ x1(mx),x2(mx),e1(my),e2(my),z1(mz),z2(mz)
      common /metric/ xi(mx),et(my),ze(mz)
c
c
      do 100 i=2,nxm1
c
c     ***** 1st & 2nd order derivative *****
c
        xx   =.5d0*(x(i+1)-x(i-1))
        xi(i)=1.d0/xx
c
        xxx=x(i+1)-2.d0*x(i)+x(i-1)
c
c     *****  metric tensor  *****
c
        x1(i) = 1.d0/(xx**2)
        x2(i) =-xxx/(xx**3)
c
  100 continue
c
      do 110 j=2,nym1
c
c     ***** 1st & 2nd order derivative *****
c
        ye   =.5d0*(y(j+1)-y(j-1))
        et(j)=1.d0/ye
c
        yee  =y(j+1)-2.d0*y(j)+y(j-1)
c
c     *****  metric tensor  *****
c
        e1(j) = 1.d0/(ye**2)
        e2(j) =-yee/(ye**3)
c
  110 continue
c
c
      do 120 k=2,nzm1
c
c     ***** 1st & 2nd order derivative *****
c
        zz    =.5d0*(z(k+1)-z(k-1))
        ze(k) =1.d0/zz
c
        zzz=z(k+1)-2.d0*z(k)+z(k-1)
c
c     *****  metric tensor  *****
c
        z1(k) = 1.d0/(zz**2)
        z2(k) =-zzz/(zz**3)
c
  120 continue
c
c
      return
      end
c
c     *****************************************
c     ** setting velocity boundary condition **
c     *****************************************
c
      subroutine boundv
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
c
c     **** no-slip condition ****
c
      do 100 j=1,ny
      do 100 i=1,nx
c
         v(i,j,nzp1,1)= v(i,j,nzm1,1)
         v(i,j,nzp1,2)=-v(i,j,nzm1,2)
         v(i,j,nzp1,3)=-v(i,j,nzm1,3)
c
         v(i,j,nz  ,1)=1.d0
         v(i,j,nz  ,2)=0.d0
         v(i,j,nz  ,3)=0.d0
c
         v(i,j,0,1)=-v(i,j,2,1)
         v(i,j,0,2)=-v(i,j,2,2)
         v(i,j,0,3)=-v(i,j,2,3)
c
         v(i,j,1,1)=0.d0
         v(i,j,1,2)=0.d0
         v(i,j,1,3)=0.d0
c
  100 continue
c
c      **** wall condition ****
c
      do 120 k=2,nzm1
      do 120 i=2,nxm1
c
         v(i,0,k,1)=-v(i,2,k,1)
         v(i,0,k,2)=-v(i,2,k,2)
         v(i,0,k,3)=-v(i,2,k,3)
c
         v(i,1,k,1)=0.d0
         v(i,1,k,2)=0.d0
         v(i,1,k,3)=0.d0
c
         v(i,ny,k,1)=0.d0
         v(i,ny,k,2)=0.d0
         v(i,ny,k,3)=0.d0
c
         v(i,nyp1,k,1)=-v(i,nym1,k,1)
         v(i,nyp1,k,2)=-v(i,nym1,k,2)
         v(i,nyp1,k,3)=-v(i,nym1,k,3)
c
  120 continue
c
c     **** side wall condition ****
c
      do 130 k=1,nz
      do 130 j=1,ny
c
         v(0,j,k,1)=-v(2,j,k,1)
         v(0,j,k,2)=-v(2,j,k,2)
         v(0,j,k,3)=-v(2,j,k,3)
c
         v(1,j,k,1)=0.d0
         v(1,j,k,2)=0.d0
         v(1,j,k,3)=0.d0
c
         v(nx,j,k,1)=0.d0
         v(nx,j,k,2)=0.d0
         v(nx,j,k,3)=0.d0
c
         v(nxp1,j,k,1)=-v(nxm1,j,k,1)
         v(nxp1,j,k,2)=-v(nxm1,j,k,2)
         v(nxp1,j,k,3)=-v(nxm1,j,k,3)
c
  130 continue
c
c
      return
      end
c
c     *****************************************
c     ** setting pressure boundary condition **
c     *****************************************
c
      subroutine boundp
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
c
c     **** on the wall condition ****
c
         di3=1.d0/3.d0
c
      do 100 j=1,ny
      do 100 i=1,nx
c
          p(i,j,1 )=(4.d0*p(i,j,2   )-p(i,j,3   ))*di3
          p(i,j,nz)=(4.d0*p(i,j,nzm1)-p(i,j,nz-2))*di3
c
  100 continue
c
c     **** on the wall condition ****
c
      do 110 k=1,nz
      do 110 i=1,nx
c
          p(i,1 ,k)=(4.d0*p(i,2   ,k)-p(i,3   ,k))*di3
          p(i,ny,k)=(4.d0*p(i,nym1,k)-p(i,ny-2,k))*di3
c
  110 continue
c
c     **** on the wall condition ****
c
      do 120 k=1,nz
      do 120 j=1,ny
c
          p(1 ,j,k)=(4.d0*p(2   ,j,k)-p(3   ,j,k))*di3
          p(nx,j,k)=(4.d0*p(nxm1,j,k)-p(nx-2,j,k))*di3
c
  120 continue
c
c
      return
      end
c
c     *****************************
c     ** compute convective term **
c     *****************************
c
      subroutine convec
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
      common /metric/ xi(mx),et(my),ze(mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
c
c     *** compute convective & viscous term ***
c
      do 110 k=2,nzm1
c
         kp2=k+2
         kp1=k+1
         km1=k-1
         km2=k-2
c
      do 110 j=2,nym1
c
         jp2=j+2
         jp1=j+1
         jm1=j-1
         jm2=j-2
c
      do 110 i=2,nxm1
c
         ip2=i+2
         ip1=i+1
         im1=i-1
         im2=i-2
c
c     **** contravariant velocity component ****
c
      uu=v(i,j,k,1)*xi(i)
      vv=v(i,j,k,2)*et(j)
      ww=v(i,j,k,3)*ze(k)
c
c     *****  convective term  *****
c
      u1=.5d0*(v(ip1,j,k,1)-v(im1,j,k,1))
      u2=0.d0
      u3=.5d0*(v(i,jp1,k,1)-v(i,jm1,k,1))
      u4=0.d0
      u5=.5d0*(v(i,j,kp1,1)-v(i,j,km1,1))
      u6=0.d0
c
      v1=.5d0*(v(ip1,j,k,2)-v(im1,j,k,2))
      v2=0.d0
      v3=.5d0*(v(i,jp1,k,2)-v(i,jm1,k,2))
      v4=0.d0
      v5=.5d0*(v(i,j,kp1,2)-v(i,j,km1,2))
      v6=0.d0
c
      w1=.5d0*(v(ip1,j,k,3)-v(im1,j,k,3))
      w2=0.d0
      w3=.5d0*(v(i,jp1,k,3)-v(i,jm1,k,3))
      w4=0.d0
      w5=.5d0*(v(i,j,kp1,3)-v(i,j,km1,3))
      w6=0.d0
c
c
      conv(i,j,k,1)=(uu*u1+dabs(uu)*u2
     &              +vv*u3+dabs(vv)*u4
     &              +ww*u5+dabs(ww)*u6)
      conv(i,j,k,2)=(uu*v1+dabs(uu)*v2
     &              +vv*v3+dabs(vv)*v4
     &              +ww*v5+dabs(ww)*v6)
      conv(i,j,k,3)=(uu*w1+dabs(uu)*w2
     &              +vv*w3+dabs(vv)*w4
     &              +ww*w5+dabs(ww)*w6)
c
  110 continue
c
c
      return
      end
c
c     **************************
c     ** compute viscous term **
c     **************************
c
      subroutine visc
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,mz)
      common /coeff1/ x1(mx),x2(mx),e1(my),e2(my),z1(mz),z2(mz)
      common /metric/ xi(mx),et(my),ze(mz)
      common /navier/ conv(mx,my,mz,3),vis(mx,my,mz,3)
c
c     *** compute convective & viscous term ***
c
      do 110 k=2,nzm1
c
         kp2=k+2
         kp1=k+1
         km1=k-1
         km2=k-2
c
      do 110 j=2,nym1
c
         jp2=j+2
         jp1=j+1
         jm1=j-1
         jm2=j-2
c
      do 110 i=2,nxm1
c
         ip2=i+2
         ip1=i+1
         im1=i-1
         im2=i-2
c
c     *****  viscous term  *****
c
      vx1=(v(ip1,j,k,1)-2.d0*v(i,j,k,1)+v(im1,j,k,1))*x1(i)
     &   +(v(i,jp1,k,1)-2.d0*v(i,j,k,1)+v(i,jm1,k,1))*e1(j)
     &   +(v(i,j,kp1,1)-2.d0*v(i,j,k,1)+v(i,j,km1,1))*z1(k)
c
      vx2=.5d0*(v(ip1,j,k,1)-v(im1,j,k,1))*x2(i)
     &   +.5d0*(v(i,jp1,k,1)-v(i,jm1,k,1))*e2(j)
     &   +.5d0*(v(i,j,kp1,1)-v(i,j,km1,1))*z2(k)
c
      vy1=(v(ip1,j,k,2)-2.d0*v(i,j,k,2)+v(im1,j,k,2))*x1(i)
     &   +(v(i,jp1,k,2)-2.d0*v(i,j,k,2)+v(i,jm1,k,2))*e1(j)
     &   +(v(i,j,kp1,2)-2.d0*v(i,j,k,2)+v(i,j,km1,2))*z1(k)
c
      vy2=.5d0*(v(ip1,j,k,2)-v(im1,j,k,2))*x2(i)
     &   +.5d0*(v(i,jp1,k,2)-v(i,jm1,k,2))*e2(j)
     &   +.5d0*(v(i,j,kp1,2)-v(i,j,km1,2))*z2(k)
c
      vz1=(v(ip1,j,k,3)-2.d0*v(i,j,k,3)+v(im1,j,k,3))*x1(i)
     &   +(v(i,jp1,k,3)-2.d0*v(i,j,k,3)+v(i,jm1,k,3))*e1(j)
     &   +(v(i,j,kp1,3)-2.d0*v(i,j,k,3)+v(i,j,km1,3))*z1(k)
c
      vz2=.5d0*(v(ip1,j,k,3)-v(im1,j,k,3))*x2(i)
     &   +.5d0*(v(i,jp1,k,3)-v(i,jm1,k,3))*e2(j)
     &   +.5d0*(v(i,j,kp1,3)-v(i,j,km1,3))*z2(k)
c
        vis(i,j,k,1)=(vx1+vx2)*rei
        vis(i,j,k,2)=(vy1+vy2)*rei
        vis(i,j,k,3)=(vz1+vz2)*rei
c
  110 continue
c
c
      return
      end
c
c     **********************************
c     ** poisson equation of pressure **
c     **********************************
c
      subroutine poiss
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      common /divide/ nxp1,nyp1,nzp1,nx,ny,nz,nxm1,nym1,nzm1
      common /conv1/ omega1,eps1,iter1
      common /tmcnd/ dt,nstep
      common /reyno/ re,rei
      common /result/ v(0:mx+1,0:my+1,0:mz+1,3),p(mx,my,mz)
      common /coeff1/ x1(mx),x2(mx),e1(my),e2(my),z1(mz),z2(mz)
      common /metric/ xi(mx),et(my),ze(mz)
      dimension rhs(mx,my,mz)
c
c     *** right hand side of poisson eqation ***
c
c
      do 110 k=2,nzm1
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nxm1
c
         ip1=i+1
         im1=i-1
c
      ux=.5d0*xi(i)*(v(ip1,j,k,1)-v(im1,j,k,1))
      vx=.5d0*xi(i)*(v(ip1,j,k,2)-v(im1,j,k,2))
      wx=.5d0*xi(i)*(v(ip1,j,k,3)-v(im1,j,k,3))
c
      uy=.5d0*et(j)*(v(i,jp1,k,1)-v(i,jm1,k,1))
      vy=.5d0*et(j)*(v(i,jp1,k,2)-v(i,jm1,k,2))
      wy=.5d0*et(j)*(v(i,jp1,k,3)-v(i,jm1,k,3))
c
      uz=.5d0*ze(k)*(v(i,j,kp1,1)-v(i,j,km1,1))
      vz=.5d0*ze(k)*(v(i,j,kp1,2)-v(i,j,km1,2))
      wz=.5d0*ze(k)*(v(i,j,kp1,3)-v(i,j,km1,3))
c
        div=ux+vy+wz
c
        divcon=ux*ux+vy*vy+wz*wz+2.d0*(vx*uy+uz*wx+vz*wy)
c
        rhs(i,j,k)=divcon-div/dt
c
  110 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,nzm1
c
         kp1=k+1
         km1=k-1
c
      do 130 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 130 i=2,nxm1
c
         ip1=i+1
         im1=i-1
c
        po=p(i,j,k)
c
      pp=x1(i)*(p(ip1,j,k)+p(im1,j,k))
     &  +x2(i)*(p(ip1,j,k)-p(im1,j,k))*.5d0
     &  +e1(j)*(p(i,jp1,k)+p(i,jm1,k))
     &  +e2(j)*(p(i,jp1,k)-p(i,jm1,k))*.5d0
     &  +z1(k)*(p(i,j,kp1)+p(i,j,km1))
     &  +z2(k)*(p(i,j,kp1)-p(i,j,km1))*.5d0+rhs(i,j,k)
c
c     **** cal. residual ****
c
        dp=.5d0*pp/(x1(i)+e1(j)+z1(k))-p(i,j,k)
c
        p(i,j,k)=po+omega1*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
      write(*,600) loop,err
  600 format(' poisson iteration & error=',i5,2x,e8.3)
c
c
      return
      end
c
c     **********************
c     ** time integration **
c     **********************
c
      subroutine euler
c
      parameter(mx=51,my=51,mz=51)
      implicit double precision(a-h,o-z)
      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(mx,my,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     *** time integration by euler forward method ***
c
c     *** velocity interchange ***
c
      do 100 k=2,nzm1
      do 100 j=2,nym1
      do 100 i=2,nxm1
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
      do 110 k=2,nzm1
c
         kp1=k+1
         km1=k-1
c
      do 110 j=2,nym1
c
         jp1=j+1
         jm1=j-1
c
      do 110 i=2,nxm1
c
         ip1=i+1
         im1=i-1
c
c     *****  pressure term  *****
c
      px=xi(i)*(p(ip1,j,k)-p(im1,j,k))*.5d0
      py=et(j)*(p(i,jp1,k)-p(i,jm1,k))*.5d0
      pz=ze(k)*(p(i,j,kp1)-p(i,j,km1))*.5d0
c
c     *** new velocity ***
c
      v(i,j,k,1)=vo(i,j,k,1)+dt*(-conv(i,j,k,1)-px+vis(i,j,k,1))
      v(i,j,k,2)=vo(i,j,k,2)+dt*(-conv(i,j,k,2)-py+vis(i,j,k,2))
      v(i,j,k,3)=vo(i,j,k,3)+dt*(-conv(i,j,k,3)-pz+vis(i,j,k,3))
c
  110 continue
c
c
      return
      end