!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 program POP_HaloTest

!----------------------------------------------------------------------
!
!  this program tests POP global reduction operations, including
!  global sums, minval, maxval, etc.  This test code relies on the
!  redistribution mod working correctly for gather/scatter operations.
!
!----------------------------------------------------------------------

   use POP_KindsMod
   use POP_ErrorMod
   use POP_CommMod
   use POP_DomainSizeMod
   use POP_BlocksMod
   use POP_DistributionMod
   use POP_RedistributeMod
   use POP_ReductionsMod
   use POP_FieldMod
   use POP_GridHorzMod
   use POP_HaloMod

   implicit none

!----------------------------------------------------------------------
!
!  local variables
!
!----------------------------------------------------------------------

   integer (POP_i4) :: &
      errorCode,       &! error flag
      istat,           &! allocate error flag
      icount,          &! counter for validation
      nProcsClinic,    &! number of processors in clinic distribution
      nBlocksClinic,   &! number of local blocks in clinic distribution
      blockID,         &! global block id
      ib,ie,jb,je,     &! beg,end indices for physical domain on block
      ig, jg,          &! global indices for array
      i,j,k,l,n         ! loop indices

   type (POP_block) :: &
      thisBlock         ! block info for local block

   type (POP_distrb) :: &
      distrbClinic       ! baroclinic distribution

   type (POP_halo) :: &
      haloDipole       ! halo data type for dipole updates

   integer (POP_i4), dimension(:), allocatable :: &
      workPerBlock  ! test values

   integer (POP_i4), dimension(:), pointer :: &
      iGlobal       ! global index for this block

   integer (POP_i4), dimension(:), pointer :: &
      jGlobal       ! global index for this block

   integer (POP_i4), &
      dimension(POP_nxBlock, POP_nyBlock, POP_maxBlocksClinic) :: &
      i4Test, i4Expect

   real (POP_r4), &
      dimension(POP_nxBlock, POP_nyBlock, POP_maxBlocksClinic) :: &
      r4Test, r4Expect  ! test values

   real (POP_r8), &
      dimension(POP_nxBlock, POP_nyBlock, POP_maxBlocksClinic) :: &
      r8Test, r8Expect  ! test values

   logical (POP_logical), &
      dimension(POP_nxBlock, POP_nyBlock, POP_maxBlocksClinic) :: &
      mask

   integer (POP_i4), dimension(POP_nxBlock, POP_nyBlock,       &
                               POP_km, POP_maxBlocksClinic) :: &
      i4Test3D

   real (POP_r4), dimension(POP_nxBlock, POP_nyBlock,       &
                            POP_km, POP_maxBlocksClinic) :: &
      r4Test3D

   real (POP_r8), dimension(POP_nxBlock, POP_nyBlock,       &
                            POP_km, POP_maxBlocksClinic) :: &
      r8Test3D

   integer (POP_i4), dimension(POP_nxBlock, POP_nyBlock,               &
                               POP_km, POP_nt, POP_maxBlocksClinic) :: &
      i4Test4D

   real (POP_r4), dimension(POP_nxBlock, POP_nyBlock,               &
                            POP_km, POP_nt, POP_maxBlocksClinic) :: &
      r4Test4D

   real (POP_r8), dimension(POP_nxBlock, POP_nyBlock,               &
                            POP_km, POP_nt, POP_maxBlocksClinic) :: &
      r8Test4D

   integer (POP_i4), dimension(POP_nxGlobal, POP_nyGlobal) :: &
      i4ArrayG

!----------------------------------------------------------------------
!
!  initialize communcation environment
!
!----------------------------------------------------------------------

   call POP_CommInitMessageEnvironment
   call POP_CommInit

   errorCode = POP_Success

!----------------------------------------------------------------------
!
!  create block decomposition and distributions
!
!----------------------------------------------------------------------

   !*** create block decomposition
 
   call POP_BlocksCreate(POP_nxGlobal, POP_nyGlobal, &
                         'cyclic', 'closed', errorCode)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error creating block decomposition')
      call POP_ErrorPrint(errorCode)
      stop
   endif

   !*** initialize a global array with land in the middle

   do j=1,POP_nyGlobal
   do i=1,POP_nxGlobal
      if (i > POP_nxGlobal/2 + POP_nxBlock .or. &
          i < POP_nxGlobal/2 - POP_nxBlock .or. & 
          j > POP_nyGlobal/2 + POP_nyBlock .or. &
          j < POP_nyGlobal/2 - POP_nyBlock) then
         i4ArrayG(i,j) = (i+j)
      else
         i4ArrayG(i,j) = 0_POP_i4
      endif
   end do
   end do

   !*** create artificial work per block with varying work and
   !*** at least one empty block 


   allocate(workPerBlock(POP_numBlocks), stat=istat)
   if (istat > 0) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error allocating work per block')
      call POP_ErrorPrint(errorCode)
      stop
   endif

   do blockID=1,POP_numBlocks
      call POP_BlocksGetBlockInfo(blockID, errorCode,         &
                                  ib=ib, ie=ie, jb=jb, je=je, &
                                  iGlobal=iGlobal, jGlobal=jGlobal)

      if (errorCode /= POP_Success) then
         call POP_ErrorSet(errorCode, &
            'HaloTest: error getting block info')
         call POP_ErrorPrint(errorCode)
         stop
      endif

      workPerBlock(blockID) = 0
      do j=jb,je
      do i=ib,ie
         ig = iGlobal(i)
         jg = jGlobal(j)
         if (ig > 0 .and. jg > 0) then
            if (i4ArrayG(iGlobal(i),jGlobal(j)) /= 0) then
               workPerBlock(blockID) = workPerBlock(blockID) + 1
            endif
         endif
      end do
      end do
   end do

   !*** create baroclinic distribution

   nprocsClinic = POP_CommGetNumProcs(POP_Communicator)

   distrbClinic = POP_DistributionCreate(POP_distribMethodRake, &
                                         nProcsClinic, &
                                         workPerBlock, errorCode)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error creating clinic distribution')
      call POP_ErrorPrint(errorCode)
      stop
   endif

!----------------------------------------------------------------------
!
!  initialize halo updates
!
!----------------------------------------------------------------------

   haloDipole = POP_HaloCreate(distrbClinic, 'closed', 'cyclic', &
                               POP_nxGlobal, errorCode)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error creating dipole halo')
      call POP_ErrorPrint(errorCode)
      stop
   endif

!----------------------------------------------------------------------
!
!  initialize test arrays by scattering global arrays to distribution
!
!----------------------------------------------------------------------

   i4Test = -999
   r4Test = -999.0_POP_r4
   r8Test = -999.0_POP_r8

   call POP_RedistributeScatter(i4Test, i4ArrayG, &
                              POP_masterTask, distrbClinic, errorCode) 

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error in i4 scatter')
   endif

   mask = (i4Test /= 0)
   r4Test = i4Test*10.0_POP_r4
   r8Test = i4Test*100.0_POP_r8

!----------------------------------------------------------------------
!
!  set up expected values for dipole case
!
!----------------------------------------------------------------------

   call POP_DistributionGet(distrbClinic, errorCode, &
                            numLocalBlocks=nBlocksClinic)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'HaloTest: error getting local block count')
      call POP_ErrorPrint(errorCode)
      stop
   endif

   do n=1,nBlocksClinic

      call POP_DistributionGetBlockID(distrbClinic, n, blockID, &
                                      errorCode)

      if (errorCode /= POP_Success) then
         call POP_ErrorSet(errorCode, &
            'HaloTest: error getting block ID')
         call POP_ErrorPrint(errorCode)
         stop
      endif

      thisBlock = POP_BlocksGetBlock(blockID, errorCode)

      if (errorCode /= POP_Success) then
         call POP_ErrorSet(errorCode, &
            'HaloTest: error getting block')
         call POP_ErrorPrint(errorCode)
         stop
      endif

      do j=1,POP_nyBlock
      do i=1,POP_nxBlock
         ig = thisBlock%iGlobal(i)
         jg = thisBlock%jGlobal(j)

         if (ig <= 0 .or. jg <= 0) then ! outside domain
            i4Expect(i,j,n) = 0_POP_i4
         else
            i4Expect(i,j,n) = i4ArrayG(ig,jg)
         endif
      end do
      end do
   end do

   r4Expect = i4Expect*10.0_POP_r4
   r8Expect = i4Expect*100.0_POP_r8

   do n=1,nBlocksClinic
   do k=1,POP_km
      i4Test3D(:,:,k,n) = i4Test(:,:,n)*k
      r4Test3D(:,:,k,n) = r4Test(:,:,n)*k
      r8Test3D(:,:,k,n) = r8Test(:,:,n)*k
   end do
   end do

   do n=1,nBlocksClinic
   do l=1,POP_nt
   do k=1,POP_km
      i4Test4D(:,:,k,l,n) = i4Test(:,:,n)*k*l
      r4Test4D(:,:,k,l,n) = r4Test(:,:,n)*k*l
      r8Test4D(:,:,k,l,n) = r8Test(:,:,n)*k*l
   end do
   end do
   end do

!----------------------------------------------------------------------
!
!  test 2d halo dipole updates
!
!----------------------------------------------------------------------

   call POP_HaloUpdate(i4Test, haloDipole,             &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0_POP_i4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 2d i4 dipole halo')
   endif

   if (count(i4Test /= i4Expect) > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 2d i4 dipole halo')
   endif

   call POP_HaloUpdate(r4Test, haloDipole,             &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 2d r4 dipole halo')
   endif

   if (count(r4Test /= r4Expect) > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 2d r4 dipole halo')
   endif

   call POP_HaloUpdate(r8Test, haloDipole,             &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r8)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 2d r8 dipole halo')
   endif

   if (count(r8Test /= r8Expect) > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 2d r8 dipole halo')
   endif

!----------------------------------------------------------------------
!
!  test 3d halo dipole updates
!
!----------------------------------------------------------------------

   call POP_HaloUpdate(i4Test3D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0_POP_i4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 3d i4 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do k=1,POP_km
      icount = icount + count(i4Test3D(:,:,k,n) /= i4Expect(:,:,n)*k)
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 3d i4 dipole halo')
   endif

   call POP_HaloUpdate(r4Test3D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 3d r4 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do k=1,POP_km
      icount = icount + count(r4Test3D(:,:,k,n) /= r4Expect(:,:,n)*k)
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 3d r4 dipole halo')
   endif

   call POP_HaloUpdate(r8Test3D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r8)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 3d r8 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do k=1,POP_km
      icount = icount + count(r8Test3D(:,:,k,n) /= r8Expect(:,:,n)*k)
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 3d r8 dipole halo')
   endif

!----------------------------------------------------------------------
!
!  test 4d halo dipole updates
!
!----------------------------------------------------------------------

   call POP_HaloUpdate(i4Test4D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0_POP_i4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 4d i4 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do l=1,POP_nt
   do k=1,POP_km
      icount = icount + count(i4Test4D(:,:,k,l,n) /= i4Expect(:,:,n)*k*l)
   end do
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 4d i4 dipole halo')
   endif

   call POP_HaloUpdate(r4Test4D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r4)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 4d r4 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do l=1,POP_nt
   do k=1,POP_km
      icount = icount + count(r4Test4D(:,:,k,l,n) /= r4Expect(:,:,n)*k*l)
   end do
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 4d r4 dipole halo')
   endif

   call POP_HaloUpdate(r8Test4D, haloDipole,           &
                       POP_gridHorzLocCenter,          &
                       POP_fieldKindScalar, errorCode, &
                       fillValue=0.0_POP_r8)

   if (errorCode /= POP_Success) then
      call POP_ErrorSet(errorCode, &
         'Halo test: error in 4d r8 dipole halo')
   endif

   icount = 0
   do n=1,nBlocksClinic
   do l=1,POP_nt
   do k=1,POP_km
      icount = icount + count(r8Test4D(:,:,k,l,n) /= r8Expect(:,:,n)*k*l)
   end do
   end do
   end do

   if (icount > 0) then
      call POP_ErrorSet(errorCode, &
         'Halo test: bad results from 4d r8 dipole halo')
   endif

!----------------------------------------------------------------------
!
!  clean up
!
!----------------------------------------------------------------------

   call POP_HaloDestroy(haloDipole, errorCode)

   call POP_ErrorPrint(errorCode)
   call POP_CommExitMessageEnvironment

!----------------------------------------------------------------------

 end program POP_HaloTest

!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||