! F90 Tuned version

SUBROUTINE acoustic3d_f90Tuned(N, niters, check)
  INTEGER, INTENT( IN ) :: N, niters
  REAL, INTENT( OUT ) :: check
  REAL, DIMENSION (N,N,N) :: P1, P2, P3, C
  INTEGER iter

  CALL acoustic3d_f90Tuned_setup(P1, P2, P3, c, N)

  IF (MODULO(niters,3) > 0) THEN
    PRINT *, 'Warning: In acoustic3d_f90Tuned: niters mod 3 != 0'
    PRINT *, 'Will do fewer iterations, benchmark result will be off'
  ENDIF

  DO iter=1, niters, 3
    CALL acoustic3d_f90Tuned_stencil(P1, P2, P3, C, N)
    CALL acoustic3d_f90Tuned_stencil(P2, P3, P1, C, N)
    CALL acoustic3d_f90Tuned_stencil(P3, P1, P2, C, N)
  END DO
    
  check = P1(N/2,N/2,N/2)

  RETURN
END

SUBROUTINE acoustic3d_f90Tuned_stencil(P1, P2, P3, C, N)
  INTEGER, INTENT( IN ) :: N
  REAL, DIMENSION (N,N,N), INTENT( IN ) :: P1, P2, C
  REAL, DIMENSION (N,N,N), INTENT( OUT ) :: P3

  P3(2:N-1,2:N-1,2:N-1) = (2-6*C(2:N-1,2:N-1,2:N-1)) &
        * P2(2:N-1,2:N-1,2:N-1) + C(2:N-1,2:N-1,2:N-1) &
        * (P2(1:N-2,2:N-1,2:N-1) + P2(3:N,2:N-1,2:N-1) &
        + P2(2:N-1,1:N-2,2:N-1) + P2(2:N-1,3:N,2:N-1) &
        + P2(2:N-1,2:N-1,1:N-2)+P2(2:N-1,2:N-1,3:N)) &
        - P1(2:N-1,2:N-1,2:N-1)
  RETURN
END