! Tuned Fortran 77 version

      SUBROUTINE acoustic3d_f77Tuned(N, niters, check)
      INTEGER N, niters, iter
      REAL check
      REAL P1(N,N,N), P2(N,N,N), P3(N,N,N), C(N,N,N)

      CALL acoust3d_f77Tuned_setup(P1,P2,P3,C,N)

      DO iter = 1, niters, 3
        CALL acoustic3d_f77Tuned_stencil(P1,P2,P3,C,N)
        CALL acoustic3d_f77Tuned_stencil(P2,P3,P1,C,N)
        CALL acoustic3d_f77Tuned_stencil(P3,P1,P2,C,N)
      END DO

      check = P1(N/2,N/2,N/2)

      RETURN
      END



      SUBROUTINE acoustic3d_f77Tuned_stencil(P1,P2,P3,C,N)
      INTEGER N
      REAL P1(N,N,N), P2(N,N,N), P3(N,N,N), C(N,N,N)
      INTEGER bi,bj,bk,ni,nj,nk,blockSize

      blockSize = 8

      DO bk=2,N-1,blockSize
        nk = min(bk+blockSize-1,N-1)
        DO bj=2,N-1,blockSize
          nj = min(bj+blockSize-1,N-1)
          DO bi=2,N-1,blockSize
            ni = min(bi+blockSize-1,N-1)
            DO k=bk,nk
              DO j=bj,nj
                DO i=bi,ni
                  P3(i,j,k) = (2-6*C(i,j,k))*P2(i,j,k) 
     .              + C(i,j,k) * (P2(i,j-1,k) + P2(i,j+1,k) 
     .              + P2(i-1,j,k) + P2(i+1,j,k) + P2(i,j,k-1)
     .              + P2(i,j,k+1)) - P1(i,j,k)
                END DO
              END DO
            END DO
          END DO
        END DO
      END DO

      RETURN
      END