      program main
      implicit none
      include 'cmax.h'

      integer dsize, s1,s2,s3,s4
      parameter (dsize=10,s1=4,s2=5,s3=6,s4=7)

      integer s_i(s1,s2,s3,s4)
      integer d1_i_1(dsize,s2,s3,s4),d2_i_1(dsize,s2,s3,s4)
      integer d1_i_2(s1,dsize,s3,s4),d2_i_2(s1,dsize,s3,s4)
      integer d1_i_3(s1,s2,dsize,s4),d2_i_3(s1,s2,dsize,s4)
      integer d1_i_4(s1,s2,s3,dsize),d2_i_4(s1,s2,s3,dsize)
      real*8 s_d(s1,s2,s3,s4)
      real*8 d1_d_1(dsize,s2,s3,s4),d2_d_1(dsize,s2,s3,s4)
      real*8 d1_d_2(s1,dsize,s3,s4),d2_d_2(s1,dsize,s3,s4)
      real*8 d1_d_3(s1,s2,dsize,s4),d2_d_3(s1,s2,dsize,s4)
      real*8 d1_d_4(s1,s2,s3,dsize),d2_d_4(s1,s2,s3,dsize)
      logical m(s1,s2,s3,s4)
      integer pad_i
      real*8 pad_d

      logical arrays_equal_i, arrays_equal_d
      
      integer i,j,k,l

      pad_i = -1
      pad_d = -1

      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

                  m(i,j,k,l) = (i .le. j) .and.
     &                         (j .le. k) .and.
     &                         (k .le. l)
               end do
            end do
         end do
      end do

      call cmax_pack_4d(d2_i_1, s_i, m, 1, pad_i,cmax_integer,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_i(d1_i_1, s_i, m, 1, pad_i,
     &                   dsize,s2,s3,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_i(d1_i_1,d2_i_1, dsize,s2,s3,s4)) then
         print *, "**** failed for integers  1 ****"
      end if

      call cmax_pack_4d(d2_i_2, s_i, m, 2, pad_i,cmax_integer,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_i(d1_i_2, s_i, m, 2, pad_i,
     &                   s1,dsize,s3,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_i(d1_i_2,d2_i_2, s1,dsize,s3,s4)) then
         print *, "**** failed for integers  2 ****"
      end if

      call cmax_pack_4d(d2_i_3, s_i, m, 3, pad_i,cmax_integer,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_i(d1_i_3, s_i, m, 3, pad_i,
     &                   s1,s2,dsize,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_i(d1_i_3,d2_i_3, s1,s2,dsize,s4)) then
         print *, "**** failed for integers  3 ****"
      end if

      call cmax_pack_4d(d2_i_4, s_i, m, 4, pad_i,cmax_integer,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_i(d1_i_4, s_i, m, 4, pad_i,
     &                   s1,s2,s3,dsize, s1,s2,s3,s4)
      if (.not. arrays_equal_i(d1_i_4,d2_i_4, s1,s2,s3,dsize)) then
         print *, "**** failed for integers  4 ****"
      end if

c now for doubles

      call cmax_pack_4d(d2_d_1, s_d, m, 1, pad_d,cmax_double,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_d(d1_d_1, s_d, m, 1, pad_d,
     &                   dsize,s2,s3,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_d(d1_d_1,d2_d_1, dsize,s2,s3,s4)) then
         print *, "**** failed for doubles  1 ****"
      end if

      call cmax_pack_4d(d2_d_2, s_d, m, 2, pad_d,cmax_double,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_d(d1_d_2, s_d, m, 2, pad_d,
     &                   s1,dsize,s3,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_d(d1_d_2,d2_d_2, s1,dsize,s3,s4)) then
         print *, "**** failed for doubles  2 ****"
      end if

      call cmax_pack_4d(d2_d_3, s_d, m, 3, pad_d,cmax_double,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_d(d1_d_3, s_d, m, 3, pad_d,
     &                   s1,s2,dsize,s4, s1,s2,s3,s4)
      if (.not. arrays_equal_d(d1_d_3,d2_d_3, s1,s2,dsize,s4)) then
         print *, "**** failed for doubles  3 ****"
      end if

      call cmax_pack_4d(d2_d_4, s_d, m, 4, pad_d,cmax_double,
     &                  dsize, s1,s2,s3,s4)
      call test_pack_d(d1_d_4, s_d, m, 4, pad_d,
     &                   s1,s2,s3,dsize, s1,s2,s3,s4)
      if (.not. arrays_equal_d(d1_d_4,d2_d_4, s1,s2,s3,dsize)) then
         print *, "**** failed for doubles  4 ****"
      end if

      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_pack_i(dest, source, mask, axis, pad, 
     &                       d1,d2,d3,d4,s1,s2,s3,s4)
      integer axis, d1,d2,d3,d4, s1,s2,s3,s4
      integer dest(d1,d2,d3,d4), source(s1,s2,s3,s4)
      logical mask(s1,s2,s3,s4)
      integer pad
      
      integer i,j,k,l, count
      
      if (axis .eq. 1) then

         do l=1,s4
            do k=1,s3
               do j=1,s2
                  count = 1
                  do i=1,s1
                     if (mask(i,j,k,l)) then
                        if (count .gt. d1) then
                           pause
                        end if
                        dest(count, j,k,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do i=count,d1
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 2) then
         
         do l=1,s4
            do k=1,s3
               do i=1,s1
                  count = 1
                  do j=1,s2
                     if (mask(i,j,k,l)) then
                        if (count .gt. d2) then
                           pause
                        end if
                        dest(i, count, k,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do j=count,d2
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 3) then
         
         do l=1,s4
            do j=1,s2
               do i=1,s1
                  count = 1
                  do k=1,s3
                     if (mask(i,j,k,l)) then
                        if (count .gt. d3) then
                           pause
                        end if
                        dest(i,j, count,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do k=count,d3
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 4) then
         
         do k=1,s3
            do j=1,s2
               do i=1,s1
                  count = 1
                  do l=1,s4
                     if (mask(i,j,k,l)) then
                        if (count .gt. d4) then
                           pause
                        end if
                        dest(i,j,k, count) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do l=count,d4
                     dest(i,j,k,l) = pad
                  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_pack_d(dest, source, mask, axis, pad, 
     &                       d1,d2,d3,d4,s1,s2,s3,s4)
      integer axis, d1,d2,d3,d4, s1,s2,s3,s4
      real*8 dest(d1,d2,d3,d4), source(s1,s2,s3,s4)
      logical mask(s1,s2,s3,s4)
      real*8 pad
      
      integer i,j,k,l, count
      
      if (axis .eq. 1) then

         do l=1,s4
            do k=1,s3
               do j=1,s2
                  count = 1
                  do i=1,s1
                     if (mask(i,j,k,l)) then
                        if (count .gt. d1) then
                           pause
                        end if
                        dest(count, j,k,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do i=count,d1
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 2) then
         
         do l=1,s4
            do k=1,s3
               do i=1,s1
                  count = 1
                  do j=1,s2
                     if (mask(i,j,k,l)) then
                        if (count .gt. d2) then
                           pause
                        end if
                        dest(i, count, k,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do j=count,d2
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 3) then
         
         do l=1,s4
            do j=1,s2
               do i=1,s1
                  count = 1
                  do k=1,s3
                     if (mask(i,j,k,l)) then
                        if (count .gt. d3) then
                           pause
                        end if
                        dest(i,j, count,l) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do k=count,d3
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      elseif (axis .eq. 4) then
         
         do k=1,s3
            do j=1,s2
               do i=1,s1
                  count = 1
                  do l=1,s4
                     if (mask(i,j,k,l)) then
                        if (count .gt. d4) then
                           pause
                        end if
                        dest(i,j,k, count) = source(i,j,k,l)
                        count = count + 1
                     end if
                  end do
                  do l=count,d4
                     dest(i,j,k,l) = pad
                  end do
               end do
            end do
         end do
         
      else

         pause "bad axis"

      endif

      return
      end
