module rnd
  implicit none
contains
  subroutine fill(a,n)
    integer, intent(out), dimension(:,:) :: a
    integer, intent(in) :: n
    real, dimension(size(a,1),size(a,2)) :: r
    call random_number(r)
    a = int(2*n*r-n)
  end subroutine fill
end module rnd
program main
  use rnd
  implicit none
  integer, parameter :: n1=400, n2=400, n3=400
  integer, dimension(n1,n2,n3) :: a, b,c
  integer :: s1, s2, s3
  integer :: dim
  integer, dimension(:,:), allocatable :: sh1, sh2, sh3
  integer, dimension(:), allocatable :: sh_shift
  integer :: sh, rsh
  integer :: n
  integer :: i,j,k,v
  real :: t1, t2

  v = 1
  do k=1,n3
     do j=1,n2
        do i=1,n1
           a(i,j,k) = v
           v = v + 1
        end do
     end do
  end do

  allocate(sh1(n2,n3))
  allocate(sh2(n1,n3))
  allocate(sh3(n1,n2))

  ! sh1 = reshape([(i,i=1,size(sh1))],shape(sh1))
  ! sh2 = reshape([(i,i=1,size(sh2))],shape(sh2))
  ! sh3 = reshape([(i,i=1,size(sh3))],shape(sh3))
  call fill(sh1,10)
  call fill(sh2,10)
  call fill(sh3,10)
!  sh1 = 1
!  sh2 = 1
!  sh3 = 1

  call cpu_time(t1)
  b = cshift(a,sh1,1)
  call cpu_time(t2)
  print *,"cpu time cshift dim=1 ", t2-t1

  n = size(a,1)
  c = 0
  call cpu_time(t1)
  do s2=1,n2
     do s3=1,n3
        sh = modulo(sh1(s2,s3), n)
        rsh = n - sh
        do i=1,rsh
           c(i,s2,s3) = a(i+sh,s2,s3)
        end do
        do i=rsh+1,n
           c(i,s2,s3) = a(i-rsh,s2,s3)
        end do
        c(:,s2,s3) = cshift(a(:,s2,s3),sh1(s2,s3))
     end do
  end do
  call cpu_time(t2)
  print *,"cpu time do loop dim=1 ", t2-t1
  if (any(b /= c)) call abort

  n = size(a,2)
  call cpu_time(t1)
  b = cshift(a,sh2,2)
  call cpu_time(t2)
  print *,"cpu time cshift dim=2 ", t2-t1

  write (10,*) sum(b)
  c = 0
  call cpu_time(t1)
  do s3=1,n3
     do s1=1,n1
        sh = modulo(sh2(s1,s3),n)
        rsh = n - sh
        do i=1,rsh
           c(s1,i,s3) = a(s1,i+sh,s3)
        end do
        do i=rsh+1,n
           c(s1,i,s3) = a(s1,i-rsh,s3)
        end do
     end do
  end do
  call cpu_time(t2)
  if (any(b /= c)) call abort
  print *,"cpu time do loop dim=2 ", t2-t1

  n = size(a,3)
  call cpu_time(t1)
  b = cshift(a,sh3,3)
  call cpu_time(t2)
  print *,"cpu time cshift dim=3 ", t2-t1

  c = 0
  call cpu_time(t1)
  do s2=1,n2
     do s1=1,n1
        sh = modulo(sh3(s1,s2),n)
        rsh = n - sh
        do i=1,rsh
           c(s1,s2,i) = a(s1,s2,i+sh)
        end do
        do i=rsh+1,n
           c(s1,s2,i) = a(s1,s2,i-rsh)
        end do
     end do
  end do
  call cpu_time(t2)
  print *,"cpu time do loop dim=3 ", t2-t1
  if (any(b /= c)) call abort
    
!9000 format (99(3(I3,1X),2X))
end program main