SUBROUTINE acoustic3d_f90_setup(P1, P2, P3, c, N)
  INTEGER, INTENT( IN ) :: N
  REAL, DIMENSION (N,N,N), INTENT( INOUT ) :: P1, P2, P3, c

  INTEGER cavityLeft, cavityRight, cavityFront, cavityBack, cavityTop, &
          cavityBottom, cavityTop2, cavityBottom2
  REAL ci, cj, ck, s2

! Set the velocity field

  c(1:N/2,:,:) = 0.05
  c(N/2+1:N,:,:) = 0.3;

  cavityLeft   = 3 * N / 7.0
  cavityRight  = 4 * N / 7.0
  cavityFront  = 3 * N / 7.0
  cavityBack   = 4 * N / 7.0
  cavityTop    = 5 * N / 7.0
  cavityBottom = 6 * N / 7.0
  cavityTop2   = 1 * N / 7.0
  cavityBottom2 = 2 * N / 7.0

  c(cavityTop:cavityBottom,cavityLeft:cavityRight, &
    cavityFront:cavityBack) = 0.02;
  c(cavityTop2:cavityBottom2,cavityLeft:cavityRight, &
    cavityFront:cavityBack) = 0.001;

! Initial pressure distribution
  P1 = 0.0
  P3 = 0.0

  ci = N/2.0
  cj = N/2.0
  ck = N/2.0
  s2 = 64.0 * 9.0 / ((N/2.0) ** 2)

  DO k=1,N
    DO j=1,N
      DO i=1,N
        P2(i,j,k) = exp(-((i-ci)**2+(j-cj)**2+(k-ck)**2) * s2)
      END DO
    END DO
  END DO

  RETURN
END