d3q81_cxDir Function

  pure function d3q81_cxDir() result(cxdir)
    ! -------------------------------------------------------------------- !
    integer :: cxDir(3,81)
    ! -------------------------------------------------------------------- !
    integer :: iX, iY, iZ
    integer :: dirPos
    integer :: length
    ! -------------------------------------------------------------------- !

    ! save the cartesian directions of the discrete velocities
    dirPos = 0

    do iX=-2,2
      do iY=-2,2
        do iZ=-2,2

          ! calculate the square of the length
          length = iX*iX + iY*iY + iZ*iZ

          ! if the actual length is smaller than 2.5
          ! this is the same as length^2 .lt. 8 for integer stencil directions
          if ((length > 0) .and. (length < 8)) then
            dirPos = dirPos + 1
            cxDir(:, dirPos) = [ iX, iY, iZ ]
          end if

        end do
      end do
    end do

    ! add the rest direction at the last position
    cxDir(:, 81) = [0, 0, 0]

  end function d3q81_cxDir