      program anb
c
************************************************************************
*                                                                      *
*     anb                                                              *
*                                ... means: anb's not best             *
*                                           -     -   -                *
*     A lattice Boltzmann CFD teaching code.                           *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             Contact: research AT bernsdorf.org                       *
*                                                                      *
*             WWW: http://bernsdorf.org/research                       *
*                                                                      *
*     Corrections and hints from Thomas Zeiser and Wolfgang Hamm       *
*     are gratefully acknowledged.                                     *
*                                                                      *
*     Copyright (C) 1998-2008  Joerg BERNSDORF                         *
*                                                                      *
*     This program is free software; you can redistribute it and/or    *
*     modify it under the terms of the GNU General Public License as   *
*     published by the Free Software Foundation; either version 2 of   *
*     the License, or (at your option) any later version.              *
*                                                                      *
*     This program is distributed in the hope that it will be useful,  *
*     but WITHOUT ANY WARRANTY; without even the implied warranty of   *
*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    *
*     GNU General Public License for more details.                     *
*                                                                      *
*     You should have received a copy of the GNU General Public        *
*     License along with this program; if not, write to the Free       *
*     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139,    *
*     USA.                                                             *
*                                                                      *
*                                                                      *
*     Last change: 2008/06/09                                          *
*                                                                      *
************************************************************************
c
c.....short introduction
c
c.....anb simulates incompressible viscous flow as governed by the 
c     Navier-Stokes equations.
c     But anb is not a Navier-Stokes solver. 
c     anb is a lattice Boltzmann solver, which means, the velocity 
c     discrete boltzmann equation is solved here, using a single time
c     relaxation (BGK) collision operator.
c
c     For more details on this method read:
c             
c                         Y.H.Quian, D.D'Humieres and P.Lallemand,
c                         Lattice BGK Models for Navier-Stokes Equation
c                         Europhys. Lett., 17 (6), pp. 479-484 (1992)
c
c     Other methods belonging to the lattice gas / lattice Boltzmann
c     familiy exist, but this one is very simple, quite good and used
c     by lots of scientists working on this topic. 
c
c     This is a very simple implementation of the lattice BGK scheme,
c     not a very efficient and not a very memory safing one. 
c     *  Do not misunderstand this as a good proposal for writing an 
c        efficient lattice Boltzmann code ! 
c     *  Do not use this code for memory- or time-intensive 
c        computations, it is not even a research code !
c     *  Don't even think about commercial application ! 
c
c     This code was written to show beginners in a simple and 
c     short way the relevant procedures of a lattice Boltzmann solver,
c     pointing on how everything works "in principle". Nearly all 
c     procedures could be implemented other (and better) as it is done
c     here, and even the algorithms used here could be changed to 
c     save memory and increase performance. But the code works correct,
c     and we hope it will be good starting point for the first steps
c     in the lattice Boltzmann field. Good luck ! 
c
      implicit none
c     
c.....parameters
c
c.....grid size in x- and y-dimension
c
      integer  lx,ly
c
c.....ENTER APPROPRIATE VALUES HERE, TAKE CARE: SIZE DOES MATTER !!
c
      parameter(lx=200,ly=100)
c
c.....variables
c
c.....fluid density per link
c
      real*8  density
c
c.....relaxation parameter
c
      real*8  omega 
c
c.....accelleration 
c
      real*8  accel
c
c.....maximum number of iterations
c
      integer  t_max
c
c.....linear dimension for Reynolds number
c
      real*8 r_rey
c
c.....iteration counter
c
      integer time
c
c.....error flag
c
      logical  error
c
c.....obstacle array
c
      logical  obst(lx,ly)
c
c.....fluid densities
c     a 9-speed lattice is used here, other geometries are possible
c
c     the densities are numbered as follows:
c
c              6   2   5
c                \ | /
c              3 - 0 - 1
c                / | \
c              7   4   8
c
c    the lattice node are numbered as follows:
c
c
c     ^
c     |
c     y
c
c        :    :    :
c  
c   3    *    *    *  ..
c      
c   2    *    *    *  ..
c      
c   1    *    *    *  ..
c                         x ->
c        1    2    3 
c
c
      real*8  node(lx,ly,0:8)
c
c.....help array for temporarely storage of fluid densities
c
      real*8  n_hlp(lx,ly,0:8)
c
c.....average velocity, computed by subroutine 'write_velocity'
c
      real*8  vel
c
c.....startup information message     
c
      write (6,*)
      write (6,*) 'anb 1.0 [ICMMES_1] (2008/06/09)'
      write (6,*) 'Copyright (C) 1998-2008 Joerg Bernsdorf /'
      write (6,*) 'anb comes with ABSOLUTELY NO WARRANTY;' 
      write (6,*) 'This is free software, and you are welcome to redistr
     &ibute it' 
      write (6,*) 'under certain conditions; see the file COPYING for de
     &tails.'
      write (6,*) 'All rights reserved.'
c
      write (6,*)
      write (6,*) '****************************************************'
      write (6,*) '***              anb starting ...                ***'
      write (6,*) '****************************************************'
      write (6,*) '*** Precompiled for lattice size lx = ',lx 
      write (6,*) '***                              ly = ',ly 
      write (6,*) '****************************************************'
      write (6,*) '***'
c
c=======================================================================
c     begin initialisation
c=======================================================================
c
c.....initialize error flag
c
      error = .false.
c
c.....read parameter file
c.....in this file you can enter all relevant parameters,  
c     only lattice size must be fixed befor compilation !
c
      call read_parametrs(error,t_max,density,accel,omega,r_rey)
c
c.....if an I/O error occurs while reading the parameter file, 
c     the "error"-flag is set "true" and the program stops.
c
      if (error) goto 990
c
c.....read obstacle file
c.....in this file you can enter the x-,and y-coordinates of 
c     any obstacles,wall boundaries are also defined here by 
c     adding single obstacles.
c
      call read_obstacles(error,obst,lx,ly)     
c
c.....if an I/O error occurs reading while the obstacle file, 
c     the "error"-flag is set "true" and the program stops.
c
      if (error) goto 990
c
      call init_density(lx,ly,density,node)
c
c.....mean sqaure flow is monitored and stored in the file'anb_qxm.out',
c     one value for each iteration. The file is opened here.
c
      open(10,file='anb_qxm.out')      
c
c=======================================================================
c     end initialisation
c=======================================================================
c=======================================================================
c     begin iterations
c=======================================================================
c
c.....main loop
c
      do 100 time = 1, t_max
c
c.......the integral fluid density is checke each t_max/10 iteration.
c       this is a good indicator, if the program is going to crash ...
c       The integral fluid density should be constant all time ...
c
        if (time .ge. 10 .and. mod(time,t_max/10) .eq. 0) then
c        
          call check_density(lx,ly,node,time)
c
        end if
c
c.......directed flow is induced by density redistribution in the first 
c       lattice column. This is not too clever, since the resulting
c       reynolds number can not be controlled and reaching stady state 
c       takes quite some time, but it is simple and it works ...
c
        call redistribute(lx,ly,obst,node,accel,density)      
c
c.......density propagation: all fluid densities are propagated from
c       non-occupied nodes along the lattice connection lines
c       to their next neighbours, periodic boundary conditiones are 
c       applied in each direction.
c
        call propagate(lx,ly,node,n_hlp)
c
c.......bounc back from obstacles: this is the no-slip boundary-
c       condition.
c       The velocity vector of all fluid densities is inverted, so all
c       the fluid densities will be sent back to the node  where they 
c       were located before the last propagation step, but with opposite
c       velocity vector 
c       ... there exist lots of other possibilities.
c
        call bounceback(lx,ly,obst,node,n_hlp)
c
c.......density relaxation: a single time relaxation with relaxation 
c       parameter omega is applied here. This step is only "local", 
c       nothing is propagated through the lattice.
c
        call relaxation(density,omega,lx,ly,node,n_hlp,obst)
c
c.......average flow velocity is computed at a cross section in the  
c       middle of the channel and written to the file "anb_qx.out"
c       ever iteration. this is a good controll for the convergence of
c       the program.
c
        call write_velocity(lx,ly,time,obst,node,vel)        
c
c.....end of the main loop
c
  100 continue
c
c.....compute fluid-velocities u,v and pressure from velocity 
c     distribution, and write to file anb_rs.out.
c
      call write_results(lx,ly,obst,node,density)
c
c.....compute reynolds number
c
      call comp_rey(lx,ly,obst,node,time,omega,density,r_rey)
c
      goto 999
c
c.....here we get only, if the "error" flag was set "true",  something
c     went wrong reading the files .The program stops with an error
c     message.
c
  990 write (6,*) '!!! error: program stopped during iteration =', time
      write (6,*) '!!!'
c      
  999 continue
c
c.....the file for mean square flow output (anb_qx.out) is closed. 
c     Look at this file before starting any other evaluating, it tells 
c     you, if you reached stady state and if you got reasonable results!
c
      close(10)
c
      write (6,*) '********************    end     ********************'
c
c      
      stop
      end	
c
c
      subroutine read_parametrs(error,t_max,density,accel,omega,r_rey)
c
************************************************************************
*                                                                      *
*     Input run-time parameters from file 'anb.par'                    *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      real*8  density,accel,omega,r_rey
c
      integer  t_max
c
      logical  error
c
c.....open parameter file
c
      open(10,file='anb.par')
c
c.......line 1: number of iterations
c
        read(10,*,err=900) t_max
c
c.......line 2: fluid density per link
c
        read(10,*,err=900) density
c
c.......line 3: density redistribution
c
        read(10,*,err=900) accel
c
c.......line 4: relaxation parameter
c
        read(10,*,err=900) omega
c
c.......line 5: linear dimension (for reynolds number)
c
        read(10,*,err=900) r_rey
c
c.....close parameter file
c
      close(10)      
c
c.....information message
c
      write (6,*) '*** Paramters read from file anb.par.'
      write (6,*) '***'
c
      goto 999
c
c.....error message: file read error
c
  900 write (6,*) '!!! Error reading file anb.par'
      write (6,*) '!!!'
c
      goto 990
c
  990 error = .true.
c
  999 continue
c
c
      return
      end
c
c
      subroutine read_obstacles(error,obst,lx,ly)
c
************************************************************************
*                                                                      *
*     Input obstacle file 'anb.obs'                                    *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1998/08/25                                          *
*                                                                      *
************************************************************************
c
      implicit none
c     
      integer  lx,ly
c
      logical  error,obst(lx,ly)
c
c.....local variables
c
      integer  x,y
c
c.....no obstacles in obstacle array
c
        do 10 y = 1, ly
          do 10 x = 1, lx
c
   10       obst(x,y) = .false.
c
c.....open obstacle file
c.....the obstacle file is defined by the x- and y- coordinates of the
c     obstacles. Each obstacle has a line of its own. Also boundaries
c     have to be defined here.
c
      open(10,file='anb.obs')
c
c.....read obstacle coordinates
c
   20 continue
c
        read(10,*,end=50,err=900) x,y
c
c.......check if obstacle inside domain boundaries
c
        if (x .le. lx .and. y .le. ly) then
c
c.......define obstacle
c
          obst(x,y) = .true.
c
        else
c
          write(6,*) '!!! Obstacle out of range, skipped'
          write(6,*) '!!! lx = ', x, ' , ly = ', y
          write(6,*) '!!!'
c
        end if
c
      goto 20
c
   50 continue
c
c.....close obstacle file
c
      close(10)
c
      write (6,*) '*** Geometry information read from file anb.obs.'
      write (6,*) '***'
c
      goto 999
c
c.....error message: file read error
c
  900 write (6,*) '!!! Error reading file anb.obs'
      write (6,*) ' '
c
      goto 990
c
  990 error = .true.
c
  999 continue
c
c
      return
      end
c
c
      subroutine init_density(lx,ly,density,node)
c
************************************************************************
*                                                                      *
*     Initialize density distribution function n with equilibrium      *
*     for zero velocity                                                *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly
c
      real*8  density,node(lx,ly,0:8)
c
c.....local variables
c
      integer  x,y
      real*8  t_0,t_1,t_2
c
c.....compute weighting factors (depending on lattice geometry)
c
      t_0 = density *  4.d0 / 9.d0
      t_1 = density /  9.d0
      t_2 = density / 36.d0
c
c.....loop over computational domain
c
      do 10 x = 1, lx
        do 10 y = 1, ly
c
c.........zero velocity density
c
          node(x,y,0) = t_0
c
c.........equilibrium densities for axis speeds
c
          node(x,y,1) = t_1
          node(x,y,2) = t_1
          node(x,y,3) = t_1
          node(x,y,4) = t_1
c
c.........equilibrium densities for diagonal speeds
c
          node(x,y,5) = t_2
          node(x,y,6) = t_2
          node(x,y,7) = t_2
          node(x,y,8) = t_2
c
   10 continue
c
c
      return
      end
c
c
      subroutine check_density(lx,ly,node,time)  
c
************************************************************************
*                                                                      *
*     compute integral density                                         *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 2000/01/14                                          *
*                                                                      *
************************************************************************
c
      implicit none
c     
      integer  lx,ly,time
c
      real*8  node(lx,ly,0:8)
c
c
c.....local variables
c
      integer  x,y,n
c
      real*8 n_sum
c
c
      n_sum = 0.d0
c
c.....loop over computational domain
c
        do 10 y = 1, ly
          do 10 x = 1, lx
c
c...........loop over all densities
c
            do 10 n = 0, 8
c
c.............sum up densities
c
   10         n_sum = n_sum + node(x,y,n)
c
      write(6,*) '*** Iteration number = ', time
      write(6,*) '*** Integral density = ', n_sum
      write(6,*) '***'
c
c
      return
      end
c
c
      subroutine redistribute(lx,ly,obst,node,accel,density)
c
************************************************************************
*                                                                      *
*     density redistribution in first lattice column                   *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly
c
      logical  obst(lx,ly)
c
      real*8   node(lx,ly,0:8),accel,density
c
c.....local variables
c
      integer  y
      real*8  t_1,t_2
c
c.....compute weighting factors (depending on lattice geometry) for 
c     increasing/decreasing inlet densities
c
      t_1 = density * accel /  9.d0
      t_2 = density * accel / 36.d0

      do 10 y = 1, ly
c
c.......accelerate flow only on non-occupied nodes
c
        if (.not. obst(1,y) .and.
c
c.........check to avoid negative densities
c
     &    node(1,y,3) - t_1 .gt. 0. .and.
     &    node(1,y,6) - t_2 .gt. 0. .and.
     &    node(1,y,7) - t_2 .gt. 0.) then 
c
c.........increase east
c
          node(1,y,1) = node(1,y,1) + t_1
c
c.........decrease west
c
           node(1,y,3) = node(1,y,3) - t_1
c
c.........increase north-east
c
          node(1,y,5) = node(1,y,5) + t_2
c
c.........decrease north-west
c
          node(1,y,6) = node(1,y,6) - t_2
c
c.........decrease south-west
c
          node(1,y,7) = node(1,y,7) - t_2
c
c.........increase south-east
c
          node(1,y,8) = node(1,y,8) + t_2
c
        end if
c
   10 continue
c
c
      return
      end
c
c
      subroutine bounceback(lx,ly,obst,node,n_hlp)
c
************************************************************************
*                                                                      *
*     Fluid densities are rotated. By the next propagation step, this  *
*     results in a bounce back from obstacle nodes.                    *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly
c
      logical  obst(lx,ly)
c
      real*8   node(lx,ly,0:8),n_hlp(lx,ly,0:8)
c
c.....local variables
c
      integer  x,y
c
c.....loop over all nodes
c
      do 10 x = 1, lx
        do 10 y = 1, ly
c
c.........consider only obstacle nodes
c
          if (obst(x,y)) then
c
c...........rotate all ensities and write back to node
c
c...........east
c
            node(x,y,1) = n_hlp(x,y,3)
c
c...........north
c
            node(x,y,2) = n_hlp(x,y,4)
c
c...........west
c
            node(x,y,3) = n_hlp(x,y,1)
c
c...........south
c
            node(x,y,4) = n_hlp(x,y,2)
c
c...........north-east
c
            node(x,y,5) = n_hlp(x,y,7)
c
c...........north-west
c
            node(x,y,6) = n_hlp(x,y,8)
c
c...........south-west
c
            node(x,y,7) = n_hlp(x,y,5)
c
c...........south-east
c
            node(x,y,8) = n_hlp(x,y,6)
c
          end if
c
   10 continue
c
c
      return
      end
c
c
      subroutine propagate(lx,ly,node,n_hlp)
c
************************************************************************
*                                                                      *
*     Propagate fluid densities to their next neighbour nodes          *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly
c
      real*8   node(lx,ly,0:8),n_hlp(lx,ly,0:8)
c
c.....local variables
c
      integer  x,y,x_e,x_w,y_n,y_s
c
c.....loop over all nodes
c
      do 10 x = 1, lx
        do 10 y = 1, ly
c
c.........compute upper and right next neighbour nodes with regard
c         to periodic boundaries
c
          y_n = mod(y,ly) + 1
          x_e = mod(x,lx) + 1
c
c.........compute lower and left next neighbour nodes with regard to
c         periodic boundaries

          y_s = ly - mod(ly + 1 - y, ly)
          x_w = lx - mod(lx + 1 - x, lx)
c
c.........density propagation
c
c.........zero: just copy
c
          n_hlp(x  ,y, 0 ) = node(x,y,0)
c
c.........east
c
          n_hlp(x_e,y, 1) = node(x,y,1)
c
c.........north
c
          n_hlp(x  ,y_n,2) = node(x,y,2)
c
c.........west
c
          n_hlp(x_w,y,3  ) = node(x,y,3)
c
c.........south
c
          n_hlp(x  ,y_s,4) = node(x,y,4)
c
c.........north-east
c
          n_hlp(x_e,y_n,5) = node(x,y,5)
c
c.........north-west
c
          n_hlp(x_w,y_n,6) = node(x,y,6)
c
c.........south-west
c
          n_hlp(x_w,y_s,7) = node(x,y,7)
c
c.........south-east
c
          n_hlp(x_e,y_s,8) = node(x,y,8)
c
   10 continue
c
c
      return
      end
c
c
      subroutine relaxation(density,omega,lx,ly,node,n_hlp,obst)
c
************************************************************************
*                                                                      *
*     One-step density relaxation process                              *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly
c
      logical  obst(lx,ly)
c
      real*8   density,omega,node(lx,ly,0:8),n_hlp(lx,ly,0:8)
c
c.....local variables
c
      integer  x,y,i
c
      real*8  c_squ,t_0,t_1,t_2,u_x,u_y,u_n(8),n_equ(0:8),u_squ,d_loc
c
c.....NEW LOCAL VARIABLES
c
      real*8 c_1a,c_1b,c_2,cpl_omega
c
c.....weighting factors (depending on lattice geometry)
c
      t_0 = 4.d0 / 9.d0
      t_1 = 1.d0 /  9.d0
      t_2 = 1.d0 / 36.d0
c
c.....complement of omega
c
      cpl_omega = 1.d0 - omega
c
c.....loop over all nodes
c.....attention: actual densities are stored after the propagation
c                step in the help-array n_hlp !
c
c.....LOOP ORDER CHANGED
c
      do 10 y = 1, ly
        do 10 x = 1, lx
c
c.........only free nodes are considered here
c
          if (.not. obst(x,y)) then
c
c...........integral local density
c
c...........UNROLLING BY HAND AND INDEX ORDER CHANGED
c
              d_loc = n_hlp(x,y,0) + n_hlp(x,y,1) + n_hlp(x,y,2)
     &              + n_hlp(x,y,3) + n_hlp(x,y,4) + n_hlp(x,y,5)
     &              + n_hlp(x,y,6) + n_hlp(x,y,7) + n_hlp(x,y,8)
c
c...........x-, and y- velocity components
c
              u_x = (n_hlp(x,y,1) + n_hlp(x,y,5) + n_hlp(x,y,8)
     &           -(n_hlp(x,y,3) + n_hlp(x,y,6) + n_hlp(x,y,7))) / d_loc
c
              u_y = (n_hlp(x,y,2) + n_hlp(x,y,5) + n_hlp(x,y,6)
     &           -(n_hlp(x,y,4) + n_hlp(x,y,7) + n_hlp(x,y,8))) / d_loc
c
c.............square velocity
c
              u_squ = u_x * u_x + u_y * u_y
c
c.............n- velocity compnents (n = lattice node connection vectors)
c.............this is only necessary for clearence, and only 3 speeds would
c.............be necessary
c
              u_n(1) =   u_x
              u_n(2) =         u_y
              u_n(3) = - u_x
              u_n(4) =       - u_y
              u_n(5) =   u_x + u_y
              u_n(6) = - u_x + u_y
              u_n(7) = - u_x - u_y
              u_n(8) =   u_x - u_y
c
c.............NEW VARIABLES TO SIMPLIFY ALGEBRAIC EXPRESSIONS
c
              c_1a = d_loc * t_1
              c_1b = d_loc * t_2
c
              c_2 = 1.5d0 * u_squ
c
c.............equilibrium densities
c.............RATHER OPTIMISED VERSION
c
c.............zero velocity density
c
              n_equ(0) = t_0 * d_loc * (1.d0 -  1.5d0 * u_squ)
c
c.............axis speeds (factor: t_1)
c
            n_equ(1) = c_1a * (1 + u_n(1) * (3.d0 + u_n(1) * 4.5d0)-c_2)
            n_equ(2) = c_1a * (1 + u_n(2) * (3.d0 + u_n(2) * 4.5d0)-c_2)
            n_equ(3) = c_1a * (1 + u_n(3) * (3.d0 + u_n(3) * 4.5d0)-c_2)
            n_equ(4) = c_1a * (1 + u_n(4) * (3.d0 + u_n(4) * 4.5d0)-c_2)
c
c...........diagonal speeds (factor: t_2)
c
            n_equ(5) = c_1b * (1 + u_n(5) * (3.d0 + u_n(5) * 4.5d0)-c_2)
            n_equ(6) = c_1b * (1 + u_n(6) * (3.d0 + u_n(6) * 4.5d0)-c_2)
            n_equ(7) = c_1b * (1 + u_n(7) * (3.d0 + u_n(7) * 4.5d0)-c_2)
            n_equ(8) = c_1b * (1 + u_n(8) * (3.d0 + u_n(8) * 4.5d0)-c_2)
c
c.............relaxation
c.............RATHER OPTIMISED VERSION
c
              node(x,y,0) = n_hlp(x,y,0) 
     &           * cpl_omega + omega * n_equ(0)
              node(x,y,1) = n_hlp(x,y,1) 
     &             * cpl_omega + omega * n_equ(1)
              node(x,y,2) = n_hlp(x,y,2) 
     &             * cpl_omega + omega * n_equ(2)
              node(x,y,3) = n_hlp(x,y,3) 
     &             * cpl_omega + omega * n_equ(3)
              node(x,y,4) = n_hlp(x,y,4) 
     &             * cpl_omega + omega * n_equ(4)
              node(x,y,5) = n_hlp(x,y,5) 
     &             * cpl_omega + omega * n_equ(5)
              node(x,y,6) = n_hlp(x,y,6) 
     &             * cpl_omega + omega * n_equ(6)
              node(x,y,7) = n_hlp(x,y,7) 
     &             * cpl_omega + omega * n_equ(7)
              node(x,y,8) = n_hlp(x,y,8) 
     &             * cpl_omega + omega * n_equ(8)
c
          end if
c
   10 continue
c
c
      return
      end
c
c
      subroutine write_velocity(lx,ly,time,obst,node,vel)        
c
************************************************************************
*                                                                      *
*     compute average velocity                                         *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly,time
c
      logical  obst(lx,ly)
c
      real*8   node(lx,ly,0:8),vel
c
c.....local variables
c
      integer  x,y,i,n_free
c
      real*8  u_x,d_loc
c
c.....loop over channel cross section at half channel length lx/2
c
      x = int(float(lx) / 2.d0)
c
c.....initialize counter
c
      n_free = 0
      u_x = 0.d0
c
      do 10 y = 1, ly
c
c.......only non-occupied nodes are considered here
c
        if(.not. obst(x,y)) then
c
c.........integral local density
c
c.........initialize variable d_loc
c
          d_loc = 0.d0
c
          do 20 i = 0, 8 
c
            d_loc = d_loc + node(x,y,i)
c
   20     continue
c
c.........x-, and y- velocity components
c
          u_x = u_x + (node(x,y,1) + node(x,y,5) + node(x,y,8)
     &               -(node(x,y,3) + node(x,y,6) + node(x,y,7))) / d_loc
c
c.........increase counter
c
          n_free = n_free + 1
c
        end if
c
   10 continue
c
c.....average velocity
c
      vel = u_x / float(n_free)
c
c.....write to file
c
      write(10,*) time, vel
c
c
      return
      end
c
c
      subroutine  write_results(lx,ly,obst,node,density)
c
************************************************************************
*                                                                      *
*     Output of rsults to file 'anb.dat'                               *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c     
      integer  lx,ly
c
      real*8  node(lx,ly,0:8),density
c
      logical  obst(lx,ly)
c
c.....local variables
c
      integer  x,y,i,obsval
c
      real*8  u_x,u_y,d_loc,press,c_squ
c
c.....square speed of sound
c
      c_squ = 1.d0 / 3.d0
c
c.....open results output file
c
      open(11,file='anb.dat')
c
c.....write header for postprocessing with TECPLOT software
c.....uncomment following line, if this header should be printed
c
c
      write(11,*) 'VARIABLES = X, Y, VX, VY, PRESS, OBST' 
      write(11,*) 'ZONE I=', lx, ', J=', ly, ', F=POINT'
c
c.....loop over all nodes
c.....attention: actual densities are stored after the propagation
c                step in the help-array n_hlp !
c
      do 10 y = 1, ly
        do 10 x = 1, lx
c
c.........if obstacle node, nothing is to do ...
c
          if (obst(x,y)) then 
c
c...........obstacle indicator
c
            obsval = 1
c
c...........velocity components = 0
c
            u_x = 0.d0
            u_y = 0.d0
c
c...........pressure = average pressure
c
            press = density * c_squ
c
          else
c
c...........integral local density
c
c...........initialize variable d_loc
c
            d_loc = 0.d0
c
            do 20 i = 0, 8 
c
              d_loc = d_loc + node(x,y,i)
c
   20       continue
c
c...........x-, and y- velocity components
c
            u_x = (node(x,y,1) + node(x,y,5) + node(x,y,8)
     &           -(node(x,y,3) + node(x,y,6) + node(x,y,7))) / d_loc
c
            u_y = (node(x,y,2) + node(x,y,5) + node(x,y,6)
     &           -(node(x,y,4) + node(x,y,7) + node(x,y,8))) / d_loc
c
c...........pressure
c
            press = d_loc * c_squ
c 
            obsval = 0
c 
          end if
c
c.........write results to file
c
          write(11,*) x, y, u_x, u_y, press, obsval
c
   10 continue
c
c.....close file 'anb.dat'
c
      close(11)
c
c
      return
      end
c
c
      subroutine comp_rey(lx,ly,obst,node,time,omega,density,r_rey)
c
************************************************************************
*                                                                      *
*     compute Reynolds number                                          *
*                                                                      *
*             Joerg BERNSDORF                                          *
*             C&C Research Laboratories, NEC Europe Ltd.               *
*             Rathausalle 10                                           *
*             D-53757 Sankt Augustin, Germany                          *
*                                                                      *
*     Last change: 1999/09/28                                          *
*                                                                      *
************************************************************************
c
      implicit none
c
      integer  lx,ly,time
c
      real*8  density,node(lx,ly,0:8),omega,r_rey
c
      logical  obst(lx,ly)
c
c.....local variables
c
      real*8  vel,visc,rey
c
c.....compute average velocity
c
      call write_velocity(lx,ly,time,obst,node,vel)        
c
c.....compute viscosity
c
      visc = 1.d0 / 6.d0 * (2.d0 / omega - 1.d0)
c
c.....compute Reynolds number
c
      rey = vel * r_rey / visc
c
c.....messages
c
       write (6,*) '*** Calculations finished, results:'
       write (6,*) '***'
       write (6,*) '*** viscosity = ', visc
       write (6,*) '*** average velocity = ', vel
       write (6,*) '*** Reynolds number = ', rey
       write (6,*) '***'
       write (6,*) '*** In the file anb.dat, you can find local'
       write (6,*) '*** information about the simulated flow.'
       write (6,*) '***'
       write (6,*) '*** In the file anb_qx.out, you can find the average
     &'
       write (6,*) '*** flow velocity plotted as a function of time.'
       write (6,*) '***'
c
c
      return
      end
