      program main
      implicit none
      include 'cmax.h'

      integer s1,s2,s3,s4, max_distance
      parameter (s1=4,s2=5,s3=6,s4=7, max_distance=12)

      integer axis, distance
      integer s_i(s1,s2,s3,s4),d1_i(s1,s2,s3,s4),d2_i(s1,s2,s3,s4)
      real*8 s_d(s1,s2,s3,s4),d1_d(s1,s2,s3,s4),d2_d(s1,s2,s3,s4)
      logical arrays_equal_i, arrays_equal_d
      
      integer i,j,k,l


      do l = 1, s4
         do k = 1,s3
            do j = 1, s2
               do i = 1, s1
                  s_i(i,j,k,l) = i + j*10 + k*100 + l*1000
                  s_d(i,j,k,l) = i + j*10 + k*100 + l*1000
               end do
            end do
         end do
      end do

      do axis = 1,4
         do distance = 0-max_distance, max_distance
c            print *, "axis = ",axis," distance = ",distance

            call cmax_cshift_4d(d2_i, s_i, axis, distance,cmax_integer,
     &                          s1,s2,s3,s4)
            call test_cshift_i(d1_i, axis, distance, s_i, s1,s2,s3,s4)
            if (.not. arrays_equal_i(d1_i,d2_i, s1,s2,s3,s4)) then
               print *, "**** failed for integers ****"
            end if

            call cmax_cshift_4d(d2_d, s_d, axis, distance,cmax_double,
     &                          s1,s2,s3,s4)
            call test_cshift_d(d1_d, axis, distance, s_d, s1,s2,s3,s4)
            if (.not. arrays_equal_d(d1_d,d2_d, s1,s2,s3,s4)) then
               print *, "**** failed for doubles ****"
            end if

         end do
      end do
      end


      logical function arrays_equal_i(arr1,arr2,s1,s2,s3,s4)
      integer s1,s2,s3,s4
      integer arr1(s1,s2,s3,s4), arr2(s1,s2,s3,s4)
      
      integer i,j,k,l
      logical result
      
      result = .true.
      
      do l = 1, s4
         do k = 1,s3
            do j = 1, s2
               do i = 1, s1
                  if (arr1(s1,s2,s3,s4) .ne. arr2(s1,s2,s3,s4)) then
                     result = .false.
                     end if
               end do
            end do
         end do
      end do
      
      arrays_equal_i = result
      
      return
      end



      subroutine test_cshift_i(dest, axis, distance, source,s1,s2,s3,s4)
      integer axis, distance, s1,s2,s3,s4
      integer dest(s1,s2,s3,s4), source(s1,s2,s3,s4)
      
      integer i,j,k,l, r_dist
      
      if (axis .eq. 1) then
         r_dist = mod(s1 + mod(distance, s1), s1)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(mod(i + r_dist - 1, s1) + 1, j, k,l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 2) then
         
         r_dist = mod(s2 + mod(distance, s2), s2)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i, mod(j + r_dist - 1, s2) + 1, k,l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 3) then
         
         r_dist = mod(s3 + mod(distance, s3), s3)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i,j, mod(k + r_dist - 1, s3) + 1, l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 4) then
         
         r_dist = mod(s4 + mod(distance, s4), s4)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i,j,k, mod(l + r_dist - 1, s4) + 1)
                  end do
               end do
            end do
         end do

      else

         pause "bad axis"

      endif

      return
      end
      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      logical function arrays_equal_d(arr1,arr2,s1,s2,s3,s4)
      integer s1,s2,s3,s4
      real*8 arr1(s1,s2,s3,s4), arr2(s1,s2,s3,s4)
      
      integer i,j,k,l
      logical result
      
      result = .true.
      
      do l = 1, s4
         do k = 1,s3
            do j = 1, s2
               do i = 1, s1
                  if (arr1(s1,s2,s3,s4) .ne. arr2(s1,s2,s3,s4)) then
                     result = .false.
                     end if
               end do
            end do
         end do
      end do
      
      arrays_equal_d = result
      
      return
      end



      subroutine test_cshift_d(dest, axis, distance, source,s1,s2,s3,s4)
      integer axis, distance, s1,s2,s3,s4
      real*8 dest(s1,s2,s3,s4), source(s1,s2,s3,s4)
      
      integer i,j,k,l, r_dist
      
      if (axis .eq. 1) then
         r_dist = mod(s1 + mod(distance, s1), s1)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(mod(i + r_dist - 1, s1) + 1, j, k,l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 2) then
         
         r_dist = mod(s2 + mod(distance, s2), s2)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i, mod(j + r_dist - 1, s2) + 1, k,l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 3) then
         
         r_dist = mod(s3 + mod(distance, s3), s3)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i,j, mod(k + r_dist - 1, s3) + 1, l)
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 4) then
         
         r_dist = mod(s4 + mod(distance, s4), s4)
         do l = 1, s4
            do k = 1,s3
               do j = 1, s2
                  do i = 1, s1
                     dest(i,j,k,l) = 
     &                    source(i,j,k, mod(l + r_dist - 1, s4) + 1)
                  end do
               end do
            end do
         end do

      else

         pause "bad axis"

      endif

      return
      end
      
