C     parameter list:
C     imax : number of gridpoints in the x-direction
C     jmax : number of gridpoints in the y-direction
C     kmax : number of gridpoints in the z-direction
C     lx   : length of the computational domain in the x-dir.
C     ly   : length of the computational domain in the y-dir.
C     lz   :    .....................................
C     Cs   : Smagorinsky constant
C     dxi=1/dx  := imax / lx         
C     dyi=1/dy  := jmax / ly
C     dzi=1/dz  := kmax / lx
C     Rei   : 1 / Re  where Re is the Reynolds number
C     nstap : Number of timesteps
C     Uold(0:imax+1,0:jmax+1,0:kmax+1) velocity in
C                  x-direction at the oldest timestep
C     Vold(0:imax+1,0:jmax+1,0:kmax+1) velocity in
C                  y-direction at the oldest timestep
C     Wold(0:imax+1,0:jmax+1,0:kmax+1) velocity in
C                  z-direction at the oldest timestep
C     Unew(.....)  velocity at the present timestep in x-direction
C     Vnew(.....)  velocity at the present timestep in y-direction
C     Wnew(.....)  velocity at the present timestep in z-direc.
C
C     dUdt(.....),dVdt(...),dWdt(....) predicted velocity at
C     the new timestep
C
C     P(imax,jmax,kmax) : Pressure array
C     visc(0:imax+1,0:jmax+1,0:kmax+1)  : viscosity array used
C                                         in LES.
C*****************************************************************
C*****************************************************************
C*****************************************************************
C*****************************************************************
C  Subroutines ::::
C  subroutine init    : set U,V,W to some initial value
C  subroutine chkdt   : set timestep   dt
C  subrouitne chkdiv  : checkes the divergence
C  subroutine momx    : integrates the momentum equation in x-dir
C  subroutine momy    : integrates the momentum equation in y-dir
C  subroutine momz    : integrates the momentum equation in z-dir
C  subroutine adams   : integrates the momentum equation in time
C  subroutine fillps  : fills the right hand side for the Poisson solver
C  subroutine correc  : corrects the velocity such that div(u,v,w)==0
C  subroutine bound   : set boundary conditions
C  subroutine solver  : Calls a Fishpak routine to solve Poisson equation 
C  subroutine submod  : subgrid model      
C**************************************************************************
C
C     Describtion of the model.
C 
C     In this model the Navier-Stokes equations are integrated
C     in space with a Finite Volume Technique on a staggerd grid.
C
C                           W          V
C                           ^         ^ 
C                           |        /
C                       ____|___________
C                     /|    |    /      /
C                   /  |    |  /      / |
C                 /    |     /      /   |
C               /      |          /     |
C              ------------*-----   -------> U  
C              |       |________|______|
C          dz  |     /          |      /  
C              |   /            |    / 
C              | /              |  / dy
C              |                |/
C              ------------------
C                     dx
C
C     U : velocity in x-direction    [ U(i,j,k) ]
C     V : velocity in y-direction    [ V(i,j,k) ]
C     W : velocity in z-direction    [ W(i,j,k) ]
C     * : Pressure point
C     dx,dy,dz  dimensions of the grid cell
C
C     For instance:
C
C     d u u    (U(i,j,k)+U(i+1,j,k))**2 - (U(i,j,k)+U(i-1,j,k))**2 
C     ----- =  ---------------------------------------------------
C      d x                         4 dx
C
C       etc.
C********************************************************************
C
C
C        ADV  = ADVECTION
C        DIFF = DIFFUSION
C        PRES = PRESSURE (GRADIENT)
C
C
C      Time integration with Adams-Bashfort / pressure correction
C
C
C       dUdt = 1.5(-ADV+DIFF)_new -0.5(-ADV+DIFF)_old
C       div grad (P) =div (dudt,dvdt,dwdt)
C       Unew = dUdt - dt grad(P)
C********************************************************************
      implicit none
      include 'param.txt'
      include 'common.txt'
      real    lx,ly,lz,Re
      integer istap,nstap
      real    time,utime
      read(5,*) Ufree
      read(5,*) lx
      read(5,*) ly
      read(5,*) lz
      read(5,*) Re
      read(5,*) nstap
********
      dxi = imax / lx
      dyi = jmax / ly
      dzi = kmax / lz
      Rei = 1.   / Re
      do k=0,k1
	do j=0,j1
	  do i=0,i1
	    visc(i,j,k)=Rei
          enddo
         enddo
      enddo
      call init
c      call loadd(0)
      call bound(Uold,Vold,Wold)
      call bound(Unew,Vnew,Wnew)
      call chkdt
*******
      time  =0.
      utime =0.
C  main loop starts here

      do istap = 1,nstap
      time = time + dt
      utime=utime + dt
      write(6,*) time,dt
      call adamsb('adams')
      call bound   (Unew,Vnew,Wnew)
      call bound   (dUdt,dVdt,dWdt)
      call fillps  (istap)    
      call solver  (p,dxi,dyi,dzi)
      call correc
      call bound   (Uold,Vold,Wold)
      call bound   (Unew,Vnew,Wnew)
      if (mod(istap,5).eq.0) then
      call chkdt
      call chkdiv
      endif
      if ( utime .gt. 0.01) then
      utime = 0.
      call avs 
      write(6,*) ' Output generated '
      endif 
      enddo
      stop
      end

      subroutine adamsb(string)
C
C
C  Integrate in time with second order Adams-Bashfort (or Euler forward)
C  dold(...),dnew(....) are dummy arrays
C
      implicit none
      include 'param.txt'
      include 'common.txt'
      real   dold(0:i1,0:j1,0:k1),dnew(0:i1,0:j1,0:k1)
      real   cof1,cof2
      character*5 string
*     Adams-Bashfort  cof1=1.5  , cof2 = -0.5
*     Euler-Forward   cof1=  1  , cof2 =    0
      if (string .eq. 'euler') then
      cof1=1.
      cof2=0.
      endif
      if (string .eq. 'adams') then
      cof1= 1.5
      cof2=-0.5
      endif
c
c   integrate momentum equations at old and new timelevel
c
      call momx(dold,uold,vold,wold,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)
      call momx(dnew,unew,vnew,wnew,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)
c    
c   calculate predicted velocity
c
      do k=1,kmax
        do j=1,jmax
          do i=1,imax
          dudt(i,j,k)=Unew(i,j,k) +
     1    dt * ( cof1 * dnew(i,j,k) + cof2 * dold(i,j,k) )
          enddo
        enddo
      enddo
      call momy(dold,uold,vold,wold,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)
      call momy(dnew,unew,vnew,wnew,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)

      do k=1,kmax
        do j=1,jmax
          do i=1,imax
          dvdt(i,j,k)=Vnew(i,j,k)+
     1    dt * ( cof1 * dnew(i,j,k) + cof2 * dold(i,j,k) )
          enddo
        enddo
      enddo
      call momz(dold,uold,vold,wold,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)
      call momz(dnew,unew,vnew,wnew,dxi,dyi,dzi,
     &          1,imax,1,jmax,1,kmax,i1,j1,k1,visc)
      do k=1,kmax
        do j=1,jmax
          do i=1,imax
          dwdt(i,j,k)=Wnew(i,j,k)+
     1    dt * ( cof1 * dnew(i,j,k) + cof2 * dold(i,j,k) )
          enddo
        enddo
      enddo
      return
      end

      subroutine fillps
      implicit none
      include 'param.txt'
      include 'common.txt'
      real   dti
      dti =1./dt
c
c
c     fill right hand side of the poisson equation
c
c     discrete divergence is
c   w(i,j,k)-w(i,j,k-1)  v(i,j,k)-v(i,j-1,k)  u(i,j,k)-u(i-1,j,k)
c   -------------------  -------------------  -------------------  = div
c          d z                  d y                  d x 
c
c  note that p is not yet the pressure
c  but only the rhs of the poisson equation
c
      do k=1,kmax
       do j=1,jmax
	  do i=1,imax
      p(i,j,k)= dti*(
     1   ( dWdt(i,j,k)-dWdt(i,j,k-1))*dzi+
     2   ( dVdt(i,j,k)-dVdt(i,j-1,k))*dyi+
     3   ( dUdt(i,j,k)-dUdt(i-1,j,k))*dxi
     e              )
          enddo
       enddo
      enddo
      return
      end


      subroutine correc
c
c
c   Corrects the velocity in such a way that div(u,v,w) is exactly zero
c
c   
      implicit none
      include 'param.txt'
      include 'common.txt'
      do k=1,kmax
      do j=1,jmax
      do i=1,imax-1
       dudt(i,j,k)= 
     1  dUdt(i,j,k) - dt * dxi * ( P(i+1,j,k)-P(i,j,k) )
      enddo
      enddo
      enddo
      do k=1,kmax
      do j=1,jmax
      dudt(imax,j,k)= 
     1  dudt(imax,j,k) -dt*  dxi *(P(1,j,k)-P(imax,j,k))
      enddo
      enddo 

      do k=1,kmax
      do j=1,jmax-1
      do i=1,imax
       dvdt(i,j,k)= 
     1  dvdt(i,j,k) - dt*dyi * ( P(i,j+1,k)-P(i,j,k) )
      enddo
      enddo
      enddo
      do k=1,kmax
      do i=1,imax
      dvdt(i,jmax,k)=
     1  dvdt(i,jmax,k)-dt*dyi*(P(i,1,k)-P(i,jmax,k))
      enddo
      enddo
      do k=1,kmax-1
      do j=1,jmax
      do i=1,imax
       dwdt(i,j,k)=
     1 dwdt(i,j,k) - dt*dzi * ( P(i,j,k+1)-P(i,j,k) )
      enddo
      enddo
      enddo
c
c
c    Update the velocity fields
c
c
      do k=1,kmax
         do j=1,jmax
            do i=1,imax
            Uold(i,j,k)=Unew(i,j,k)
            Unew(i,j,k)=dUdt(i,j,k)
            Vold(i,j,k)=Vnew(i,j,k)
            Vnew(i,j,k)=dVdt(i,j,k)
            Wold(i,j,k)=Wnew(i,j,k)
            Wnew(i,j,k)=dWdt(i,j,k)
            enddo
        enddo
      enddo
      return
      end

      subroutine avs
c
c     output routine for AVS
c
      include 'param.txt'
      include 'common.txt'
      real vor(imax,kmax)
      open(30,file='tecdata')
      write(30,*) 'VARIABLES ="X", "Y", "U-vel","V-vel","W-vel","P","VORT"'
      write(30,*) ' ZONE I =       ',imax,'  J =         ',kmax ,'   F = POINT'

      do k=1,kmax
	 do i=1,imax
	 write(30,'(10E18.7)') i/dxi,k/dzi,
     ^	 Unew(i,j,k),Vnew(i,j,k),Wnew(i,j,k), p(i,j,k),
     ^ (Unew(i,j,k+1)-Unew(i,j,k))*dzi  - (Wnew(i+1,j,k)-Wnew(i,j,k))*dxi
        enddo
      enddo
      close(30)

c
c output for Gnuplot
c 
      open (24,file='congnu')
      do k=1,kmax
	do i=1,imax
	  write(24,111) 1.*i/dxi,1.*k/dzi,vor(i,k)
        enddo
      enddo
111   FORMAT(3F12.6)
      close(24)
      Return
      end
c
c
      subroutine chkdiv
c
c checks the discrete divergence (should be of order of machine prec.)  
c
      implicit none
      include 'param.txt'
      include 'common.txt'
      real div,divmax,divtot
      divmax = -9999.9999
      divtot = 0.
      do k=1,kmax
         do j=1,jmax
            do i=1,imax
            div = (
     1   ( Wnew(i,j,k)-Wnew(i,j,k-1) ) * dzi +
     2   ( Vnew(i,j,k)-Vnew(i,j-1,k) ) * dyi +
     3   ( Unew(i,j,k)-Unew(i-1,j,k) ) * dxi )
           divtot = divtot + div
           divmax = max ( div,divmax)
c      write(6,*) i,j,k,div
           enddo
         enddo
      enddo
      write(6,111) divtot,divmax
111   FORMAT('Mass loss/gain : Tot =',e13.6,' Max = ',e13.6)
      if (divtot .gt. 10e-3) stop 'Fatal Error **DIVERGENCE TO LARGE**'
      return
      end

      subroutine chkdt
c
c     Calculates the timestep dt.
c
      implicit none
      include 'param.txt'
      include 'common.txt'
      real    Courant,dtemp,dtmax
      Courant  = .5
      dtmax = -9999.9999
      do k=1,kmax
         do j=1,jmax
            do i=1,imax
            dtemp =
     1      abs(Unew(i,j,k) * dxi) +
     2      abs(Vnew(i,j,k) * dyi) +
     3      abs(Wnew(i,j,k) * dzi) +
     4       visc(i,j,k) * (dxi*dxi + dyi*dyi + dzi*dzi)
            dtmax = max ( dtmax , dtemp )
            enddo
         enddo
      enddo
      dt = Courant / dtmax
      return
      end

      subroutine bound(U,V,W)
      include 'param.txt'
      include 'common.txt'
      real U(0:i1,0:j1,0:k1),V(0:i1,0:j1,0:k1),
     $     W(0:i1,0:j1,0:k1)
c   
c   boundary conditions at z=zmin and z =zmax
c
      do i=0,i1
        do j=0,j1
         W(i,j,0    ) =      0. 
         W(i,j,kmax ) =      0.
         V(i,j,0    ) =  V(i,j,1   )
         V(i,j,k1   ) =  V(i,j,kmax)
         U(i,j,0    ) =  U(i,j,1   )
         U(i,j,k1   ) =  U(i,j,kmax)
        enddo
      enddo
c
c  periodic boundary conditions in the x direction
c
      do k=0,k1
        do j=0,j1
          U(0 ,j,k)=U(imax,j,k)
          V(0 ,j,k)=V(imax,j,k)
          W(0 ,j,k)=W(imax,j,k)
          U(i1,j,k)=U(1   ,j,k)
          V(i1,j,k)=V(1   ,j,k)
          W(i1,j,k)=W(1   ,j,k)
         enddo
       enddo
c
c periodic boundary conditions in the y-direction
c
       do k=0,k1
         do i=0,i1
           U(i,0,k )=U(i,jmax,k)
           V(i,0,k )=V(i,jmax,k)
           W(i,0,k )=W(i,jmax,k)
           U(i,j1,k)=U(i,1   ,k)
           V(i,j1,k)=V(i,1   ,k)
           W(i,j1,k)=W(i,1   ,k)
         enddo
       enddo
       return
       end


      subroutine momx(dudt,u,v,w,dxi,dyi,dzi,ib,ie,jb,je,kb,ke,i1,j1,k1,
     $                visc)
      implicit none
      integer ib,ie,jb,je,kb,ke,i1,j1,k1
      integer ip,im,jp,jm,kp,km,i,j,k
      real dudt(0:i1,0:j1,0:k1),
     $     u   (0:i1,0:j1,0:k1),
     $     v   (0:i1,0:j1,0:k1),
     $     w   (0:i1,0:j1,0:k1),
     $    dxi,dyi,dzi,visc(0:i1,0:j1,0:k1)
      real uuip ,uuim ,uvjp ,uvjm ,uwkp ,uwkm,vzp,vzm,vyp,vym
      real du1ip,du1im,du1jp,du1jm,du1kp,du1km
      do k=kb,ke
        do j=jb,je
          do i=ib,ie
*
          ip = i + 1
          jp = j + 1
          kp = k + 1
          im = i - 1
          jm = j - 1
          km = k - 1
          uuip  = 0.25 * ( U(ip,j,k)+U(i,j,k) )*( U(ip,j,k)+U(i,j,k)  )
          uuim  = 0.25 * ( U(im,j,k)+U(i,j,k) )*( U(im,j,k)+U(i,j,k)  )
          uvjp  = 0.25 * ( U(i,jp,k)+U(i,j,k) )*( V(ip,j,k)+V(i,j,k)  )
          uvjm  = 0.25 * ( U(i,jm,k)+U(i,j,k) )*( V(ip,jm,k)+V(i,jm,k))
          uwkp  = 0.25 * ( U(i,j,kp)+U(i,j,k) )*( W(ip,j,k) +W(i,j,k) )
          uwkm  = 0.25 * ( U(i,j,km)+U(i,j,k) )*( W(ip,j,km)+W(i,j,km))
          du1ip = 2.*( dxi * (U(ip,j,k)-U(i ,j,k) ))
          du1im = 2.*( dxi * (U(i ,j,k)-U(im,j,k) ))
          du1jp = dyi*(U(i,jp,k)-U(i, j,k))+dxi*(V(ip,j, k)-V(i ,j ,k))
          du1jm = dyi*(U(i,j ,k)-U(i,jm,k))+dxi*(V(ip,jm,k)-V(i ,jm,k))
          du1kp = dzi*(U(i,j,kp)-U(i ,j,k))+dxi*(W(ip,j, k)-W(i ,j ,k))
          du1km = dzi*(U(i,j ,k)-U(i,j,km))+dxi*(W(ip,j,km)-W(i ,j,km))
      vzp =
     &   0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,j,kp)+visc(i,j,kp))
      vzm =
     &   0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,j,km)+visc(i,j,km))
      vyp =
     &   0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,jp,k)+visc(i,jp,k))
      vym =
     &   0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,jm,k)+visc(i,jm,k))

          dudt(i,j,k) =
     1   dxi*( -uuip + visc(ip,j,k)*du1ip+uuim - visc(i,j,k)*du1im )+
     2   dyi*( -uvjp + vyp*du1jp+uvjm - vym*du1jm )+
     3   dzi*( -uwkp + vzp*du1kp+uwkm - vzm*du1km )
            enddo
         enddo
      enddo
      return
      end

      subroutinemomy(dvdt,u,v,w,dxi,dyi,dzi,ib,ie,jb,je,kb,ke,i1,j1,k1,
     $                visc)
      implicit none
      integer ib,ie,jb,je,kb,ke,i1,j1,k1
      integer im,ip,jm,jp,km,kp,i,j,k
      real dvdt(0:i1,0:j1,0:k1),
     $     u   (0:i1,0:j1,0:k1),
     $     v   (0:i1,0:j1,0:k1),
     $     w   (0:i1,0:j1,0:k1),
     $    dxi,dyi,dzi,visc(0:i1,0:j1,0:k1)
      real uvip ,uvim ,vvjp ,vvjm ,wvkp ,wvkm
      real dv1ip,dv1im,dv1jp,dv1jm,dv1kp,dv1km
      real vxp,vxm,vzp,vzm
      do k=kb,ke
        do j=jb,je
          do i=ib,ie
*
          ip = i + 1
          jp = j + 1
          kp = k + 1
          im = i - 1
          jm = j - 1
          km = k - 1
          uvip  = 0.25 * (U(i ,j,k)+U(i ,jp,k))*( V(i,j,k)+V(ip,j,k) )
          uvim  = 0.25 * (U(im,j,k)+U(im,jp,k))*( V(i,j,k)+V(im,j,k) )
          vvjp  = 0.25 * (V(i,j,k )+V(i,jp,k) )*( V(i,j,k)+V(i,jp,k) )
          vvjm  = 0.25 * (V(i,j,k )+V(i,jm,k) )*( V(i,j,k)+V(i,jm,k) )
          wvkp  = 0.25 * (W(i,j,k )+W(i,jp,k) )*( V(i,j,kp)+V(i,j,k) )
          wvkm  = 0.25 * (W(i,j,km)+W(i,jp,km))*( V(i,j,km)+V(i,j,k) )

          dv1ip = dxi*(V(ip,j,k)-V(i ,j,k))+dyi*(U(i,jp ,k)-U(i ,j,k))
          dv1im = dxi*(V(i ,j,k)-V(im,j,k))+dyi*(U(im,jp,k)-U(im,j,k))
          dv1jp = 2.*dyi*(V(i,jp,k)-V(i,j ,k))
          dv1jm = 2.*dyi*(V(i,j ,k)-V(i,jm,k))
          dv1kp = dzi*(V(i,j,kp)-V(i, j,k))+dyi*(W(i,jp,k )-W(i ,j,k))
          dv1km = dzi*(V(i,j,k )-V(i,j,km))+dyi*(W(i,jp,km)-W(i,j,km))
      vxp =
     & 0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,jp,k)+visc(i,jp,k))
      vxm =
     & 0.25*(visc(i,j,k)+visc(im,j,k)+visc(im,jp,k)+visc(i,jp,k))    
      vzp =
     & 0.25*(visc(i,j,k)+visc(i,jp,k)+visc(i,jp,kp)+visc(i,j,kp))
      vzm =
     & 0.25*(visc(i,j,k)+visc(i,jp,k)+visc(i,jp,km)+visc(i,j,km))

          dvdt(i,j,k)=
     1 dxi*((-uvip + vxp         *dv1ip)-(-uvim + vxm        *dv1im))+
     2 dyi*((-vvjp + visc(i,jp,k)*dv1jp)-(-vvjm + visc(i,j,k)*dv1jm))+
     3 dzi*((-wvkp + vzp         *dv1kp)-(-wvkm + vzm        *dv1km))
          enddo
        enddo
      enddo
      return
      end
      subroutinemomz(dwdt,u,v,w,dxi,dyi,dzi,ib,ie,jb,je,kb,ke,i1,j1,k1,
     $                visc)
      implicit none
      integer ib,ie,jb,je,kb,ke,i1,j1,k1
      integer im,jm,km,ip,jp,kp,i,j,k
      real dwdt(0:i1,0:j1,0:k1),
     $     u   (0:i1,0:j1,0:k1),
     $     v   (0:i1,0:j1,0:k1),
     $     w   (0:i1,0:j1,0:k1),
     $    dxi,dyi,dzi,visc(0:i1,0:j1,0:k1)
      real uwip ,uwim ,vwjp ,vwjm ,wwkp ,wwkm
      real dw1ip,dw1im,dw1jp,dw1jm,dw1kp,dw1km
      real  vxm,vym,vxp,vyp
      do k=kb,ke
        do j=jb,je
          do i=ib,ie
*
          ip = i + 1
          jp = j + 1
          kp = k + 1
          im = i - 1
          jm = j - 1
          km = k - 1
          uwip  = 0.25 * ( W(i,j,k)+W(ip,j,k))*(U(i ,j,k)+U(i ,j,kp) )
          uwim  = 0.25 * ( W(i,j,k)+W(im,j,k))*(U(im,j,k)+U(im,j,kp) )
          vwjp  = 0.25 * ( W(i,j,k)+W(i,jp,k))*(V(i ,j,k)+V(i,j ,kp) )
          vwjm  = 0.25 * ( W(i,j,k)+W(i,jm,k))*(V(i,jm,k)+V(i,jm,kp) )
          wwkp  = 0.25 * ( W(i,j,k)+W(i,j,kp))*(W(i ,j,k)+W(i,j,kp ) )
          wwkm  = 0.25 * ( W(i,j,k)+W(i,j,km))*(W(i ,j,k)+W(i,j,km ) )

          dw1ip = dxi*(W(ip,j,k)-W(i ,j,k))+dzi*(U(i ,j,kp)-U(i ,j,k))
          dw1im = dxi*(W(i ,j,k)-W(im,j,k))+dzi*(U(im,j,kp)-U(im,j,k))
          dw1jp = dyi*(W(i,jp,k)-W(i, j,k))+dzi*(V(i ,j,kp)-V(i ,j,k))
          dw1jm = dyi*(W(i,j ,k)-W(i,jm,k))+dzi*(V(i,jm,kp)-V(i,jm,k))
          dw1kp = 2*dzi*(W(i,j,kp)-W(i,j ,k))
          dw1km = 2*dzi*(W(i,j,k )-W(i,j,km))
          vxp =
     & 0.25*(visc(i,j,k)+visc(ip,j,k)+visc(ip,j,kp)+visc(i,j,kp))
          vxm =
     & 0.25*(visc(i,j,k)+visc(im,j,k)+visc(im,j,kp)+visc(i,j,kp))
          vyp =
     & 0.25*(visc(i,j,k)+visc(i,jp,k)+visc(i,jp,kp)+visc(i,j,kp))
          vym =
     & 0.25*(visc(i,j,k)+visc(i,jm,k)+visc(i,jm,kp)+visc(i,j,kp))
          dwdt(i,j,k) =
     1dxi*((-uwip + vxp         *dw1ip)-(-uwim + vxm        *dw1im))+
     2dyi*((-vwjp + vyp         *dw1jp)-(-vwjm + vym        *dw1jm))+
     3dzi*((-wwkp + visc(i,j,kp)*dw1kp)-(-wwkm + visc(i,j,k)*dw1km))
            enddo
         enddo
      enddo
      return
      end


      subroutine init
      include 'param.txt'
      include 'common.txt'
      real delta,z,x  
      delta = 1/28.
      do k=1,kmax
         do j=1,jmax
           do i=1,imax
           Uold(i,j,k)=0.
           Vold(i,j,k)=0.
           Wold(i,j,k)=0.
           Unew(i,j,k)=Ufree*tanh(28.* ((1.*k/kmax) -.5  ))
           Uold(i,j,k)=Unew(i,j,k)
           Vnew(i,j,k)=0.
           Wnew(i,j,k)=0.
           dUdt(i,j,k)=0.
           dVdt(i,j,k)=0.
           dWdt(i,j,k)=0.
           enddo
         enddo
      enddo
      do k=1,kmax
        z= (1.*k/kmax)-0.5
	do j=1,jmax
           do i=1,imax
           x = 2.*i/imax 
	     Unew(i,j,k)=Unew(i,j,k)+0.1*Ufree*(-z/delta)*
     &    exp( -(z/delta)**2)*cos(0.44/delta*x) 
             Wnew(i,j,k)=Wnew(i,j,k)+0.1*Ufree*      
     &    exp( -(z/delta)**2)*sin(0.44/delta*x)
           enddo
        enddo
      enddo

      return
      end
      subroutine loadd(in)
      implicit none
      include 'param.txt'
      include 'common.txt'
      integer reclengte,in
      reclengte=8
      if(in.eq.0) then
      open(15,file='start',access='direct',recl=6*(imax+2)*(jmax+2)*
     ^     reclengte)
      do 100 k=0,k1
      read(15,rec=(k+1))
     &  ( (Unew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Vnew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Wnew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Uold(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Vold(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Wold(i,j,k) , j=0,j1) , i=0,i1)
100   continue
      close(15)
      endif
      if(in.eq.1) then
      open(25,file='start',access='direct',recl=6*(imax+2)*(jmax+2)*
     ^      reclengte)
      do 200 k=0,k1
      write(25,rec=(k+1))
     &  ( (Unew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Vnew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Wnew(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Uold(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Vold(i,j,k) , j=0,j1) , i=0,i1),
     &  ( (Wold(i,j,k) , j=0,j1) , i=0,i1)
200   continue
      close(25)
      endif
      return
      end


      subroutine solver(p,dxi,dyi,dzi)
      implicit none
      include 'param.txt'
      integer  ier
      real     p(imax,jmax,kmax)
      real     dxi,dyi,dzi,dzi2,c1,c2,a(kmax),b(kmax),c(kmax)
      real     rhs(kmax,jmax,imax),d(kmax)
      real     save(k1*jmax*imax+3*kmax+3*jmax+4*imax+1000)
      real     w(i1*j1*k1)
      c1   = dxi*dxi
      c2   = dyi*dyi
      dzi2 = dzi*dzi
      do k=1,kmax
      a(k)=dzi2
      c(k)=dzi2
      b(k)=-(a(k)+c(k))
      d(k)=c2
      enddo
      b(1)=b(1)+a(1)
      a(1)=0.
      b(kmax)=b(kmax)+c(kmax)
      c(kmax)=0.
      callpois3d(0,imax,c1,0,jmax,c2,1,kmax,a,b,c,imax,jmax,p,ier,w)
      if (ier .ne. 0) write(6,*) ier
      If (ier .ne. 0) stop 'Fatal Error **Poisson Solver**'
      return
      end